{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Models of Inverse Beta Decay\n",
    "\n",
    "The interactions module in ASTERIA contains two implementations of inverse beta decay:\n",
    "\n",
    "1. A simple parametric model.\n",
    "2. A tabular model based on a full calculation of the interaction.\n",
    "\n",
    "Both calculations can be found in [Strumia and Vissani, Phys. Lett. B 564:42, 2003](https://www.sciencedirect.com/science/article/pii/S0370269303006166)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from asteria import set_rcparams\n",
    "from asteria.interactions import InvBetaPar, InvBetaTab\n",
    "from snewpy.neutrino import Flavor\n",
    "\n",
    "import astropy.units as u\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup styles for Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_rcparams(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare Cross Sections and Mean Energies\n",
    "\n",
    "For the energies of interest to IceCube, the parametric and tabular models differ by a few percent in cross section and mean energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(12, 6), sharex=True,\n",
    "                         gridspec_kw={'height_ratios':[3,1], 'hspace':0, 'wspace': 0.25})\n",
    "\n",
    "xs_old = None\n",
    "lep_old = None\n",
    "enu = np.linspace(0., 200., 101) * u.MeV\n",
    "\n",
    "for ibd, style, lab in zip([InvBetaPar(), InvBetaTab()],\n",
    "                           ['-', '.'],\n",
    "                           ['Parametric (Eq. 25)', 'Table 1']):\n",
    "    xs = ibd.cross_section(Flavor.NU_E_BAR, enu)\n",
    "    lep = ibd.mean_lepton_energy(Flavor.NU_E_BAR, enu)\n",
    "    axes[0][0].plot(enu, xs, style, label=lab)\n",
    "    axes[0][1].plot(enu, lep, style, label=lab)\n",
    "\n",
    "    # Plot residuals\n",
    "    if xs_old is not None and lep_old is not None:\n",
    "        cut = xs_old != 0.\n",
    "        res = (xs[cut] - xs_old[cut])/xs_old[cut]\n",
    "        axes[1][0].plot(enu[cut], 1e2*res, style)\n",
    "        cut = lep_old != 0.\n",
    "        res = (lep[cut] - lep_old[cut])/lep_old[cut]\n",
    "        axes[1][1].plot(enu[cut], 1e2*res, style)\n",
    "    xs_old = xs\n",
    "    lep_old = lep\n",
    "\n",
    "axes[0,0].set_ylabel( r'$\\sigma(\\bar{\\nu}_e p\\rightarrow e^{+}n)$ [cm$^2$]', horizontalalignment='right', y=1)\n",
    "axes[0,0].set(yscale='log')\n",
    "axes[1,0].set_xlabel( r'$E_\\nu$ [MeV]', horizontalalignment='right', x=1)\n",
    "axes[1,0].set_ylabel( r'$\\Delta\\sigma/\\sigma$', horizontalalignment='right', y=1)\n",
    "axes[1,0].set(ylim=[-1.2,1.2])\n",
    "axes[0,1].set_ylabel( r'$\\langle E_e\\rangle$ [MeV]', horizontalalignment='right', y=1)\n",
    "axes[1,1].set_xlabel( r'$E_\\nu$ [MeV]', horizontalalignment='right', x=1)\n",
    "axes[1,1].set_ylabel( r'$\\Delta\\langle E\\rangle/\\langle E\\rangle$', horizontalalignment='right', y=1)\n",
    "\n",
    "axes[1,1].set(ylim=[-5,15])\n",
    "\n",
    "leg = axes[0,0].legend()\n",
    "fig.suptitle(r'Inv. $\\beta$ Models from Strumia and Vissani, PLB 564, 2003')\n",
    "fig.subplots_adjust(left=0.025, right=0.975)"
   ]
  }
 ],
 "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.11.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}