{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1b4c5f7f",
   "metadata": {},
   "source": [
    "# Synthetic Charged Particle Radiographs with Multiple Detector Layers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d77e602c",
   "metadata": {},
   "source": [
    "[charged_particle_radiography]: ../../ad/diagnostics/charged_particle_radiography/index.rst\n",
    "\n",
    "[Bragg peak]:\n",
    "https://en.wikipedia.org/wiki/Bragg_peak\n",
    "    \n",
    "Charged particle radiographs are often recorded on detector stacks with multiple layers (often either radiochromic film or CR39 plastic). Since charged particles deposit most of their energy at a depth related to their energy (the [Bragg peak]), each detector layer records a different energy band of particles. The energy bands recorded on each layer are further discriminated by including layers of filter material (eg. thin strips of aluminum) between the active layers. In order to analyze this data, it is necessary to calculate what range of particle energies are stopped in each detector layer.\n",
    "\n",
    "The [charged_particle_radiography] module includes a tool for calculating the energy bands for a given stack of detectors. This calculation is a simple approximation (more accurate calculations can be made using Monte Carlo codes), but is sufficient in many cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d24d2476",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tempfile import gettempdir\n",
    "\n",
    "import astropy.units as u\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from plasmapy.diagnostics.charged_particle_radiography.detector_stacks import (\n",
    "    Layer,\n",
    "    Stack,\n",
    ")\n",
    "from plasmapy.utils.data.downloader import Downloader"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91186291",
   "metadata": {},
   "source": [
    "[charged_particle_radiography]: ../../ad/diagnostics/charged_particle_radiography/index.rst\n",
    "\n",
    "[Layer]:\n",
    "../../api/plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Layer.rst#plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Layer\n",
    "\n",
    "[Stack]:\n",
    "../../api/plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Stack.rst#plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Stack\n",
    "\n",
    "[astropy Quantity]:\n",
    "https://docs.astropy.org/en/stable/api/astropy.units.quantity.Quantity.html\n",
    "\n",
    "[PSTAR]:\n",
    "https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html\n",
    "\n",
    "[ESTAR]:\n",
    "https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html\n",
    "\n",
    "The [charged_particle_radiography] module represents a detector stack with a [Stack] object, which contains an ordered list of [Layer] objects. In this example, we will create a list of [Layer] objects, use them to initialize a [Stack], then compute the energy bands.\n",
    "\n",
    "Each [Layer] is defined by the following required properties:\n",
    "- `thickness`: an [astropy Quantity] defining the thickness of the layer.\n",
    "- `active`: a boolean flag which indicates whether the layer is active (detector) or inactive (eg. a filter).\n",
    "- `energy_axis`: An [astropy Quantity] array of energies.\n",
    "- `stopping_power`: An [astropy Quantity] array containing the product of stopping power and material density at each energy in `energy_axis`. \n",
    "\n",
    "The stopping power in various materials for protons and electrons is tabulated by NIST in the [PSTAR] and [ESTAR] databases. The stopping powers are tabulated in units of MeV cm$^2$ / g, so the product of stopping power and density has units of MeV/cm. \n",
    "\n",
    "For this demonstration, we will load two stopping power tables, downloaded from [PSTAR], for aluminum and a human tissue equivalent that is a good approximation for many radiochromic films."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "073dad43",
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_dir = gettempdir()  # Get a temporary directory to save the files\n",
    "dl = Downloader(directory=temp_dir, validate=False)\n",
    "tissue_path = dl.get_file(\"NIST_PSTAR_tissue_equivalent.txt\")\n",
    "aluminum_path = dl.get_file(\"NIST_PSTAR_aluminum.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a4de6ef",
   "metadata": {},
   "source": [
    "[PSTAR]:\n",
    "https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html\n",
    "\n",
    "These files can be directly downloaded from [PSTAR], and the contents look like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "748834f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(tissue_path) as f:\n",
    "    for i in range(15):\n",
    "        print(f.readline(), end=\"\")\n",
    "print(\"...\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08fd3990",
   "metadata": {},
   "source": [
    "Now we will load the energy and stopping power arrays from the files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffc2bb16",
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = np.loadtxt(tissue_path, skiprows=8)\n",
    "energy_axis = arr[:, 0] * u.MeV\n",
    "tissue_density = 1.04 * u.g / u.cm**3\n",
    "tissue_stopping_power = arr[:, 1] * u.MeV * u.cm**2 / u.g * tissue_density\n",
    "\n",
    "arr = np.loadtxt(aluminum_path, skiprows=8)\n",
    "aluminum_density = 2.7 * u.g / u.cm**3\n",
    "aluminum_stopping_power = arr[:, 1] * u.MeV * u.cm**2 / u.g * aluminum_density"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59d7daf8",
   "metadata": {},
   "source": [
    "If we plot the stopping powers as a function of energy on a log-log scale, we see that they are non-linear, with much higher values at lower energies. This is why particles deposit most of their energy at a characteristic depth. As particles propagate into material, they lose energy slowly at first. But, as they loose energy, the effective stopping power increases expontentially. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6f59d48",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(energy_axis.value, tissue_stopping_power, label=\"Tissue equivalent\")\n",
    "ax.plot(energy_axis.value, aluminum_stopping_power, label=\"Aluminum\")\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xlabel(\"Energy (MeV)\")\n",
    "ax.set_ylabel(\"Stopping power (MeV/cm)\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fbc08a9",
   "metadata": {},
   "source": [
    "[HD-V2]:\n",
    "http://www.gafchromic.com/documents/gafchromic-hdv2.pdf\n",
    "Before creating a Stack, we will start by creating a list that represents a common type of radiochromic film, [HD-V2]. HD-V2 consists of two layers: a 12 $\\mu$m active layer deposited on a 97 $\\mu$m substrate. The stopping powers of both layers are similar to human tissue (by design, HD-V2 is commonly used for medical research)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e991901",
   "metadata": {},
   "outputs": [],
   "source": [
    "HDV2 = [\n",
    "    Layer(12 * u.um, energy_axis, tissue_stopping_power, active=True),\n",
    "    Layer(97 * u.um, energy_axis, tissue_stopping_power, active=False),\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a460c1ba",
   "metadata": {},
   "source": [
    "[Stack]:\n",
    "../../api/plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Stack.rst#plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Stack\n",
    "\n",
    "We will now create a list of layers of HDV2 separated by aluminum filter layers of various thickness, then use this list to instantiate a [Stack] object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc58debf",
   "metadata": {},
   "outputs": [],
   "source": [
    "layers = [\n",
    "    Layer(100 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(100 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(100 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(100 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(500 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(500 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(1000 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(1000 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(2000 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(2000 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "    *HDV2,\n",
    "    Layer(2000 * u.um, energy_axis, aluminum_stopping_power, active=False),\n",
    "]\n",
    "\n",
    "stack = Stack(layers)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10726a37",
   "metadata": {},
   "source": [
    "The number of layers, active layers, and the total thickness of the film stack are available as properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "559a1ff3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Number of layers: {stack.num_layers}\")\n",
    "print(f\"Number of active layers: {stack.num_active}\")\n",
    "print(f\"Total stack thickness: {stack.thickness:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a94123a",
   "metadata": {},
   "source": [
    "[deposition_curves()]: ../../api/plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Stack.rst#plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Stack.deposition_curves\n",
    "\n",
    "The curves of deposited energy per layer for a given array of energies can then be calculated using the [deposition_curves()] method. The stopping power for a particle can change significantly within a layer, so the stopping power is numerically integrated within each layer. The spatial resolution of this integration can be set using the `dx` keyword. Setting the `return_only_active` keyword means that deposition curves will only be returned for the active layers (inactive layers are still included in the calculation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e984e152",
   "metadata": {},
   "outputs": [],
   "source": [
    "energies = np.arange(1, 60, 0.1) * u.MeV\n",
    "deposition_curves = stack.deposition_curves(\n",
    "    energies, dx=1 * u.um, return_only_active=True\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b62ada4",
   "metadata": {
    "nbsphinx-thumbnail": {
     "output-index": 0,
     "tooltip": "Charged particle radiography film stacks"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "ax.set_title(\"Energy deposition curves\")\n",
    "ax.set_xlabel(\"Energy (MeV)\")\n",
    "ax.set_ylabel(\"Normalized energy deposition curve\")\n",
    "for layer in range(stack.num_active):\n",
    "    ax.plot(energies, deposition_curves[layer, :], label=f\"Layer {layer + 1}\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d3e83b5",
   "metadata": {},
   "source": [
    "We can look at the deposition curves including all layers (active and inactive) to see where particles of a given energy are deposited."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96609ac5",
   "metadata": {},
   "outputs": [],
   "source": [
    "deposition_curves_all = stack.deposition_curves(energies, return_only_active=False)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "ax.set_title(\"10 MeV Particles per layer\")\n",
    "ax.set_xlabel(\"Layer #\")\n",
    "ax.set_ylabel(\"Fraction of particles\")\n",
    "E0 = np.argmin(np.abs(10 * u.MeV - energies))\n",
    "ax.plot(np.arange(stack.num_layers) + 1, deposition_curves_all[:, E0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3ed8684",
   "metadata": {},
   "source": [
    "The sharp peaks correspond to the aluminum filter layers, which stop more particles than the film layers.\n",
    "\n",
    "The deposition curve is normalized such that each value represents the fraction of particles with a given energy that will be stopped in that layer. Since all particles are stopped before the last layer, the sum over all layers (including the inactive layers) for a given energy is always unity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b231fee1",
   "metadata": {},
   "outputs": [],
   "source": [
    "integral = np.sum(deposition_curves_all, axis=0)\n",
    "print(np.allclose(integral, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d96105ad",
   "metadata": {},
   "source": [
    "[FWHM]:\n",
    "https://en.wikipedia.org/wiki/Full_width_at_half_maximum\n",
    "\n",
    "[energy_bands()]: ../../api/plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Stack.rst#plasmapy.diagnostics.charged_particle_radiography.detector_stacks.Stack.energy_bands\n",
    "\n",
    "We can quantify the range of initial particle energies primarily recorded on each film layer by calculating the full with at half maximum ([FWHM]) of the deposition curve for that layer. This calculation is done by the [energy_bands()] method, which returns an array of the $\\pm$ FWHM values for each film layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db3e8a81",
   "metadata": {},
   "outputs": [],
   "source": [
    "ebands = stack.energy_bands([0.1, 60] * u.MeV, 0.1 * u.MeV, return_only_active=True)\n",
    "for layer in range(stack.num_active):\n",
    "    print(\n",
    "        f\"Layer {layer + 1}: {ebands[layer, 0].value:.1f}-{ebands[layer, 1].value:.1f} MeV\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b27de99",
   "metadata": {},
   "source": [
    "From this information, we can conclude see that the radiograph recorded on Layer 5 is primarily created by protons bwetween 12.5 and 12.8 MeV. This energy information can then be used to inform analysis of the images."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dev312",
   "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.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}