{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "115a4e1c-1ab6-40f9-ba0d-1be45e0bbd4d",
   "metadata": {},
   "source": [
    "# Detector Hits in IceCube-Gen2\n",
    "\n",
    "Demonstrate detector hits in the IC86 + Gen2 configuration of IceCube.\n",
    "\n",
    "This includes both mDOMs and WLS modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2eb397f-92c2-40d2-932f-5816b85cf46b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from snewpy.neutrino import Flavor\n",
    "\n",
    "from asteria.simulation import Simulation\n",
    "from asteria import set_rcparams\n",
    "from asteria import interactions\n",
    "\n",
    "from importlib.resources import files\n",
    "\n",
    "set_rcparams(verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7af47266-c61c-45a7-b572-b865a1092864",
   "metadata": {},
   "source": [
    "## Configure an mDOM Gen2 Simulation\n",
    "\n",
    "Set up a progenitor 10 kpc from Earth and generate simulated hits.\n",
    "\n",
    "To simulate the Gen2 detector, add the argument `detector_scope='Gen2'` when initializing the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dff7931-b55f-4563-ab3b-12cf5e78e2bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = {'name': 'Nakazato_2013',\n",
    "         'param':{\n",
    "             'progenitor_mass': 13 * u.Msun,\n",
    "             'revival_time': 300 * u.ms,\n",
    "             'metallicity': 0.004,\n",
    "             'eos': 'shen'}\n",
    "         }\n",
    "\n",
    "sim = Simulation(model=model,\n",
    "                 distance=10 * u.kpc, \n",
    "                 Emin=0*u.MeV, Emax=100*u.MeV, dE=1*u.MeV,\n",
    "                 tmin=-1*u.s, tmax=1*u.s, dt=1*u.ms,\n",
    "                 detector_scope='Gen2')\n",
    "sim.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3052237c-8160-44b7-ac0e-bb2c4e986ee2",
   "metadata": {},
   "source": [
    "## Plot the Energy Deposit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8252ebd0-52b8-4504-915e-3f36189f5419",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(8,6), tight_layout=True)\n",
    "\n",
    "for flavor in sim.flavors:\n",
    "    ax.plot(sim.time, sim.E_per_V[flavor], label=flavor.to_tex())\n",
    "ax.legend()\n",
    "ax.set(xlabel=r'$t-t_\\mathrm{bounce}$ [s]',\n",
    "       ylabel='energy deposit [Mev m$^{-3}$])',\n",
    "       xlim=(-0.15, 0.75));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe116039-a071-467a-8bd6-06aa2a812468",
   "metadata": {},
   "source": [
    "## Plot Detector Response\n",
    "\n",
    "### Expected Signal from Each Subdetector\n",
    "\n",
    "Set a time resolution `dt`. Using the `sim.detector_signal()` function we can read out the detector signal for each subdetector class. Separately plot hits from the main IceCube strings and DeepCore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d90d09cf-5ecc-475c-bbc5-178c95e4e2c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 2 * u.ms\n",
    "t, sim_i3 = sim.detector_signal(subdetector='i3', dt=dt)\n",
    "t, sim_dc = sim.detector_signal(subdetector='dc', dt=dt)\n",
    "t, sim_md = sim.detector_signal(subdetector='md', dt=dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52d3c595-17e0-4df0-9746-af9d47e77a35",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(8,6), tight_layout=True)\n",
    "ax.plot(t, sim_i3, label='IceCube DOM')\n",
    "ax.plot(t, sim_dc, label='HQE DOM (DeepCore)')\n",
    "ax.plot(t, sim_md, label='mDOM (Gen2)')\n",
    "ax.legend(fontsize=14)\n",
    "ax.set(xlabel=r'$t-t_\\mathrm{bounce}$ [s]',\n",
    "       ylabel=f'detector hits',\n",
    "       xlim=(-0.15, 0.75),\n",
    "       ylim=(0,10000));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5854b46f-9889-48fd-bdef-d69678b5c147",
   "metadata": {},
   "source": [
    "### Generated Hits from Signal Only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c29ae3a0-b978-4ada-8580-a0a167d89dba",
   "metadata": {},
   "outputs": [],
   "source": [
    "t, hits_i3 = sim.detector_hits(subdetector='i3', dt=dt)\n",
    "t, hits_dc = sim.detector_hits(subdetector='dc', dt=dt)\n",
    "t, hits_md = sim.detector_signal(subdetector='md', dt=dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb6a1409-3da3-4071-ab48-e9a8d5bfdfef",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(8,6), tight_layout=True)\n",
    "ax.plot(t, hits_i3, label='IceCube DOM')\n",
    "ax.plot(t, hits_dc, label='HQE DOM (DeepCore)')\n",
    "ax.plot(t, hits_md, label='mDOM (Gen2)')\n",
    "\n",
    "ax.legend(fontsize=14)\n",
    "ax.set(xlabel=r'$t-t_\\mathrm{bounce}$ [s]',\n",
    "       ylabel=f'detector hits',\n",
    "       xlim=(-0.15, 0.75),\n",
    "       ylim=(0,10000));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f3f7030-fb90-4a89-b98f-d0d0e9787446",
   "metadata": {},
   "source": [
    "### Generated Hits from Signal + Background\n",
    "\n",
    "Separately compute the background hits and signal from each subdetector and add them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2f39f58-504e-42ca-8f0f-8a72f7dd8eb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_i3 = sim.detector.i3_bg(dt, size=hits_i3.size)\n",
    "bkg_dc = sim.detector.dc_bg(dt, size=hits_dc.size)\n",
    "bkg_md = sim.detector.md_bg(dt, size=hits_md.size)\n",
    "\n",
    "bkg = bkg_i3 + bkg_dc + bkg_md\n",
    "hits = hits_i3 + hits_dc + hits_md"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4df03ef-e5f6-40a4-a2cd-1f0c94c51d46",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(8,6), tight_layout=True)\n",
    "ax.plot(t, hits_i3 + bkg_i3, label='IceCube DOM')\n",
    "ax.plot(t, hits_dc + bkg_dc, label='HQE DOM (DeepCore)')\n",
    "ax.plot(t, hits_md + bkg_md, label='mDOM (Gen2)')\n",
    "ax.plot(t, hits + bkg, label='Total hits')\n",
    "\n",
    "ax.legend(loc='center right', fontsize=14)\n",
    "ax.set(xlabel=r'$t-t_\\mathrm{bounce}$ [s]',\n",
    "       ylabel=f'detector hits',\n",
    "       xlim=(-0.15, 0.75),\n",
    "       ylim=(0,65000)\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e03521b3-ec42-4bd9-87ce-d289bd8a1591",
   "metadata": {},
   "source": [
    "## Add WLS Modules\n",
    "\n",
    "In the Gen2 detector scope, using the `add_wls=True` argument will add wavelength-shifting modules to the Gen2 detector response.\n",
    "\n",
    "The hits can be accessed as the `ws` subdetector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cebd2b2-7533-4d1f-a074-d5299a2f00df",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_ws = Simulation(model=model,\n",
    "                    distance=10 * u.kpc, \n",
    "                    Emin=0*u.MeV, Emax=100*u.MeV, dE=1*u.MeV,\n",
    "                    tmin=-1*u.s, tmax=1*u.s, dt=1*u.ms,\n",
    "                    detector_scope='Gen2',\n",
    "                    add_wls=True)\n",
    "sim_ws.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c8de85b-f01a-41f1-99cd-7cd79a296ad0",
   "metadata": {},
   "outputs": [],
   "source": [
    "t, hits_i3 = sim_ws.detector_hits(subdetector='i3', dt=dt)\n",
    "t, hits_dc = sim_ws.detector_hits(subdetector='dc', dt=dt)\n",
    "t, hits_md = sim_ws.detector_signal(subdetector='md', dt=dt)\n",
    "t, hits_ws = sim_ws.detector_signal(subdetector='ws', dt=dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00ac3f74-787d-4098-bd3b-24f2d394eb39",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(8,6), tight_layout=True)\n",
    "ax.plot(t, hits_i3, label='IceCube DOM')\n",
    "ax.plot(t, hits_dc, label='HQE DOM (DeepCore)')\n",
    "ax.plot(t, hits_md, label='mDOM (Gen2)')\n",
    "ax.plot(t, hits_ws, label='mDOM+WLS (Gen2)')\n",
    "\n",
    "ax.legend(fontsize=14)\n",
    "ax.set(xlabel=r'$t-t_\\mathrm{bounce}$ [s]',\n",
    "       ylabel=f'detector hits',\n",
    "       xlim=(-0.15, 0.75),\n",
    "       ylim=(0,10000));"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "286dd290-371f-4672-92d5-0a54eb7fb77b",
   "metadata": {},
   "source": [
    "### Plot the Sum of Signal + Background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c2d27df-8173-45c6-84ad-5c63214a681c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_i3 = sim.detector.i3_bg(dt, size=hits_i3.size)\n",
    "bkg_dc = sim.detector.dc_bg(dt, size=hits_dc.size)\n",
    "bkg_md = sim.detector.md_bg(dt, size=hits_md.size)\n",
    "bkg_ws = sim.detector.ws_bg(dt, size=hits_ws.size)\n",
    "\n",
    "bkg = bkg_i3 + bkg_dc + bkg_md + bkg_ws\n",
    "hits = hits_i3 + hits_dc + hits_md + hits_ws"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a6fe116-4718-46b8-8991-a2d6cb0db220",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(8,6), tight_layout=True)\n",
    "ax.plot(t, hits_i3 + bkg_i3, label='IceCube DOM')\n",
    "ax.plot(t, hits_dc + bkg_dc, label='HQE DOM (DeepCore)')\n",
    "ax.plot(t, hits_md + bkg_md, label='mDOM (Gen2)')\n",
    "ax.plot(t, hits_ws + bkg_ws, label='mDOM+WLS (Gen2)')\n",
    "ax.plot(t, hits + bkg, label='Total hits')\n",
    "\n",
    "ax.legend(loc='center', ncol=2, fontsize=14)\n",
    "ax.set(xlabel=r'$t-t_\\mathrm{bounce}$ [s]',\n",
    "       ylabel=f'detector hits',\n",
    "       xlim=(-0.15, 0.75),\n",
    "       ylim=(0,75000)\n",
    ");"
   ]
  }
 ],
 "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
}