{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "60b03c4d-6663-43da-9727-7a300d1e9b46",
   "metadata": {},
   "source": [
    "# Interaction Contributions\n",
    "\n",
    "Loop through the set of standard IceCube interactions and plot the contributions to the signal from each under various flavor transformation scenarios and hierarchies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40c44d7e-9206-45d7-a946-f28d7ea8b53c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from asteria.simulation import Simulation\n",
    "from asteria import set_rcparams\n",
    "from asteria import interactions\n",
    "\n",
    "from snewpy.neutrino import Flavor\n",
    "from scipy.interpolate import PchipInterpolator\n",
    "import astropy.units as u\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "\n",
    "set_rcparams(verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91d41f5c-c7e4-4608-9de4-cab7aaf87c1a",
   "metadata": {},
   "source": [
    "## In-Ice Interactions\n",
    "\n",
    "See details in [R. Abbasi+, A&A 535:A109, 2011](https://arxiv.org/abs/1108.0171). We'll plot the interactions in reverse order of their relative importance.\n",
    "\n",
    "1. $^{18}\\mathrm{O}$ CC interaction.\n",
    "2. $^{16}\\mathrm{O}$ NC interaction.\n",
    "3. $^{16}\\mathrm{O}$ CC interaction.\n",
    "4. Electron elastic scattering.\n",
    "5. Inverse beta decay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "572fd112-839f-440d-b351-0c4ca795157b",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = {'name': 'Nakazato_2013',\n",
    "         'param':{\n",
    "             'progenitor_mass': 13 * u.Msun,\n",
    "             'revival_time': 300 * u.ms,\n",
    "             'metallicity': 0.02,\n",
    "             'eos': 'shen'}\n",
    "         }\n",
    "sims = []\n",
    "xform = 'NoTransformation'\n",
    "nmo = 'Normal'\n",
    "\n",
    "#- Loop through individual interactions and compute the signal from each.\n",
    "processes = (interactions.Oxygen18,\n",
    "             interactions.Oxygen16NC,\n",
    "             interactions.Oxygen16CC,\n",
    "             interactions.ElectronScatter,\n",
    "             interactions.InvBetaPar)\n",
    "\n",
    "procnames = ('$^{18}\\mathrm{O}$ CC',\n",
    "             '$^{16}\\mathrm{O}$ NC',\n",
    "             '$^{16}\\mathrm{O}$ CC',\n",
    "             'Electron elastic scattering',\n",
    "             'Inverse beta decay')\n",
    "\n",
    "for proc in processes:\n",
    "    sim = Simulation(model=model,\n",
    "                     interactions=[proc],\n",
    "                     distance=10 * u.kpc,\n",
    "                     Emin=0*u.MeV, Emax=100*u.MeV, dE=1*u.MeV,\n",
    "                     tmin=-10*u.s,  tmax=20*u.s, dt=1*u.ms,\n",
    "                     mixing_scheme=xform, hierarchy=nmo)\n",
    "\n",
    "    sim.run()\n",
    "    sims.append(sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a3d7e14-8e2e-4697-888e-83c43d24c248",
   "metadata": {},
   "source": [
    "## Plot the Contributions of Each Interaction\n",
    "\n",
    "Make a stacked plot fo the total signal as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf7420b7-afd2-4448-81a8-e9b35d9851c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "binsize = 500*u.ms\n",
    "bg = None\n",
    "\n",
    "signal_hits = []\n",
    "for sim in sims:\n",
    "    print(sim.metadata['interactions'])\n",
    "\n",
    "    t, dmu = sim.detector_hits(binsize)\n",
    "    signal_hits.append(dmu)\n",
    "    \n",
    "    if bg is None:\n",
    "        bg = sim.detector.i3_bg(binsize, size=dmu.size) + sim.detector.dc_bg(binsize, size=dmu.size)\n",
    "    hits = (dmu + bg)\n",
    "    print(np.sum(dmu), np.sum(bg))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "763d1ef5-ad99-4a3c-9ab2-46544acb2f72",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,5), tight_layout=True)\n",
    "ax.stackplot(t, np.vstack(signal_hits), step='post', labels=procnames)\n",
    "ax.set(xlim=(0,20),\n",
    "       xlabel='time [s]',\n",
    "       ylim=(1,1e6),\n",
    "       yscale='log',\n",
    "       ylabel='signal hits')\n",
    "ax.legend(fontsize=12);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bc7fc05-4374-43a0-adf2-0183f5f9106e",
   "metadata": {},
   "source": [
    "## Divide Hits into Phases of the CCSN\n",
    "\n",
    "Plot hits showing the distinct phases of the explosion.\n",
    "\n",
    "Consider also different oscillation scenarios and neutrino mass hierarchies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec4d7256-a082-4a97-9a80-b0f0f70d950d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Time limits for phases:\n",
    "limits = [\n",
    "    (-0.025, 0.100),   # Deleptonization, <0.1 s\n",
    "    (0.1, 0.6),        # Accretion, 0.1 - 0.6 s\n",
    "    (1, 10),         # Cooling, 0.6 - 10 s\n",
    "]\n",
    "\n",
    "scale = 1e4\n",
    "binnings = [4e-3, 10e-3, 0.25] * u.s"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf27ba52-3afb-4839-b073-6f1549aeb5ef",
   "metadata": {},
   "source": [
    "### No Flavor Transformations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9537cb4-52b3-4324-bcee-70d4d5a5c56e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,3, figsize = (17,5))\n",
    "\n",
    "bbox_style = {'boxstyle': 'round', 'edgecolor': 'gray', 'facecolor': 'white', 'alpha': 0.75}\n",
    "\n",
    "for ax, binsize, xlim in zip(axes, binnings, limits):\n",
    "    bg = None\n",
    "    signal_hits = []\n",
    "    \n",
    "    # Simulations iterate by mixing scheme (see cell 2)\n",
    "    for sim in sims:\n",
    "        label='x'\n",
    "        color='k'\n",
    "    # for sim, label, color in zip(sims, osc_labels, colors):\n",
    "\n",
    "        # Generate Signal hits\n",
    "        t, dmu = sim.detector_hits(binsize)\n",
    "\n",
    "        signal_hits.append(dmu / scale)\n",
    "        # ax.step(t, hits, label=label, c=color)\n",
    "\n",
    "        # Add 1-sigma band around the expected hits, assuming Poisson uncertainties.\n",
    "        # hits_up = ((dmu + bg) + np.sqrt(dmu + bg))/scale\n",
    "        # hits_lo = ((dmu + bg) - np.sqrt(dmu + bg))/scale\n",
    "        # ax.fill_between(t, hits_lo, hits_up, step='pre', color=color, alpha=0.25)\n",
    "\n",
    "    ax.stackplot(t, np.vstack(signal_hits), step='post', labels=procnames)\n",
    "    ax.set(ylim=(0, None))\n",
    "    \n",
    "    # Normalized to single dom rate in Hz\n",
    "    # ax.step(t, bg/5160/binsize.to(u.s).value, label='Background', c='k', alpha=0.75)\n",
    "    ax.set(xlim=xlim)\n",
    "    if binsize <= 100 * u.ms:\n",
    "        scaled_binsize = binsize.to(u.ms)\n",
    "        annotation = f'{scaled_binsize.value} {scaled_binsize.unit} bins'\n",
    "    else:\n",
    "       annotation = f'{binsize.value} {binsize.unit} bins'\n",
    "    ax.text(0.05, 0.925, annotation, bbox=bbox_style, horizontalalignment='left',\n",
    "            verticalalignment='center', transform = ax.transAxes, fontsize=16)\n",
    "\n",
    "# Plot background\n",
    "axes[0].set_title('Deleptonization')\n",
    "axes[0].set_ylabel(r'signal hits [$\\times10^4$]', horizontalalignment='right', y=1.0)\n",
    "axes[0].set_ylim(0, 0.9)\n",
    "axes[1].set_title('Accretion')\n",
    "axes[1].set_ylim(0, 1.5)\n",
    "axes[2].set_title('PNS Cooling')\n",
    "axes[2].set_xlabel(r'$t-t_\\mathrm{bounce}$ [s]')\n",
    "axes[2].set_ylim(0, 10)\n",
    "axes[2].legend(loc='upper right', ncol=1, fontsize = 12);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ae0f60c-dead-4743-842b-41a14d9d3236",
   "metadata": {},
   "source": [
    "### Adiabatic MSW Effect + Normal Ordering\n",
    "\n",
    "Run a new simulation and regenerate the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "555b055b-df25-4480-babe-ebee464aacd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "xform = 'AdiabaticMSW'\n",
    "nmo = 'Normal'\n",
    "\n",
    "sims = []\n",
    "for proc in processes:\n",
    "    sim = Simulation(model=model,\n",
    "                     interactions=[proc],\n",
    "                     distance=10 * u.kpc,\n",
    "                     Emin=0*u.MeV, Emax=100*u.MeV, dE=1*u.MeV,\n",
    "                     tmin=-10*u.s,  tmax=20*u.s, dt=1*u.ms,\n",
    "                     mixing_scheme=xform, hierarchy=nmo)\n",
    "\n",
    "    sim.run()\n",
    "    sims.append(sim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb3c5a23-2b7d-4ec1-90ec-04553dd9acd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,3, figsize = (17,5))\n",
    "\n",
    "bbox_style = {'boxstyle': 'round', 'edgecolor': 'gray', 'facecolor': 'white', 'alpha': 0.75}\n",
    "\n",
    "for ax, binsize, xlim in zip(axes, binnings, limits):\n",
    "    bg = None\n",
    "    signal_hits = []\n",
    "    \n",
    "    # Simulations iterate by mixing scheme (see cell 2)\n",
    "    for sim in sims:\n",
    "        label='x'\n",
    "        color='k'\n",
    "    # for sim, label, color in zip(sims, osc_labels, colors):\n",
    "\n",
    "        # Generate Signal hits\n",
    "        t, dmu = sim.detector_hits(binsize)\n",
    "\n",
    "        signal_hits.append(dmu / scale)\n",
    "        # ax.step(t, hits, label=label, c=color)\n",
    "\n",
    "        # Add 1-sigma band around the expected hits, assuming Poisson uncertainties.\n",
    "        # hits_up = ((dmu + bg) + np.sqrt(dmu + bg))/scale\n",
    "        # hits_lo = ((dmu + bg) - np.sqrt(dmu + bg))/scale\n",
    "        # ax.fill_between(t, hits_lo, hits_up, step='pre', color=color, alpha=0.25)\n",
    "\n",
    "    ax.stackplot(t, np.vstack(signal_hits), step='post', labels=procnames)\n",
    "    ax.set(ylim=(0, None))\n",
    "    \n",
    "    # Normalized to single dom rate in Hz\n",
    "    # ax.step(t, bg/5160/binsize.to(u.s).value, label='Background', c='k', alpha=0.75)\n",
    "    ax.set(xlim=xlim)\n",
    "    if binsize <= 100 * u.ms:\n",
    "        scaled_binsize = binsize.to(u.ms)\n",
    "        annotation = f'{scaled_binsize.value} {scaled_binsize.unit} bins'\n",
    "    else:\n",
    "       annotation = f'{binsize.value} {binsize.unit} bins'\n",
    "    ax.text(0.05, 0.925, annotation, bbox=bbox_style, horizontalalignment='left',\n",
    "            verticalalignment='center', transform = ax.transAxes, fontsize=16)\n",
    "\n",
    "# Plot background\n",
    "axes[0].set_title('Deleptonization')\n",
    "axes[0].set_ylabel(r'signal hits [$\\times10^4$]', horizontalalignment='right', y=1.0)\n",
    "axes[0].set_ylim(0, 0.9)\n",
    "axes[1].set_title('Accretion')\n",
    "axes[1].set_ylim(0, 1.5)\n",
    "axes[2].set_title('PNS Cooling')\n",
    "axes[2].set_xlabel(r'$t-t_\\mathrm{bounce}$ [s]')\n",
    "axes[2].set_ylim(0, 10)\n",
    "axes[2].legend(loc='upper right', ncol=1, fontsize = 12);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa08ffe0-5793-46cf-846e-21eb29c80aa7",
   "metadata": {},
   "source": [
    "### Adiabatic MSW Effect + Inverted Ordering\n",
    "\n",
    "Run a new simulation and regenerate the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f851857-31f5-4f9a-afc6-2cb0a6dfcb8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "xform = 'AdiabaticMSW'\n",
    "nmo = 'Inverted'\n",
    "\n",
    "sims = []\n",
    "for proc in processes:\n",
    "    sim = Simulation(model=model,\n",
    "                     interactions=[proc],\n",
    "                     distance=10 * u.kpc,\n",
    "                     Emin=0*u.MeV, Emax=100*u.MeV, dE=1*u.MeV,\n",
    "                     tmin=-10*u.s,  tmax=20*u.s, dt=1*u.ms,\n",
    "                     mixing_scheme=xform, hierarchy=nmo)\n",
    "\n",
    "    sim.run()\n",
    "    sims.append(sim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73655883-429f-4134-87c7-6ed1b46d95c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1,3, figsize = (17,5))\n",
    "\n",
    "bbox_style = {'boxstyle': 'round', 'edgecolor': 'gray', 'facecolor': 'white', 'alpha': 0.75}\n",
    "\n",
    "for ax, binsize, xlim in zip(axes, binnings, limits):\n",
    "    bg = None\n",
    "    signal_hits = []\n",
    "    \n",
    "    # Simulations iterate by mixing scheme (see cell 2)\n",
    "    for sim in sims:\n",
    "        label='x'\n",
    "        color='k'\n",
    "    # for sim, label, color in zip(sims, osc_labels, colors):\n",
    "\n",
    "        # Generate Signal hits\n",
    "        t, dmu = sim.detector_hits(binsize)\n",
    "\n",
    "        signal_hits.append(dmu / scale)\n",
    "        # ax.step(t, hits, label=label, c=color)\n",
    "\n",
    "        # Add 1-sigma band around the expected hits, assuming Poisson uncertainties.\n",
    "        # hits_up = ((dmu + bg) + np.sqrt(dmu + bg))/scale\n",
    "        # hits_lo = ((dmu + bg) - np.sqrt(dmu + bg))/scale\n",
    "        # ax.fill_between(t, hits_lo, hits_up, step='pre', color=color, alpha=0.25)\n",
    "\n",
    "    ax.stackplot(t, np.vstack(signal_hits), step='post', labels=procnames)\n",
    "    ax.set(ylim=(0, None))\n",
    "    \n",
    "    # Normalized to single dom rate in Hz\n",
    "    # ax.step(t, bg/5160/binsize.to(u.s).value, label='Background', c='k', alpha=0.75)\n",
    "    ax.set(xlim=xlim)\n",
    "    if binsize <= 100 * u.ms:\n",
    "        scaled_binsize = binsize.to(u.ms)\n",
    "        annotation = f'{scaled_binsize.value} {scaled_binsize.unit} bins'\n",
    "    else:\n",
    "       annotation = f'{binsize.value} {binsize.unit} bins'\n",
    "    ax.text(0.05, 0.925, annotation, bbox=bbox_style, horizontalalignment='left',\n",
    "            verticalalignment='center', transform = ax.transAxes, fontsize=16)\n",
    "\n",
    "# Plot background\n",
    "axes[0].set_title('Deleptonization')\n",
    "axes[0].set_ylabel(r'signal hits [$\\times10^4$]', horizontalalignment='right', y=1.0)\n",
    "axes[0].set_ylim(0, 0.9)\n",
    "axes[1].set_title('Accretion')\n",
    "axes[1].set_ylim(0, 1.5)\n",
    "axes[2].set_title('PNS Cooling')\n",
    "axes[2].set_xlabel(r'$t-t_\\mathrm{bounce}$ [s]')\n",
    "axes[2].set_ylim(0, 10)\n",
    "axes[2].legend(loc='upper right', ncol=1, fontsize = 12);"
   ]
  }
 ],
 "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": 5
}