LegPy_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "z3LOpluQhNkM"
   },
   "source": [
    "# LegPy tutorial\n",
    "\n",
    "#### Fernando Arqueros, Jaime Rosado, Victor Moya\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GBwnegUHapmU"
   },
   "source": [
    "LegPy (Low energy gamma-ray simulation in Python) is a Monte Carlo (MC) simulation code for the transportation of gamma rays and electrons with energies up to few MeVs (depending of the medium).\n",
    "\n",
    "Several geometries are supported (cylinder, orthohedron and sphere). Photons or electrons are produced as paralell (or divergent) beams or from isotropic sources with energies following several spectral distributions (monoenergetic, flat, exponential, etc.) including those of realistic X-ray and gamma ray sources. Particles can be transported either through a single material or in geometries filled with two different media.\n",
    "\n",
    "The description of the MC algorithm and the validation of the various approximations are described in https://doi.org/10.1016/j.radmeas.2023.107029"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "E9YHsRxa_q9v"
   },
   "source": [
    "This **tutorial** explains the procedure to run a simulation with the LegPy package. It will be shown that only five code lines in this notebook are sufficient. Later a number of examples will be presented as well as procedures to visualize and analyze the results of the simulation.         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "maXNyDYLbqUo"
   },
   "source": [
    "## Previous steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "y5mmA9wL_q9v"
   },
   "source": [
    "If you have downloaded LegPy in you PC and are runnig this notebook with Jupyter or a similar application, ignore this step. But if you are running this notebook from Google Colab, first you should clone the last version of the LegPy repository each time you inizialite Colab and add the directory where LegPy is located to the Python search path:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 1157,
     "status": "ok",
     "timestamp": 1714128306916,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "4VFq0S4y_q9w",
    "outputId": "f3177d8d-4291-4708-b3d1-7ff31ae7c57c"
   },
   "outputs": [],
   "source": [
    "#Only for Colab users!!\n",
    "!git clone https://github.com/JaimeRosado/LegPy.git\n",
    "import sys\n",
    "sys.path.insert(0,'/content/LegPy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TRQNk7zC_q9w"
   },
   "source": [
    "Then, import the necessary modules to run this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 1928,
     "status": "ok",
     "timestamp": 1714128315602,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "zm_QEL5L_q9w"
   },
   "outputs": [],
   "source": [
    "import LegPy as lpy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_G_eNS8bt1nK"
   },
   "source": [
    "The user has to \"construct\" four objects that are the main ingredients of the MC simulation:\n",
    "\n",
    "- The composition of the medium and their physical properties at microscopic level.\n",
    "- The geometry of the medium.\n",
    "- The energy spectrum of the beam crossing the medium.\n",
    "- The geometry of the beam and type of particles (photons or electrons)\n",
    "\n",
    "Let's start."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "W5mjzSpmhNkO"
   },
   "source": [
    "## Construction of the medium I: Composition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "piCc52MuALhl"
   },
   "source": [
    "LegPy provides the necessary data related to the interaction of photons and electrons with a large number of atoms, compounds and mixtures from the National Institute of Standards and Technology (NIST Database) https://www.nist.gov/.\n",
    "\n",
    "Select your media from this list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 302,
     "status": "ok",
     "timestamp": 1714128324436,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "AHhmw0eE_q9y",
    "outputId": "16c1770c-a0c9-48da-c9e2-bee203bd2e0a"
   },
   "outputs": [],
   "source": [
    "lpy.List_Media()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "k-tKx2TV_q9y"
   },
   "source": [
    "You also may search a medium name containing a string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 112
    },
    "executionInfo": {
     "elapsed": 338,
     "status": "ok",
     "timestamp": 1714128354331,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "hIDlQSz-_q9y",
    "outputId": "b1070eda-bb27-491a-ca06-54ed2661a523"
   },
   "outputs": [],
   "source": [
    "lpy.List_Media('water')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3RdovAuXnNWb"
   },
   "source": [
    "**A) Construct the medium providing the following data:**\n",
    "\n",
    "- `name` (as written in the above list): 'Al', 'Ge', 'Bone, Cortical (Icrp)', 'Sodium Iodide', 'C-552 Air-Equivalent Plastic', 'Water, Liquid'...\n",
    "- `density` (g/cm$^3$): Optional. Default value is stored in data files. Different values can be used for testing purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 312,
     "status": "ok",
     "timestamp": 1714128374811,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "M3OaMFC4A31p"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Aluminum Oxide')\n",
    "#medium = lpy.Medium(name='Al', density=2.7)\n",
    "#medium = lpy.Medium(name='Bone, Cortical (Icrp)', density=1.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "X47Mdq6M_q9z"
   },
   "source": [
    "In case of a simulation on a **two-media** geometry, both have to be constructed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 464,
     "status": "ok",
     "timestamp": 1714129735943,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "_cdDQfqT_q9z"
   },
   "outputs": [],
   "source": [
    "medium_1 = lpy.Medium(name='Aluminum Oxide')\n",
    "medium_2 = lpy.Medium(name='Al', density=2.7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xi3S1RIYbqz1"
   },
   "source": [
    "### Plot attenuation coefficients vs E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gHrluL3cb-o1"
   },
   "source": [
    "You can plot the components of the attenuation coeficients vs energy of a medium. In the next command you have to provide:\n",
    "\n",
    "- `energies` : array of energies in MeV, better in log scale.\n",
    "- `l_style` : line style for plotting ('', ':', etc.), defalut '' (solid line).\n",
    "- `ph` : True or False, default to True.\n",
    "- `inc` : True or False, default to True.\n",
    "- `coh` : True or False, default to True. Only for NIST option.\n",
    "- `pair` : True or False, default to True. Only for NIST option.\n",
    "- `tot` : True or False, default to True.\n",
    "\n",
    "Several media can be plotted in the same figure so you can define several media above (with different names) and compare attenuation coefficients.\n",
    "\n",
    "This step is optional so you can skip it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 476
    },
    "executionInfo": {
     "elapsed": 1201,
     "status": "ok",
     "timestamp": 1714128409901,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "jcJlsgcWpXkM",
    "outputId": "08ee78a7-5cb7-4c7a-9a7d-56bbe1b3ddf0",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "E1 = 0.01 # MeV\n",
    "E2 = 20. # MeV\n",
    "energy_range = np.logspace(np.log10(E1), np.log10(E2), num=150) # 150 points in a log-scale E(MeV) between E1 and E2\n",
    "medium.plot_mu(energies=energy_range)\n",
    "#medium2.plot_mu(energies=energy_range, l_style=':', tot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xnhjsMQD_q90"
   },
   "source": [
    "### Plot CSDA ranges vs E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "86H8C99s_q90"
   },
   "source": [
    "You can plot the continuous slowing down approximation (CSDA) range vs energy of the selected media. In the next command you have to provide:\n",
    "\n",
    "- `energies` : array of energies in MeV, better in log scale.\n",
    "- `l_style` : line style for plotting ('', ':', etc.), defalut '' (solid line).\n",
    "- `units` : 'cm' or 'gcm2'. Default is 'cm'.\n",
    "\n",
    "Several media can be plotted in the same figure so you can define several media above (with different names) and compare ranges (in the same units).\n",
    "\n",
    "This step is optional so you can skip it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 476
    },
    "executionInfo": {
     "elapsed": 1651,
     "status": "ok",
     "timestamp": 1714128421097,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "O8c5UsuU_q90",
    "outputId": "39675b92-48bd-42e4-be73-23cfe1f95245"
   },
   "outputs": [],
   "source": [
    "E1 = 0.01 # MeV\n",
    "E2 = 20. # MeV\n",
    "energy_range = np.logspace(np.log10(E1), np.log10(E2), num=150) # 150 points in a log-scale E(MeV) between E1 and E2\n",
    "medium.plot_R(energies=energy_range)\n",
    "#medium2.plot_R(energies=energy_range, l_style = ':')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4baNPkd-hNkR"
   },
   "source": [
    "## Construction of the medium II: Geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QWEfDqSAUmYq"
   },
   "source": [
    "Several geometries are available: cylinder, orthohedron and sphere.\n",
    "\n",
    "**Note**: *This geometry can be filled either with **one** single **medium** or contained **two media** splitted by an interface surface. The location of this surface has to be given for the construction of the geometry object*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SHSkFDEg_q90"
   },
   "source": [
    "**B) Select your geometry**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "z19PCvmQBt73"
   },
   "source": [
    "### Cylinder\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bBc-7XbDqUNS"
   },
   "source": [
    "Cylinder oriented with its axis along the z axis and its base at z=0. You have to provide:\n",
    "\n",
    "- `r` or `diam` (cm): radius or diameter.\n",
    "- `z` (cm): height.\n",
    "\n",
    "In the **two-media** case the interface can be a horizontal surface at a given depth:\n",
    "\n",
    "- `z_int` (cm): depth of the interface, that is, medium_1 in [0, z_int] and medium_2 in [z_int, z].\n",
    "\n",
    "or a cylinder with a given radius:\n",
    "\n",
    "- `r_int` (cm): medium_1 in the radial interval [0, r_int] and medium_2 in [r_int, r]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MD8zBO1Wywzg"
   },
   "source": [
    "For this geometry, you may choose either cylindrical or cartesian voxelization for the energy deposit matrix. Cylindrical voxelization is appropriate for parallel beams along the z axis and isotropic sources located at the z axis. In this case, you have to input the number of intervals along the coordinates r and z:\n",
    "\n",
    "- `n_r`.\n",
    "- `n_z`.\n",
    "\n",
    "Cartesian voxelization can also be applied in any situation and medium geometry. Here, you have to provide:\n",
    "\n",
    "- `n_x`.\n",
    "- `n_y`.\n",
    "- `n_z`.\n",
    "\n",
    "Choose your option and construct the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 390,
     "status": "ok",
     "timestamp": 1714129807931,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "xg_NkIl5SfM2"
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='cylinder', z=2., r=1., n_z=50, n_r=1) # Cylindrical voxelization\n",
    "#geometry = lpy.Geometry(name='cylinder', z=3., r=2., r_int = 0.5, n_x=10, n_y=10, n_z=10) # Cartesian voxelization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Hzj1Nd_XWNqB"
   },
   "source": [
    "### Orthohedron"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sj4xnK6uy6fy"
   },
   "source": [
    "Rectangular parallelepiped oriented with its longitudinal axes parallel to the x, y, z axes. The center of bottom side is assumed to be at the origin of coordinates (0,0,0). In this geometry, only the cartesian voxelization is supported. You have to provide the dimensions of the orthohedron and the number of intervals along each axis:\n",
    "\n",
    "- `x` (cm).\n",
    "- `y` (cm).\n",
    "- `z` (cm).\n",
    "\n",
    "In the **two-media** case the interface is a horizontal surface at a given depth:\n",
    "\n",
    "- `z_int` (cm): depth of the interface, that is, medium_1 in [0, z_int] and medium_2 in [z_int, z].\n",
    "\n",
    "Cartesian voxelization is used in this geometry. Here, you have to provide:\n",
    "- `n_x`.\n",
    "- `n_y`.\n",
    "- `n_z`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 538,
     "status": "ok",
     "timestamp": 1714129560425,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "6R4XvE8gDsgt"
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='orthohedron', x=10., y=10., z=10., z_int = 3., n_x=10, n_y=10, n_z=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ctlYVumbYCCo"
   },
   "source": [
    "### Sphere\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Kv3HZ51rYIwe"
   },
   "source": [
    "Sphere centered at the origin of coordinates (0,0,0). Both cartesian and spherical voxelization can be chosen. So you have to provide either:\n",
    "\n",
    "- `r` or `diam` (cm).\n",
    "- `n_r`.\n",
    "\n",
    "or:\n",
    "\n",
    "- `r` (cm).\n",
    "- `n_x`.\n",
    "- `n_y`.\n",
    "- `n_z`.\n",
    "\n",
    "In the **two-media** case the interface can be a flat surface at a given depth:\n",
    "\n",
    "- `z_int` (cm): z location of the interface, that is, medium_1 in the z interval [-r, z_int] and medium_2 in [z_int, r].\n",
    "\n",
    "or a spherical surface with a given radius:\n",
    "\n",
    "- `r_int` (cm): medium_1 in the sphere [0, r_int] and medium_2 in [r_int, r]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 2,
     "status": "ok",
     "timestamp": 1714128963933,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "fQSk-0dqYtKd"
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='sphere', r=15.72, n_r=15) # Spherical voxelization\n",
    "#geometry = lpy.Geometry(name='sphere', diam=10., r_int = 3., n_x=10, n_y=10, n_z=10) # Cartesian voxelization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ioqRp_pWZXEh"
   },
   "source": [
    "#### Plot the geometry\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "d0kuAxnRazD9"
   },
   "source": [
    "Plot the geometry in the reference coordinate system. This step is optional."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 418
    },
    "executionInfo": {
     "elapsed": 2401,
     "status": "ok",
     "timestamp": 1714128968178,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "oFLFtkX9hNkT",
    "outputId": "d963ae18-9afa-42a6-e95a-34e4ccb81910"
   },
   "outputs": [],
   "source": [
    "geometry.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RMgOPU02hNkU"
   },
   "source": [
    "## Construction of the beam I:  Energy spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Hp_iJsra_q92"
   },
   "source": [
    "**C) Select the energy spectrum of the particle beam**\n",
    "\n",
    "from one of the following options:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ElTUCsPicN9g"
   },
   "source": [
    "\n",
    "\n",
    "### Monoenergetic\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sfi5hy6UhYjZ"
   },
   "source": [
    "Input parameters:\n",
    "- `E` (MeV)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 227,
     "status": "ok",
     "timestamp": 1714128993117,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "uq1Ja2EzceUk"
   },
   "outputs": [],
   "source": [
    "spectrum = lpy.Spectrum(name = 'mono', E = 1.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TWfX_KuhcmDC"
   },
   "source": [
    "### Multi-monoenergetic\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hCdDmNovhuJj"
   },
   "source": [
    "Input parameters:\n",
    "- `E_w`: energies (MeV) and their corresponding weights in a numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 367,
     "status": "ok",
     "timestamp": 1714128998650,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "bAP7MlLqdKnN"
   },
   "outputs": [],
   "source": [
    "E_w = np.array([[0.511, .80], [1.25, 0.20]]) # [[E1, w1], [E2, w2],....]\n",
    "spectrum = lpy.Spectrum(name = 'multi_mono', E_w = E_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QFByI27SdiWU"
   },
   "source": [
    "\n",
    "\n",
    "### Flat\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Aa3ot_kwiITJ"
   },
   "source": [
    "Input parameters:\n",
    "- `E_min` (MeV).\n",
    "- `E_max` (MeV)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 257,
     "status": "ok",
     "timestamp": 1714129002947,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "2to4vC6YeN1q"
   },
   "outputs": [],
   "source": [
    "spectrum = lpy.Spectrum(name = 'flat', E_min = 0.1, E_max = 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SlBr3spZegV7"
   },
   "source": [
    "\n",
    "\n",
    "### Gaussian profile.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2NDOyGW5iRTx"
   },
   "source": [
    "Input parameters:\n",
    "- `E_mean` (MeV).\n",
    "- `E_sigma` (MeV).\n",
    "\n",
    "Internal cut: 2 x E_mean > E > 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 1,
     "status": "ok",
     "timestamp": 1714129005820,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "THzC4uRUe46j"
   },
   "outputs": [],
   "source": [
    "spectrum = lpy.Spectrum(name = 'gaussian', E_mean = 0.5, E_sigma = 0.03)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Nvy01AJYepFV"
   },
   "source": [
    "\n",
    "\n",
    "### Exponential\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wRtmtrwpifb6"
   },
   "source": [
    "\n",
    "$I(E) \\propto  e^{-E/E_{ch}}$, with E_min < E < E_max.\n",
    "\n",
    "Input parameters:\n",
    "\n",
    "- `E_min` (MeV).\n",
    "- `E_max` (MeV).\n",
    "- `E_ch` (MeV)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 2,
     "status": "ok",
     "timestamp": 1714129008895,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "nShl5yRDfi50"
   },
   "outputs": [],
   "source": [
    "spectrum = lpy.Spectrum(name = 'exponential', E_min = 0.1, E_max = 1.0, E_ch = 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oSNwzgeJelu4"
   },
   "source": [
    "\n",
    "\n",
    "### Reciprocal\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oyiaVZ9Xirbn"
   },
   "source": [
    "$ I(E) \\propto \\frac{1}{E} $, with E_min < E < E_max.\n",
    "\n",
    "Input parameters:\n",
    "- `E_min` (MeV).\n",
    "- `E_max` (MeV)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 394,
     "status": "ok",
     "timestamp": 1714129014418,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "-liB0ZZgf5Gf"
   },
   "outputs": [],
   "source": [
    "spectrum = lpy.Spectrum(name = 'reciprocal', E_min = 0.01, E_max = 15.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Pze164hHejOV"
   },
   "source": [
    "\n",
    "\n",
    "### From a file\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8jXnDDzRjEtV"
   },
   "source": [
    "A number of realistic spectra are provided in LegPy\n",
    "\n",
    "- X-Ray spectra generated with SpekPy Web 2.08: For example 80kV_3.0Al is for a tube potential of 80kV on Tungsten (default) with a filter of 3.0 mm Aluminum. In 60kV-Mo_0.04Mo both target and filter are Molybdenum.\n",
    "- Gamma-ray spectra of clinical accelerators: For example 10MV_VARIAN is a photon spectrum of a VARIAN accelerator of 10MV, 6MeV_VARIAN is an electron spectrum of 6MeV. More details of these spectra are given in the header of the corresponding .txt file inside the folder LegPy/beam_spectra."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8jXnDDzRjEtV"
   },
   "source": [
    "Check the available spectra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 175
    },
    "executionInfo": {
     "elapsed": 337,
     "status": "ok",
     "timestamp": 1714129019539,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "wE3YmiR8_q95",
    "outputId": "16bf4717-fabb-465b-f31d-d9b2c093a5bd"
   },
   "outputs": [],
   "source": [
    "lpy.List_Spectra()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QNRyE_5__q95"
   },
   "source": [
    "and select one of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 234,
     "status": "ok",
     "timestamp": 1714129629413,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "l6J2fe7-_q95"
   },
   "outputs": [],
   "source": [
    "spectrum = lpy.Spectrum(name = 'from_file', file = '6MV_Elekta')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gIr8_kUY24O5"
   },
   "source": [
    "If you want to use your own spectrum file, you should pass its path to the parameter `file`. The file must be a txt file with two columns: Energy (MeV)   ------   Relative Intensity (au)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "twadoIsqO8v4"
   },
   "outputs": [],
   "source": [
    "#This only works if there is a spectrum file called my_spectrum.txt in the working directory\n",
    "#spectrum = lpy.Spectrum(name = 'from_file', file = 'test')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xUhWwRffhCes"
   },
   "source": [
    "#### Plot the spectrum\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kBUnKUoqjjfF"
   },
   "source": [
    "You can plot the energy spectrum of incident beam. Again, just to check it is OK.\n",
    "\n",
    "A number of photons are generated randomly following the requested spectrum in logaritmic scale. You should input:\n",
    "- `n_part` : number of particles, default to 10^5.\n",
    "- `n_bin` : number of intervals, default to 50."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 476
    },
    "executionInfo": {
     "elapsed": 1043,
     "status": "ok",
     "timestamp": 1714129125539,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "nyFE-jCvhNkW",
    "outputId": "f7ce1ea4-d8d0-4e2d-8286-c75fede3b646"
   },
   "outputs": [],
   "source": [
    "spectrum.plot(n_part = 100000, n_bin = 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZoKl1zjohNkX"
   },
   "source": [
    "## Construction of the beam II:  Geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7kWebZ6CHH6X"
   },
   "source": [
    "**D) Select the particle (photon or electron) and the geometrical properties of the beam.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TcM0zGho_q96"
   },
   "source": [
    "Input parameter in all beam geometries:\n",
    "- `particle`= 'photon' or 'electron'. Default is `particle` = 'photon'."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mHeYsE2q_q96"
   },
   "source": [
    "NOTE: In order not to waste computing time the beam geometry has to be defined in such a way that all particles reach the medium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GUgNKbLZlRbp"
   },
   "source": [
    "### Parallel beam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KtTGumdNmGlu"
   },
   "source": [
    "Parallel beam with entrance plane perpendicular to z axis. In general not applicable for the sphere.\n",
    "\n",
    "Input parameters:\n",
    "\n",
    "- `theta` (degrees) : zenith angle from z axis, default to 0.\n",
    "- `phi` (degrees) : azimuth angle from x axis, default to 0.\n",
    "- `p_in` (cm) : numpy array with the coordinates of the center of the beam cross section at the entrance plan, default to (0,0,0).\n",
    "- `diam` (cm) : beam diameter, default to 0. (i.e., pencil beam)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 297,
     "status": "ok",
     "timestamp": 1714129830750,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "nZBzLDafmpPQ"
   },
   "outputs": [],
   "source": [
    "beam = lpy.Beam(name = 'parallel')\n",
    "#beam = lpy.Beam(particle='electron', name = 'parallel')\n",
    "#beam = lpy.Beam(name = 'parallel', theta = 15.0, phi = 30.0, p_in = np.array([0.1, -0.1, 0.0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TXmzG5rKlq0-"
   },
   "source": [
    "###  Divergent beam / Isotropic source"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "F3Lp7jfo_q96"
   },
   "source": [
    "Three options are available:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bUmNjg7e_q97"
   },
   "source": [
    "1) **Divergent beam** with the focus located on the z(<0) axis. Not applicable for the sphere."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iGgWOlOa_q97"
   },
   "source": [
    "Input parameters:\n",
    "\n",
    "- `p_in` (cm): numpy array with the coordinates of the focus location below the XY plane (z<0). Default to x=y=0 (source on the z axis). Small (x,y) values are also possible but might increase the computing time.      \n",
    "- size of the field located on the XY plane, centered at (0,0):\n",
    "    - if circular, provide `diam` (cm)\n",
    "    - if rectangular, provide `x_ap`, `y_ap` (cm) size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 240,
     "status": "ok",
     "timestamp": 1714129321504,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "Z_ZLMiq4_q97"
   },
   "outputs": [],
   "source": [
    "z = 25. # cm\n",
    "diam = 1. # cm\n",
    "beam = lpy.Beam(name = 'isotropic', diam = diam, p_in = np.array ([0., 0., -z]))\n",
    "#beam = lpy.Beam(name = 'isotropic',  particle='electron', x_ap = 1., y_ap = 0.5, p_in = np.array ([1., 0., -z]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sFuygSTK_q97"
   },
   "source": [
    "2) Isotropic **source** located **outside the medium**. Not applicable for the sphere. Currently only available for point sources. \n",
    "\n",
    "It is geometrically equals to the **divergent beam**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZskAcS1M_q97"
   },
   "source": [
    "Input parameters:\n",
    "- `p_in` (cm): numpy array with the coordinates of the source location below the XY plane (z<0). Default to x=y=0 (source on the z axis). Small (x,y) values are also possible but might increase the computing time.\n",
    "- size of the entrance aperture located on the XY plane, centered at (0,0):\n",
    "    - if circular, provide `diam` (cm)\n",
    "    - if rectangular, provide `x_ap`, `y_ap` (cm) size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 1,
     "status": "ok",
     "timestamp": 1714129680110,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "Jzz7Opbh_q97"
   },
   "outputs": [],
   "source": [
    "z = .5 # cm\n",
    "x, y = 0.05, -0.01 # cm\n",
    "beam = lpy.Beam(name = 'isotropic', diam = 1., p_in = np.array ([x, y, -z]))\n",
    "#beam = lpy.Beam(name = 'isotropic', particle='electron', x_ap = 1., y_ap = 0.5, p_in = np.array ([0., 0., -z]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6nZoVqnG_q97"
   },
   "source": [
    "3) Isotropic **source** located **inside the medium**.\n",
    "\n",
    "Both point sources and uniformly distributed in a sphere, cylinder or orthohedron are allowed. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iGgWOlOa_q97"
   },
   "source": [
    "Input parameters:\n",
    "\n",
    "- `p_in` (cm): numpy array with the coordinates of the location of a point source (or the center of a finite size  source). Default to p_in = [0., 0., 0.] (source on the coordinates origin).\n",
    "- size of the source. Default is point source\n",
    "    - if orthohedron, provide `x_s`, `y_s`, `z_s` (cm) size. Only one of these parameters are allowed to be 0.\n",
    "    - if cylinder, provide `r_s`, `z_s` (cm) size. Null radius `r_s` is allowed.  \n",
    "    - if sphere, provide `r_s`, (cm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 430,
     "status": "ok",
     "timestamp": 1714129640916,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "TCvMVHyX_q97"
   },
   "outputs": [],
   "source": [
    "x, y, z = 0.3, -0.3, 0.5 # cm\n",
    "beam = lpy.Beam(name = 'isotropic', particle='photon', p_in = np.array ([x, y, z]), r_s = 5.)\n",
    "#beam = lpy.Beam(name = 'isotropic')\n",
    "#beam = lpy.Beam(name = 'isotropic', particle='electron', p_in = np.array ([0., 0., 0.]), z_s = 1., r_s = 5.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pLaOZ3pWEKlX"
   },
   "source": [
    "#### Plot a few tracks\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5lnvTX95EaN_"
   },
   "source": [
    "Check a few (50) beam tracks into the medium with the geometry you have just constructed. Optional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 418
    },
    "executionInfo": {
     "elapsed": 1317,
     "status": "ok",
     "timestamp": 1714129682781,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "RW9f9Kw7_q98",
    "outputId": "4c1ebec3-3c65-4b40-9ea1-a35ac272248b"
   },
   "outputs": [],
   "source": [
    "lpy.Plot_beam(medium, geometry, spectrum, beam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jNfsfv73E3CV"
   },
   "source": [
    "## Monte Carlo Simulation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KIUJYEbs_q98"
   },
   "source": [
    "It transports the particle beam through the medium (or media) defined above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iCnVe-ia4Uy3"
   },
   "source": [
    "Parameters to be provided:\n",
    "\n",
    "For both particle beams, **photons and electrons**:  \n",
    "\n",
    "- `n_part`: number of particles (simulation cases). Default is `n_part` = 100.\n",
    "- `E_cut`: Energy cut in MeV. Default is `E_cut` = 0.01 MeV.\n",
    "- `tracks`: Plot tracks of the particle beam (not advised for > 100 photons). Default is `tracks` = False.\n",
    "\n",
    "In case of a **photon** beam:\n",
    "\n",
    "- `n_ang`: Number of intervals to construct angular histograms of photons leaving the medium. Default is 20.\n",
    "- `n_E`: Number of energy intervals to construct several histograms. It applies to histogram of photons leaving the medium, spectra of absorbed energy and spectral fluence. Default is 20.\n",
    "- `n_zloc` : Number of locations along the z axis to evaluate the photon fluence. Default is 20. Note that `n_zloc` is different from `n_z` defined in the `Geometry` object.\n",
    "- `points` : Plot interaction points (not advised for > 100 photons). Default is `points` = False.\n",
    "- `gamma_out` : If True, a plot angle vs energy for outgoing photons is made. Default is `gamma_out` = False. Do not use with > 10000 beam photons.\n",
    "- `x_ray`: If False, X-rays from atomic de-excitation after photoelectric absorption are neglected and the whole photon energy is absorbed in the interaction point. Default is `x_ray` = True. Only K$\\alpha$ and K$\\beta$ photons are generated.\n",
    "- `e_transport`: If True, all secondary electrons (photoelectric absorption and Compton scattering) are transported. By default (`e_transport` = False) electron energy is absorbed at the interaction point.\n",
    "\n",
    "In case of an **electron** beam:\n",
    "\n",
    "- `n_ang`: Number of intervals to construct histograms of angular distribution of emerging electrons. Default is 20.\n",
    "- `n_zloc`: Number of z positions/depths to construct transmission plots. Default is 20.\n",
    "- `e_length`, `e_K` : The electron transport is performed by steps of either same length (`e_length`) or same energy loss fraction (`e_K`). By default, e_length is obtained from the voxel size, but these parameters may also be specified, if desired. The parameter `e_length` is in $\\mu$m."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jG-yMed8_q98"
   },
   "source": [
    "Next some examples:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0V3fERwH_q98"
   },
   "source": [
    "Simulation on **one** single **medium**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 4801,
     "status": "ok",
     "timestamp": 1714129715670,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "_yAH61S4_q98",
    "outputId": "99e5d74a-a2cf-402a-e53f-08ccc60df506"
   },
   "outputs": [],
   "source": [
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part=1000, n_E=50)\n",
    "#result = lpy.MC(medium, geometry, spectrum, beam, n_part=10000, e_transport = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "caNeUF5w_q98"
   },
   "source": [
    "In case of a **two-media** geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 34571,
     "status": "ok",
     "timestamp": 1714129883926,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "aaLmDEa2_q99",
    "outputId": "87908a34-165c-457e-ce9e-d98135206d43"
   },
   "outputs": [],
   "source": [
    "result = lpy.MC([medium_1, medium_2], geometry, spectrum, beam, n_part=1000)\n",
    "#result = lpy.MC([medium_1, medium_2], geometry, spectrum, beam, n_part=int(1.e5), n_zloc = 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RXdxUCH-_q99"
   },
   "source": [
    "Return:\n",
    "\n",
    "An object (dubbed **result**) containing\n",
    "\n",
    "If **photon** beam:\n",
    "\n",
    "- the spatial distribution of deposited energy\n",
    "- the angular and energy spectra of outgoing photons\n",
    "- the histogram of absorbed energy\n",
    "- the spectral fluence along the z axis.\n",
    "\n",
    "If **electron** beam:\n",
    "\n",
    "- the spatial distribution of deposited energy\n",
    "- the histogram of electron ranges\n",
    "- the angular histogram of backscattered electrons.\n",
    "\n",
    "Once this object is generated, this information can be plotted or stored in files (see examples below)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9Yy9VvYm_q9-"
   },
   "source": [
    "# Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "luabfS64_q9-"
   },
   "source": [
    "Once the user gets familiar with the construction of the above objects and the MC parameters, running a simulation can be easily performed. Next we show a few **examples**:  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GSfMUSu7_q9-"
   },
   "source": [
    "- If not yet you should import several modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 400,
     "status": "ok",
     "timestamp": 1714130430395,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "O-1Ogy8ZhNkY",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import LegPy as lpy\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 8253,
     "status": "ok",
     "timestamp": 1714130440227,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "Kq3_CAL7crfb",
    "outputId": "704a4f9d-75ce-4e49-f8e6-4eab5ee5944e"
   },
   "outputs": [],
   "source": [
    "#Only for Colab users!!\n",
    "!pip install xlsxwriter\n",
    "import xlsxwriter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AYBZWEhF_q9-"
   },
   "source": [
    "## 1) Photon beam - one medium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lBtQT0U-_q9-"
   },
   "source": [
    "### 1.1) Display of photons tracks and/or interaction points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 412
    },
    "executionInfo": {
     "elapsed": 6936,
     "status": "ok",
     "timestamp": 1714130542602,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "LEe4YUXY_q9-",
    "outputId": "9934fc82-cafb-414e-9e26-19cabc17a10e",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Sodium Iodide')\n",
    "geometry = lpy.Geometry(name='cylinder', diam=50., z=20., n_x=30, n_y=30, n_z=100)\n",
    "beam = lpy.Beam(name='parallel', theta=0.0, phi=0.0, diam=4., p_in=(0.,0.,0.))\n",
    "spectrum = lpy.Spectrum(name='exponential', E_min=0.1, E_max=2., E_ch=0.3)\n",
    "lpy.Plot_beam(medium, geometry, spectrum, beam, n_part=100, points=True) #tracks=True by default\n",
    "#Alternative:\n",
    "#result = lpy.MC(medium, geometry, spectrum, beam, n_part=100, tracks=True, points=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IOJ-G_Su_q9-"
   },
   "source": [
    "### 1.2) Plot of angle vs energy of outgoing photons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 476
    },
    "executionInfo": {
     "elapsed": 773,
     "status": "ok",
     "timestamp": 1714130547065,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "1WFbsJsq_q9-",
    "outputId": "d07dca42-4f72-4f73-8cb7-a093100636ac",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Bismuth Germanium Oxide')\n",
    "geometry = lpy.Geometry(name='cylinder', diam=50., z=1., n_x=30, n_y=30, n_z=100)\n",
    "beam = lpy.Beam(name='parallel', theta=0.0, phi=0.0)\n",
    "spectrum = lpy.Spectrum(name='mono', E=1.)\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part = 1000, gamma_out=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SwHm3n9f_q9-"
   },
   "source": [
    "### 1.3) Histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 2475,
     "status": "ok",
     "timestamp": 1714130583500,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "os8EInBW_q9-",
    "outputId": "5c13837e-b69d-4fff-8d17-19385977cac0",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Water, Liquid')\n",
    "geometry = lpy.Geometry(name='cylinder', diam=50., z=10., n_x=30, n_y=30, n_z=100)\n",
    "beam = lpy.Beam(name='parallel', theta=0.0, phi=0.0)\n",
    "spectrum = lpy.Spectrum(name='mono', E=1.)\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part=10000, n_ang=50, n_E=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "npoAWJSV_q9-"
   },
   "source": [
    "**Histograms of angular and energy distributions of outgoing photons as well as the spectrum of energy absorption can be**\n",
    "\n",
    "- plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 316
    },
    "executionInfo": {
     "elapsed": 4085,
     "status": "ok",
     "timestamp": 1714130590227,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "IR6d8RON_q9-",
    "outputId": "05331364-f940-458a-ba2a-94287a354e80"
   },
   "outputs": [],
   "source": [
    "result.plot_hists()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6Ch7MTZe_q9-"
   },
   "source": [
    "- stored in dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 206
    },
    "executionInfo": {
     "elapsed": 337,
     "status": "ok",
     "timestamp": 1714130601553,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "n5v-pR9b_q9-",
    "outputId": "9f705934-1b15-4005-dfd4-dc2333363dac"
   },
   "outputs": [],
   "source": [
    "ang_out_df = result.ang_out()\n",
    "ang_out_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 206
    },
    "executionInfo": {
     "elapsed": 3,
     "status": "ok",
     "timestamp": 1714130604775,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "vfbITV5E_q9_",
    "outputId": "7733b90e-dc9a-4053-c4aa-38e5f7fe0024"
   },
   "outputs": [],
   "source": [
    "E_out_df = result.E_out()\n",
    "E_out_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 206
    },
    "executionInfo": {
     "elapsed": 413,
     "status": "ok",
     "timestamp": 1714130609395,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "l91B4qXi_q9_",
    "outputId": "5478420c-18b4-44a8-9d59-7a083d7a9a57"
   },
   "outputs": [],
   "source": [
    "E_ab_df = result.E_ab()\n",
    "E_ab_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8WoNp6CS_q9_"
   },
   "source": [
    "- written in excel files (a sheet per histogram)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 300
    },
    "executionInfo": {
     "elapsed": 327,
     "status": "error",
     "timestamp": 1714130615845,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "1f00DApA_q9_",
    "outputId": "855341ef-6b29-4744-f874-b09576a586e3"
   },
   "outputs": [],
   "source": [
    "result.hists_to_excel(\"my_excel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nPoLiYn3_q9_"
   },
   "source": [
    "### 1.4) Energy deposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-bmjgb8F_q9_"
   },
   "source": [
    "#### 1.4.1) Cylindral symmetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 11867,
     "status": "ok",
     "timestamp": 1714130671041,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "briIIxYS_q9_",
    "outputId": "8c146d1f-5912-412d-ee4c-7d1aebabfe17"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Tissue, Soft (Icru Four-Component)')\n",
    "geometry = lpy.Geometry(name='cylinder', diam=50., z=100., n_r=30, n_z=100)\n",
    "beam = lpy.Beam(name='parallel', theta=0.0, phi=0.0, diam=10.)\n",
    "spectrum = lpy.Spectrum(name='exponential', E_min=0.1, E_max=2., E_ch=0.3)\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part = 10000, e_transport=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_q-4W2RM_q9_"
   },
   "source": [
    "#### Spatial distribution can be:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "u31gwUx0_q9_"
   },
   "source": [
    "- plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 4915,
     "status": "ok",
     "timestamp": 1714130675952,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "51a-YQqq_q9_",
    "outputId": "2e77f15b-4a5f-461d-f8a6-670336bfac80"
   },
   "outputs": [],
   "source": [
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "YRzIr7HO_q9_"
   },
   "source": [
    "- stored in dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 287
    },
    "executionInfo": {
     "elapsed": 314,
     "status": "ok",
     "timestamp": 1714130679336,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "XhNfejIE_q9_",
    "outputId": "d2afe419-0e3c-41fb-81b7-90b8f1855e4d"
   },
   "outputs": [],
   "source": [
    "E_dep_df = result.Edep_to_df()\n",
    "E_dep_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EUYZmv0X_q9_"
   },
   "source": [
    "- written in excel file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 244,
     "status": "ok",
     "timestamp": 1714130682911,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "2a41Yq4K_q9_",
    "outputId": "d317db34-c1f6-40f4-fb16-f04082acfdb3"
   },
   "outputs": [],
   "source": [
    "result.Edep_to_excel('my_excel')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZQqOlmBa_q9_"
   },
   "source": [
    "Note that above results are averaged over axial angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Egh8Ls15_q9_"
   },
   "source": [
    "#### 1.4.2) Cartesian symmetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 10122,
     "status": "ok",
     "timestamp": 1714131121503,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "yff90YP7_q-E",
    "outputId": "fa830f8f-811c-4fd6-8daf-08791a2d2c1a"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Tissue, Soft (Icru Four-Component)')\n",
    "geometry = lpy.Geometry(name='cylinder', diam=50., z=100., n_x=30, n_y=30, n_z=100)\n",
    "beam = lpy.Beam(name='parallel', theta=0.0, phi=0.0, diam=10.)\n",
    "spectrum = lpy.Spectrum(name='exponential', E_min=0.1, E_max=2., E_ch=0.3)\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part = 10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zQAksqAP_q-E"
   },
   "source": [
    "Spatial distribution can be:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JJUM5D_O_q-E"
   },
   "source": [
    "- plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 286
    },
    "executionInfo": {
     "elapsed": 1343,
     "status": "ok",
     "timestamp": 1714131122842,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "yYjUvx3O_q-E",
    "outputId": "1f947b16-9243-43e3-dff5-7cc0adb105b6"
   },
   "outputs": [],
   "source": [
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UbqYmERg_q-E"
   },
   "source": [
    "- stored in python binary file with extension .npy. This 3D energy distribution cannot be stored in an excel file or a dataframe, but can be saved in this file for later use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 5,
     "status": "ok",
     "timestamp": 1714131122842,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "8kd8TkXM_q-E"
   },
   "outputs": [],
   "source": [
    "result.Edep_to_npy('my_3D_Edep')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "D6Uc72xxKYFk"
   },
   "source": [
    "Total energy deposited in the medium can be calculated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 5,
     "status": "ok",
     "timestamp": 1714131122842,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "7AvLM0RWKYFo",
    "outputId": "cce305e1-13cb-431d-ddf3-488e532a3f9f"
   },
   "outputs": [],
   "source": [
    "Vp = geometry.delta_v # pixel volume (cm^3)\n",
    "Ed = result.Edep.sum() * Vp\n",
    "print('total energy deposit =', round(Ed, 3), 'keV')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "trk9AEj_KYFq"
   },
   "source": [
    "Personalized plots of the spatial distribution of Edep can be done with this tool:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 237
    },
    "executionInfo": {
     "elapsed": 2216,
     "status": "ok",
     "timestamp": 1714131134673,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "ehzV500_KYFq",
    "outputId": "cf3feaf1-1be6-4f9d-aae6-8954664140ad"
   },
   "outputs": [],
   "source": [
    "z_ind = [0, 4, 9, 13, 16] # index of z layers.\n",
    "prof_lev = [3.5, 4., 5., 6., 7.] # dosis levels (adjust according the color bar)\n",
    "result.plot_Edep_layers(axis=\"z\", indexes=z_ind, c_profiles=True, lev=prof_lev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 216
    },
    "executionInfo": {
     "elapsed": 2295,
     "status": "ok",
     "timestamp": 1714131182707,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "ks9p9ITEKYFr",
    "outputId": "67dc4457-02b5-4353-9316-ac585c2c4e7d"
   },
   "outputs": [],
   "source": [
    "x_ind = [2, 9, 10, 15, 20] # index of z layers.\n",
    "prof_lev = [3.5, 4., 5., 6., 7.] # dosis levels (adjust according the color bar)\n",
    "result.plot_Edep_layers(axis=\"x\", indexes=x_ind, c_profiles=True, lev=prof_lev)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fAoR-x3o_q-E"
   },
   "source": [
    "### 1.5) Fluence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 9120,
     "status": "ok",
     "timestamp": 1714131204962,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "Q6u56BhR_q-E",
    "outputId": "1e136741-187c-4def-b79f-7fef1d5b92c3"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Water, Liquid')\n",
    "geometry = lpy.Geometry(name='cylinder', diam=50., z=50., n_r=30, n_z=50)\n",
    "beam = lpy.Beam(name='parallel', theta=0.0, phi=0.0)\n",
    "spectrum = lpy.Spectrum(name='exponential', E_min=0.1, E_max=2., E_ch=0.3)\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part = 10000, n_zloc=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4YWRVNgT_q-E"
   },
   "source": [
    "Total and spectral photon fluence is evaluated along the z axis in n_zloc+1 locations (from z = 0 to z = maximum depth) by counting photons traversing surfaces of the size of the cross section of the voxels, corrected for the direction of each photon. In this way it is calculated for all these z location the number of photons per unit area perpendicular to the surface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gf7i4Mz4_q-E"
   },
   "source": [
    "- Fluence plot\n",
    "\n",
    "    - First histogram: total fluence as function of depth\n",
    "    - Second histogram: spectral fluence at the entrance (z=0)\n",
    "    - Third histogram: spectral fluence at the exit (maximum z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 316
    },
    "executionInfo": {
     "elapsed": 2152,
     "status": "ok",
     "timestamp": 1714131207112,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "ep8W9KYi_q-F",
    "outputId": "83228536-6ae3-4193-8c51-028b144b5e86"
   },
   "outputs": [],
   "source": [
    "result.plot_fluence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Ugct1xHh_q-F"
   },
   "source": [
    "- Fluence data can be stored in a dataframe or in an excel file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 287
    },
    "executionInfo": {
     "elapsed": 11,
     "status": "ok",
     "timestamp": 1714131207112,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "J2QxEOQa_q-F",
    "outputId": "18cf47e1-0d1a-4782-d157-b700308188bc"
   },
   "outputs": [],
   "source": [
    "fluence_df = result.fluence_to_df()\n",
    "fluence_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 667,
     "status": "ok",
     "timestamp": 1714131812589,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "QGPJd3o5_q-F",
    "outputId": "b4d3553d-7592-45c4-d628-f040dca87935"
   },
   "outputs": [],
   "source": [
    "result.fluence_to_excel(\"my_excel\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XPk5MiuJ_q-F"
   },
   "source": [
    "## 2) Electron beam - one medium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JfRfA6xF_q-F"
   },
   "source": [
    "### 2.1) Display of electron tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 418
    },
    "executionInfo": {
     "elapsed": 16092,
     "status": "ok",
     "timestamp": 1714131831741,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "hp487Jpx_q-F",
    "outputId": "ceb252ce-0159-4c95-f7ef-dc5c7f159670"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Al')\n",
    "geometry = lpy.Geometry(name='orthohedron', x =.2, y = 0.2, z=.2, n_x=30, n_y=30, n_z=100)\n",
    "spectrum = lpy.Spectrum(name='mono', E = 1.)\n",
    "beam = lpy.Beam(particle='electron', name='parallel', theta=30.0, phi=20.0, diam=.02, p_in=(-0.03, 0., 0.))\n",
    "lpy.Plot_beam(medium, geometry, spectrum, beam, n_part=50)\n",
    "#Alternative:\n",
    "#result = lpy.MC(medium, geometry, spectrum, beam, n_part=50, tracks=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HBZ2rDAe_q-F"
   },
   "source": [
    "### 2.2) Electron range and backscattering studies ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "euRKQ5Is_q-F"
   },
   "source": [
    "First, check the CSDA range for an appropriate choice of the size of the medium and the `e_length` (`e_K`) parameters   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 331,
     "status": "ok",
     "timestamp": 1714132135394,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "9htBTbUQ_q-F",
    "outputId": "7fea671b-2e13-487b-c0ad-326dad9285bd"
   },
   "outputs": [],
   "source": [
    "E = 1.0 #MeV\n",
    "spectrum = lpy.Spectrum(name='mono', E=E)\n",
    "medium = lpy.Medium(name='Al') # NIST\n",
    "e_data = medium.e_data\n",
    "CSDA = np.interp(E, e_data.E_ref, e_data.R_ref) # cm\n",
    "CSDAum = CSDA * 1.e4 # um\n",
    "print('CSDA = ', round(CSDA, 3), 'cm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "inVt-UhG_q-F"
   },
   "source": [
    "A step length of CSDA/100 (or `e_K` = 0.95) should be enough to get accurate results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 342,
     "status": "ok",
     "timestamp": 1714132137811,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "X7qmq1Rw_q-F"
   },
   "outputs": [],
   "source": [
    "e_length = CSDAum * 0.01 # um (1/100 of CSDA)\n",
    "#e_K = 0.95"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 109595,
     "status": "ok",
     "timestamp": 1714132249969,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "3DgAO_Yj_q-F",
    "outputId": "793cd0f3-96a1-4b42-eb59-a9c35449cddc"
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='cylinder', diam=.25, z=.25, n_x=20, n_y=20, n_z=20)\n",
    "beam = lpy.Beam(particle = 'electron', name='parallel')\n",
    "n_part = 3000\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part=n_part, n_zloc=50, e_length=e_length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "H1umQW4l_q-F"
   },
   "source": [
    "- plot histograms of range, transmission and backscattered angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 353
    },
    "executionInfo": {
     "elapsed": 1519,
     "status": "ok",
     "timestamp": 1714132251484,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "-aM3DZik_q-F",
    "outputId": "29764ef9-6543-42de-d264-64f63c909157"
   },
   "outputs": [],
   "source": [
    "result.plot_hists()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qGY3oPO0_q-F"
   },
   "source": [
    "#### 2.2.1) Electron range ####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kYeQLtEa_q-F"
   },
   "source": [
    "Depending on the practical case, two different definitions of electron range R are used:\n",
    "\n",
    "- The depth z of the electron at the end of its path. This definition is used for the first histogram of the above figure.\n",
    "- The maximum depth reached by the electron. This definition is associated to the transmission curve obtained experimentally when the number of electrons traversing layers of several depths z are measured and is used for the second histogram of the above figure.\n",
    "\n",
    "Differences between both ranges might be non-negligible in cases with strong backscattering.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WTUNJRPt_q-G"
   },
   "source": [
    "Both definitions of range can be computed and stored in a dataframe with three components: x = R(cm), y = number of electrons, z = fraction of electrons.\n",
    "\n",
    "Check its shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 206
    },
    "executionInfo": {
     "elapsed": 327,
     "status": "ok",
     "timestamp": 1714132264486,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "5HUtRN_R_q-G",
    "outputId": "7cdd207c-6f5f-4f28-e68e-0c0fbe270ac8"
   },
   "outputs": [],
   "source": [
    "range_df = result.final_z()\n",
    "#range_df = result.max_z()\n",
    "range_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "96aV9DaI_q-G"
   },
   "source": [
    "Plot of R distribution and the corresponding integral function of the electron fraction vs R in either definition. Similar to to the plots obtained above but you can personalized them.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 881
    },
    "executionInfo": {
     "elapsed": 771,
     "status": "ok",
     "timestamp": 1714132269164,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "YB2gBRVK_q-G",
    "outputId": "c3c83da5-e128-455b-d191-e4c7990a85cb"
   },
   "outputs": [],
   "source": [
    "range_df.plot(kind='scatter', x=0, y=1);\n",
    "range_df.plot(kind='scatter', x=0, y=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6YeIkQxS_q-G"
   },
   "source": [
    "Calculation of the **extrapolated range** and other parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 505
    },
    "executionInfo": {
     "elapsed": 659,
     "status": "ok",
     "timestamp": 1714132275361,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "z8XEao2h_q-G",
    "outputId": "0f64da06-a96a-4d65-bfed-f9beb20eef8b"
   },
   "outputs": [],
   "source": [
    "ext_R, mode, av = lpy.ext_range(range_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mZJUkgfg_q-G"
   },
   "source": [
    "This funcion is also incorporated in the result object. The range `definition` should be set to \"final\" (default) or \"max\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 505
    },
    "executionInfo": {
     "elapsed": 350,
     "status": "ok",
     "timestamp": 1714132278820,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "NmbchCM5_q-G",
    "outputId": "7debafb9-71d6-4fad-a029-cc4b22159641"
   },
   "outputs": [],
   "source": [
    "ext_R, mode, av = result.ext_range(definition=\"max\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FtlhEBk8_q-H"
   },
   "source": [
    "#### 2.2.2) Electron Backscattering ####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UWqURoFY_q-H"
   },
   "source": [
    "The angle distribution of backscattered electrons can be computed and stored in a dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 206
    },
    "executionInfo": {
     "elapsed": 252,
     "status": "ok",
     "timestamp": 1714132310043,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "USxoIzW0_q-H",
    "outputId": "c853a1e2-3e06-4d66-8bfa-5d08482bbb13"
   },
   "outputs": [],
   "source": [
    "back = result.backscattering()\n",
    "back.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "L9cVfvEm_q-H"
   },
   "source": [
    "The backscattering coefficiente, b, is the fraction of backscattered electrons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 316,
     "status": "ok",
     "timestamp": 1714132313283,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "Wdkz61Yj_q-H",
    "outputId": "28d8cf7e-1c0e-4682-e60b-42cea6895688"
   },
   "outputs": [],
   "source": [
    "b = back.sum()[1]/n_part\n",
    "print('b = ', round(100.*b, 2), '%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5PCet4bn_q-H"
   },
   "source": [
    "Plot the angular distribution of backscattered electrons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 449
    },
    "executionInfo": {
     "elapsed": 638,
     "status": "ok",
     "timestamp": 1714132320178,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "h6D9ky5I_q-H",
    "outputId": "411f6e55-ca4c-461e-c968-b50d08a72345"
   },
   "outputs": [],
   "source": [
    "back.plot(kind='scatter', x=0, y=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DxciNpw3_q-H"
   },
   "source": [
    "### 2.3) Energy deposition ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aLSv_QE0_q-H"
   },
   "source": [
    "#### 2.3.1) Cylindrical symmetry   ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 24525,
     "status": "ok",
     "timestamp": 1714132348655,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "CIHGMtQK_q-H",
    "outputId": "62f14b1b-fb32-41bd-8e3d-96929edf88d9"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Al')\n",
    "geometry = lpy.Geometry(name='cylinder', r =.2, z=.2, n_r=20, n_z=20)\n",
    "spectrum = lpy.Spectrum(name='mono', E=1.)\n",
    "beam = lpy.Beam(particle='electron', name='parallel', diam=0.02)\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9B7xsEBg_q-H"
   },
   "source": [
    "Spatial distribution can be:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Fvxk48xP_q-H"
   },
   "source": [
    "- plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 4214,
     "status": "ok",
     "timestamp": 1714132352867,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "7TLQL8L6_q-H",
    "outputId": "87e8b5e8-3be6-4e2c-d377-a602fff5c4c3"
   },
   "outputs": [],
   "source": [
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3NukWFwl_q-H"
   },
   "source": [
    "- stored in dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 258
    },
    "executionInfo": {
     "elapsed": 1397,
     "status": "ok",
     "timestamp": 1714132396885,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "K5sTjF8M_q-H",
    "outputId": "461c4193-b849-4dfb-8bf3-57efc79a1f10"
   },
   "outputs": [],
   "source": [
    "E_dep_df = result.Edep_to_df()\n",
    "E_dep_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6OOsvk8l_q-H"
   },
   "source": [
    "- written in excel file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 315,
     "status": "ok",
     "timestamp": 1714132400673,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "nCxpD-cs_q-H",
    "outputId": "8953c514-a2d9-42ff-d5a1-36c38d4d776d"
   },
   "outputs": [],
   "source": [
    "result.Edep_to_excel('my_excel')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5_Rl4WGB_q-H"
   },
   "source": [
    "Note that above results are averaged over axial angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zGw0adUK_q-I"
   },
   "source": [
    "#### 2.3.2) Cartesian symmetry   ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 28152,
     "status": "ok",
     "timestamp": 1714132466715,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "dpTIvI7T_q-I",
    "outputId": "8f3ccbdd-d4b3-4f81-999e-bf9eddcc2e70"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Al')\n",
    "geometry = lpy.Geometry(name='cylinder', r=0.2, z=0.2, n_x=40, n_y=40, n_z=20)\n",
    "spectrum = lpy.Spectrum(name='mono', E=1.)\n",
    "beam = lpy.Beam(particle='electron', name='parallel', diam=0.02)\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_r7OGHRW_q-I"
   },
   "source": [
    "Spatial distribution can be:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "o6xspCoS_q-I"
   },
   "source": [
    "- plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 274
    },
    "executionInfo": {
     "elapsed": 1976,
     "status": "ok",
     "timestamp": 1714132468690,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "eITpGVNa_q-I",
    "outputId": "5814ec6a-5d4c-4b18-9147-3815d854f829"
   },
   "outputs": [],
   "source": [
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "w1yO3fQt_q-I"
   },
   "source": [
    "- stored in python binary file with extension .npy. Note that this 3D energy distribution cannot be stored in an excel file or a dataframe, but can be saved in this file for later use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 3,
     "status": "ok",
     "timestamp": 1714132468690,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "-ert928B_q-I"
   },
   "outputs": [],
   "source": [
    "result.Edep_to_npy('my_3D_Edep')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1CZVdlC5_q-I"
   },
   "source": [
    "Total energy deposited in the medium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 2,
     "status": "ok",
     "timestamp": 1714132468690,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "QqbDvmP4_q-I",
    "outputId": "d2bfd158-1aa0-40c4-f730-2819f94ccd95"
   },
   "outputs": [],
   "source": [
    "Vp = geometry.delta_v # pixel volume (cm^3)\n",
    "Ed = result.Edep.sum() * Vp\n",
    "print('total energy deposit =', round(Ed, 3), 'keV')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EBewd0sq_q-I"
   },
   "source": [
    "Personalized plots of the spatial distribution of Edep can be done with this tool:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 232
    },
    "executionInfo": {
     "elapsed": 2227,
     "status": "ok",
     "timestamp": 1714132470916,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "9uwfjmfL_q-I",
    "outputId": "037890b0-9714-463a-c7e4-628942f9639b"
   },
   "outputs": [],
   "source": [
    "z_ind = [0, 4, 9, 13, 16] # index of z layers.\n",
    "prof_lev = [3.5, 4., 5., 6., 7.] # dosis levels (adjust according the color bar)\n",
    "result.plot_Edep_layers(axis=\"z\", indexes=z_ind, c_profiles=True, lev=prof_lev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 316
    },
    "executionInfo": {
     "elapsed": 2164,
     "status": "ok",
     "timestamp": 1714132473078,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "7BbO92vi_q-I",
    "outputId": "a3eff8b1-0451-4bda-ee98-dffc7602bb99"
   },
   "outputs": [],
   "source": [
    "x_ind = [4, 14, 20, 30, 35] # index of z layers.\n",
    "prof_lev = [3.5, 4., 5., 6., 7.] # dosis levels (adjust according the color bar)\n",
    "result.plot_Edep_layers(axis=\"x\", indexes=x_ind, c_profiles=True, lev=prof_lev)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vfgUfxzK_q-I"
   },
   "source": [
    "### 2.4) Backscattered energy ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WMCwUQpS_q-I"
   },
   "source": [
    "Similar to the above case (both cylndrical or cartesian voxelization) but forcing the electron beam to start its path in some point inside the medium     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zx9M7EzC_q-I"
   },
   "source": [
    "#### 2.4.1) Cylindrical symmetry   ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 29058,
     "status": "ok",
     "timestamp": 1714132731238,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "BCcZ8kCy_q-I",
    "outputId": "6383be6e-7a54-4633-846f-da62d414c317"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Al')\n",
    "geometry = lpy.Geometry(name='cylinder', r=.2, z=.30, n_r=20, n_z=30)\n",
    "spectrum = lpy.Spectrum(name='mono', E=1.)\n",
    "z_in = 0.1\n",
    "beam = lpy.Beam(particle='electron', name='parallel', diam=0.02, p_in=np.array([0., 0., z_in]))\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 4794,
     "status": "ok",
     "timestamp": 1714132736027,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "tzZR2-5R_q-I",
    "outputId": "7088f474-a2c6-4261-b2b1-5b6ef4a3fa17"
   },
   "outputs": [],
   "source": [
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "agRfGkEw_q-J"
   },
   "source": [
    "#### 2.4.2) Cartesian symmetry   ####"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 26587,
     "status": "ok",
     "timestamp": 1714132763421,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "G7CU4R-F_q-J",
    "outputId": "ca420bab-4508-42eb-ba94-e496463d4417"
   },
   "outputs": [],
   "source": [
    "medium = lpy.Medium(name='Al')\n",
    "geometry = lpy.Geometry(name='cylinder', r=0.2, z=0.3, n_x=40, n_y=40, n_z=30)\n",
    "spectrum = lpy.Spectrum(name='mono', E=1.)\n",
    "z_in = 0.1\n",
    "beam = lpy.Beam(particle='electron', name='parallel', diam=0.02, p_in=np.array([0., 0., z_in]))\n",
    "result = lpy.MC(medium, geometry, spectrum, beam, n_part=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 274
    },
    "executionInfo": {
     "elapsed": 1841,
     "status": "ok",
     "timestamp": 1714132765256,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "9qXWWaf-_q-J",
    "outputId": "e80c9fde-4d19-4d3b-8b45-e48ae6c80abe"
   },
   "outputs": [],
   "source": [
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "OIqFcNWQ_q-J"
   },
   "source": [
    "Personalized plots of the spatial distribution of Edep can be done with this tool:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 232
    },
    "executionInfo": {
     "elapsed": 3151,
     "status": "ok",
     "timestamp": 1714132768405,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "jZT6V_cM_q-J",
    "outputId": "dad2f850-b3af-4e9d-83b7-c75faa632f78"
   },
   "outputs": [],
   "source": [
    "z_ind = [4, 9, 10, 13, 25] # index of z layers.\n",
    "prof_lev = [3.5, 4., 5., 6., 7.] # dosis levels (adjust according the color bar)\n",
    "result.plot_Edep_layers(axis=\"z\", indexes=z_ind, c_profiles=True, lev=prof_lev)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 267
    },
    "executionInfo": {
     "elapsed": 2389,
     "status": "ok",
     "timestamp": 1714132770791,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "51-bIudt_q-J",
    "outputId": "599f4145-5d34-4584-f0a6-4857122def72"
   },
   "outputs": [],
   "source": [
    "x_ind = [4, 14, 20, 30, 35] # index of z layers.\n",
    "prof_lev = [3.5, 4., 5., 6., 7.] # dosis levels (adjust according the color bar)\n",
    "result.plot_Edep_layers(axis=\"x\", indexes=x_ind, c_profiles=True, lev=prof_lev)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KNoccaWE_q-J"
   },
   "source": [
    "The energy deposited at z < z_in is the backscattered energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 5,
     "status": "ok",
     "timestamp": 1714132770791,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "cvISe4Vs_q-J",
    "outputId": "950ab0d6-c888-40fc-9961-032fef3ce717"
   },
   "outputs": [],
   "source": [
    "Vp = geometry.delta_v # pixel volume (cm^3)\n",
    "Edep = result.Edep\n",
    "Ed = Edep.sum() * Vp #total Edep\n",
    "Eb = Edep[:, :, 0:10].sum() * Vp #back Edep\n",
    "Ef = Edep[:, :, 10:30].sum() * Vp #for Edep\n",
    "\n",
    "print('Fraction of backsacttered energy =', round(100. * Eb / Ed, 2), '%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RhKryCoR_q-J"
   },
   "source": [
    "## 3) Simulations on two-media geometries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VLrRz85x_q-J"
   },
   "source": [
    "The simulation on geometries with two media can be performed in a similar way but you have to:\n",
    "\n",
    "- Define (construct) both media.\n",
    "- Give the position of the interface surface in your geometry.\n",
    "\n",
    "Note that a common voxelization in used in both media. The result can be analyze with the same tools described  above.\n",
    "\n",
    "Next you can find some examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WKSn0r6e_q-J"
   },
   "source": [
    "### 3.1) A divergent photon beam on a two-media cylinder ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 395,
     "status": "ok",
     "timestamp": 1714132808537,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "hxqLyU4C_q-J"
   },
   "outputs": [],
   "source": [
    "water = lpy.Medium(name='Water, Liquid') # NIST\n",
    "gold = lpy.Medium(name='Au') # NIST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 290,
     "status": "ok",
     "timestamp": 1714132810624,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "hAfuD37N_q-J"
   },
   "outputs": [],
   "source": [
    "E = 1. # MeV\n",
    "spectrum = lpy.Spectrum(name='mono', E = E)\n",
    "beam = lpy.Beam(name = 'isotropic', diam = 2.0, p_in = np.array ([0., 0., -5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IFq4yqgx_q-J"
   },
   "source": [
    "#### 3.1.1) z interface   ####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "s7M9_RG0_q-J"
   },
   "source": [
    "Interface: horizontal plane at z = z_int; water in [0, z_int] and gold in [z_int, z]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 6617,
     "status": "ok",
     "timestamp": 1714132821588,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "oQN_5m8r_q-J",
    "outputId": "f3f7a4d0-ba9c-419c-ada4-6fedc86a793a",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='cylinder', z=2., r=2., n_z=20, n_r=20, z_int=1.5)\n",
    "result = lpy.MC([water, gold], geometry, spectrum, beam, n_part=10000)\n",
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rnDbxAn4_q-J"
   },
   "source": [
    "#### 3.1.2) r interface   ####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zI5fD3_9_q-J"
   },
   "source": [
    "Interface: Inner cylinder of radius equals r_int; water in [0, r_int] and gold in [r_int, r]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 6084,
     "status": "ok",
     "timestamp": 1714132832394,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "OGnN5AG6_q-J",
    "outputId": "d7d4eea9-433c-46e1-a890-7dd163308f7b",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='cylinder', z=2., r=2., n_z=20, n_r=20, r_int=1.)\n",
    "result = lpy.MC([water, gold], geometry, spectrum, beam, n_part=10000)\n",
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9c-aOqE0_q-K"
   },
   "source": [
    "### 3.2) A paralell electron beam on a two-media orthohedron ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 242,
     "status": "ok",
     "timestamp": 1714132841811,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "nCPUxg6w_q-K"
   },
   "outputs": [],
   "source": [
    "water = lpy.Medium(name='Water, Liquid') # NIST\n",
    "bone = lpy.Medium(name='Bone, Compact (Icru)') # NIST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 300,
     "status": "ok",
     "timestamp": 1714132844891,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "0z0KTOvz_q-K"
   },
   "outputs": [],
   "source": [
    "E = .75 # MeV\n",
    "spectrum = lpy.Spectrum(name='mono', E=E)\n",
    "beam = lpy.Beam(name='parallel')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 379
    },
    "executionInfo": {
     "elapsed": 5886,
     "status": "ok",
     "timestamp": 1714132853595,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "ywLf4iGe_q-K",
    "outputId": "c8f895d2-897d-4692-ba98-b58acc5f365b",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='orthohedron', x=1., y=1., z=5., n_x=19, n_y=19, n_z=25, z_int=3.)\n",
    "result = lpy.MC([water, gold], geometry, spectrum, beam, n_part=10000)\n",
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "B7edvhie_q-K"
   },
   "source": [
    "### 3.3) An isotropic photon-source in the center of a two media sphere. ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 264,
     "status": "ok",
     "timestamp": 1714132857836,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "GpGtYA-n_q-K"
   },
   "outputs": [],
   "source": [
    "water = lpy.Medium(name='Water, Liquid') # NIST\n",
    "gold = lpy.Medium(name='Au') # NIST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "executionInfo": {
     "elapsed": 1,
     "status": "ok",
     "timestamp": 1714132858573,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "JiGAu_0o_q-K"
   },
   "outputs": [],
   "source": [
    "E = 1. # MeV\n",
    "beam = lpy.Beam(name='isotropic', p_in=np.array([0., 0., 0.]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "twXS-Kd9_q-K"
   },
   "source": [
    "#### 3.3.1) r interface   ####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qJ8YVwhf_q-K"
   },
   "source": [
    "Interface: inner sphere of radius r_int(cm); medium_1 in the sphere [0, r_int] and medium_2 in [r_int, r]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 600
    },
    "executionInfo": {
     "elapsed": 6094,
     "status": "ok",
     "timestamp": 1714132888645,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "gOM__HRg_q-K",
    "outputId": "8622316f-6b89-4bb1-ce1c-6f128be3cfa0",
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='sphere', r=5., n_r=50, r_int=3.)\n",
    "result = lpy.MC([water, gold], geometry, spectrum, beam, n_part=10000)\n",
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EfqtRLzG_q-K"
   },
   "source": [
    "#### 3.3.2) z interface   ####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "G-8HDUp2_q-K"
   },
   "source": [
    "Interface: horizontal plane at z = z_int; water in -r > z > z_int and gold in z_int > z > r."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 600
    },
    "executionInfo": {
     "elapsed": 4168,
     "status": "ok",
     "timestamp": 1714132897930,
     "user": {
      "displayName": "FERNANDO ARQUEROS MARTINEZ",
      "userId": "11315362602337364761"
     },
     "user_tz": -120
    },
    "id": "C9k8vvfu_q-K",
    "outputId": "92f81c52-bbf3-405f-f3a2-8a282f1c2eb6",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "geometry = lpy.Geometry(name='sphere', diam=10., n_r=50, n_z=10, z_int=2.5)\n",
    "result = lpy.MC([water, gold], geometry, spectrum, beam, n_part=10000)\n",
    "result.plot_Edep()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}