{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "dd469178-bab9-4515-b3e5-ee0fae371906",
   "metadata": {},
   "source": [
    "# Working with the `Reaction` helper\n",
    "\n",
    "This notebook demonstrates the convenience methods exposed by `Reaction` for building reaction systems and combining model ingredients. It is especially useful if you want a simpler interface than wiring every component together by hand.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "18b024d7-3a06-419b-93fa-9af470d4010d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jitr.reactions as rx\n",
    "import jitr.utils.constants as constants\n",
    "import jitr.utils.mass as mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4f31683c-a746-4c67-8169-780f143b24a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['BMA',\n",
       " 'UNEDF2',\n",
       " 'BCPM',\n",
       " 'dz31',\n",
       " 'hfb31',\n",
       " 'ws4rbf',\n",
       " 'FRDM2012',\n",
       " 'SLY4',\n",
       " 'SKP',\n",
       " 'ame2020',\n",
       " 'SKM',\n",
       " 'UNEDF0',\n",
       " 'D1M',\n",
       " 'SVMIN',\n",
       " 'UNEDF1',\n",
       " 'HFB24']"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mass.mass_models()  # default is ame2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "74f28c46-a21f-4ca3-b195-629ec2edcf54",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "48-Ca(p,n)48-Sc\n"
     ]
    }
   ],
   "source": [
    "mass_model = \"ame2024\"\n",
    "rxn = rx.Reaction(\n",
    "    target=(48, 20),\n",
    "    projectile=(1, 1),\n",
    "    product=(1, 0),\n",
    "    residual=(48, 21),\n",
    "    mass_kwargs={\"model\": \"ame2020\"},\n",
    ")\n",
    "print(rxn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cab65e63-db45-43c6-b8d8-996177e47444",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle ^{48} \\rm{Ca}(p,n)^{48} \\rm{Sc}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import Math, display\n",
    "\n",
    "display(Math(rxn.reaction_latex))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c4897eee-e550-48a7-8b6b-5a4bbe7add62",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q-value: -0.5036 MeV\n"
     ]
    }
   ],
   "source": [
    "print(f\"Q-value: {rxn.Q:1.4f} MeV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "732e9e4c-54bc-493d-96f3-f437e7197d5f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ChannelKinematics(Elab=35, Ecm=34.27976496842483, mu=np.float64(951.4664814091528), k=np.float64(1.28287146975579), eta=np.float64(0.5485537668413036))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 35 MeV proton beam incident on 48Ca\n",
    "Elab = 35\n",
    "entrance_kinematics = rxn.kinematics(Elab)\n",
    "entrance_kinematics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4b39f4ec-6757-4607-bdd0-bff1375e88ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ChannelKinematics(Elab=27.671337282951978, Ecm=27.10121706917766, mu=np.float64(945.9090824120628), k=np.float64(1.1394155566936675), eta=np.float64(0.0))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the excitation energy of the state in 48Sc that is the isobaric analog to the 48Ca G.S.s\n",
    "# (above the 48Sc G.S.)\n",
    "# from https://www.sciencedirect.com/science/article/pii/S0092640X97907403\n",
    "Ex_IAS = 6.675\n",
    "exit_kinematics = rxn.kinematics_exit(\n",
    "    entrance_kinematics, residual_excitation_energy=Ex_IAS\n",
    ")\n",
    "exit_kinematics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bcfb9701-e656-4403-8ffd-01c42dc31911",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "assert np.isclose(exit_kinematics.Ecm, 27.101217)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4ec0471b-bb2c-4ace-b0c7-bb36abc44249",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proton mass: 938.2717  MeV/c^2 ~ 938.2717  MeV/c^2\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"proton mass: {rxn.projectile.m0:1.4f}  MeV/c^2 ~ {constants.MASS_P:1.4f}  MeV/c^2\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "532dbfbd-5522-4b3a-87a5-af3a5c15cba5",
   "metadata": {},
   "source": [
    "This small difference will be the result of some of the differences below when manually calculating vs using the mass model passed into `Reaction`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bb4d702c-5798-4331-a51a-9ac86d836307",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "target mass: 44657.2720  MeV/c^2 = 44657.2720  MeV/c^2\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"target mass: {rxn.target.m0:1.4f}  MeV/c^2 = {mass.mass(*rxn.target)[0]:1.4f}  MeV/c^2\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "823dff13-0fee-4977-ab1f-cda62759eb5b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reaction threshold: 15.8014 MeV = proton separation energy in 48-Ca: 15.8014 MeV\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"reaction threshold: {rxn.threshold:1.4f} MeV = proton separation energy \"\n",
    "    f\"in {rxn.target}: {mass.proton_separation_energy(*rxn.target)[0]:1.4f} MeV\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c4e9eeca-e6f7-44a9-845d-a5def13ac805",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "49-Sc"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rxn.compound_system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "63cf15ff-81f8-4c5a-8204-38d6e42c9967",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "threshold for p-removal from 49-Sc: 9.6261 MeV = proton separation energy in 49-Sc: 9.6261 MeV\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"threshold for {rxn.projectile}-removal from {rxn.compound_system}: {rxn.compound_system_threshold:1.4f} MeV = \"\n",
    "    f\"proton separation energy in {rxn.compound_system}:\"\n",
    "    f\" {mass.proton_separation_energy(*rxn.compound_system)[0]:1.4f} MeV\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "eb948788-5447-4312-aea0-63075df810fb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "48-Ca proton Fermi energy: -12.7138 MeV = -12.7138 MeV = -12.7138 MeV\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"{rxn.target} proton Fermi energy: {rxn.Ef:1.4f} MeV = \"\n",
    "    f\"{rxn.target.Efp:1.4f} MeV = \"\n",
    "    f\"{mass.proton_fermi_energy(*rxn.target)[0]:1.4f} MeV\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d693290d-5364-4a65-ba13-ec0d6a69e136",
   "metadata": {},
   "outputs": [],
   "source": [
    "rxn = rx.Reaction(\n",
    "    target=(50, 20),\n",
    "    projectile=(4, 2),\n",
    "    product=(2, 1),\n",
    "    mass_kwargs={\"model\": \"ame2020\"},\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a53be2bb-a752-48ee-9864-ba2a876411b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50-Ca(alpha,d)52-Sc"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rxn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0fd5e06e-08f6-4937-a8e8-341a102b861a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle ^{50} \\rm{Ca}(\\alpha,d)^{52} \\rm{Sc}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Math(rxn.reaction_latex))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "46c148fc-34a1-462c-beb4-599f6cb23209",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q value: -9.7765 MeV\n"
     ]
    }
   ],
   "source": [
    "print(f\"Q value: {rxn.Q:1.4f} MeV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c8764195-9ce8-4a94-94f7-32777e302ca7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "54-Ti"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rxn.compound_system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "4320ac52-7405-4fbf-9406-896461a3655e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "threshold for alpha from compound system: 8.5795 MeV = alpha separation energy in 54-Ti: 8.5795 MeV\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    f\"threshold for {rxn.projectile} from compound system: {rxn.compound_system_threshold:1.4f} MeV = \"\n",
    "    f\"{rxn.projectile} separation energy in {rxn.compound_system}:\"\n",
    "    f\" {rx.cluster_separation_energy(rxn.compound_system, rxn.projectile):1.4f} MeV\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1c678e3-4527-480c-92d2-7981344e6694",
   "metadata": {},
   "source": [
    "## Mix mass models when needed\n",
    "For example some of the models don't have some light particle masses, but `ame2020` doesn't have a `60-Ca` mass. We can use the `Nucleus` class to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "cef38ecb-37cd-47b5-9f20-55ba5f65ce09",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle ^3_1t$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================\n",
      "BMA     : nan\n",
      "UNEDF2  : nan\n",
      "BCPM    : nan\n",
      "dz31    : nan\n",
      "hfb31   : nan\n",
      "ws4rbf  : nan\n",
      "FRDM2012: nan\n",
      "SLY4    : nan\n",
      "SKP     : nan\n",
      "ame2020 : 2808.921118\n",
      "SKM     : nan\n",
      "UNEDF0  : nan\n",
      "D1M     : nan\n",
      "SVMIN   : nan\n",
      "UNEDF1  : nan\n",
      "HFB24   : nan\n"
     ]
    }
   ],
   "source": [
    "nuc = rx.Nucleus(3, 1)\n",
    "display(Math(f\"^{nuc.A}_{nuc.Z}\" + nuc.latex()))\n",
    "print(\"=====================\")\n",
    "\n",
    "for m in mass.__MASS_MODELS__:\n",
    "    print(f\"{m:8}: {mass.mass(*nuc, model=m)[0]:1.6f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "f2faca5e-5fad-4f70-84f0-0f3b2284fd6a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle ^{60} \\rm{Ca}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================\n",
      "BMA     : 55887.086308\n",
      "UNEDF2  : 55891.426141\n",
      "BCPM    : 55883.586141\n",
      "dz31    : 55891.798141\n",
      "hfb31   : 55889.466141\n",
      "ws4rbf  : 55891.702141\n",
      "FRDM2012: 55890.176141\n",
      "SLY4    : 55885.506141\n",
      "SKP     : 55882.246141\n",
      "ame2020 : nan\n",
      "SKM     : 55878.546141\n",
      "UNEDF0  : 55885.826141\n",
      "D1M     : 55891.736141\n",
      "SVMIN   : 55885.796141\n",
      "UNEDF1  : 55891.066141\n",
      "HFB24   : 55890.236141\n"
     ]
    }
   ],
   "source": [
    "nuc = rx.Nucleus(60, 20)\n",
    "display(Math(nuc.latex()))\n",
    "print(\"=====================\")\n",
    "for m in mass.__MASS_MODELS__:\n",
    "    print(f\"{m:8}: {mass.mass(*nuc, model=m)[0]:1.6f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6ff413a4-9bfd-45c6-803e-8e350e46cd2f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60-Ca(d,t)59-Ca"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rxn = rx.Reaction(\n",
    "    target=rx.Nucleus(60, 20, mass_kwargs={\"model\": \"UNEDF2\"}),\n",
    "    projectile=(2, 1),\n",
    "    product=(3, 1),\n",
    "    mass_kwargs={\"model\": \"ame2020\"},  # AME2020 for everything else\n",
    ")\n",
    "rxn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "3225e9f5-b7cc-4b79-b974-bceda890395b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.2259119999994255"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rxn.Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c3a36fe0-fc52-40ca-90a6-26ad543536bd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ChannelKinematics(Elab=240, Ecm=232.2075441296029, mu=np.float64(2027.0754691804402), k=np.float64(4.779677811174583), eta=np.float64(0.31367517605122724))"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "entrance_kinematics = rxn.kinematics(240)\n",
    "entrance_kinematics"
   ]
  }
 ],
 "metadata": {
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}