{ "cells": [ { "cell_type": "markdown", "id": "1baaab8f-1a4a-4616-83f9-6821b0a811bc", "metadata": {}, "source": [ "# Analytic Model Initialization\n", "\n", "Generate a model file based on the old LLNL 1D simulation tuned to SN 1987A." ] }, { "cell_type": "code", "execution_count": null, "id": "41490e7b", "metadata": {}, "outputs": [], "source": [ "import os\n", "\n", "from astropy.table import Table\n", "import astropy.units as u\n", "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "import numpy as np\n", "\n", "from snewpy import model_path\n", "from snewpy.neutrino import Flavor\n", "from snewpy.models.ccsn import Analytic3Species\n", "\n", "from asteria.simulation import Simulation\n", "from asteria import set_rcparams\n", "\n", "set_rcparams()" ] }, { "cell_type": "code", "execution_count": null, "id": "2790c4a0", "metadata": {}, "outputs": [], "source": [ "filename = \"AnalyticFluenceExample.dat\"\n", "model_folder = f\"{model_path}/AnalyticFluence/\"\n", "\n", "if not os.path.exists(model_folder):\n", " os.makedirs(model_folder)\n", "file_path = os.path.join(model_folder, filename)" ] }, { "cell_type": "markdown", "id": "1ae4d76e", "metadata": {}, "source": [ "## Creating a SN model file modeled after the Livermore model\n", "\n", "This code was taken from [SNEWS2/snewpy](https://github.com/SNEWS2/snewpy) repository, from the [AnalyticFluence.ipynb](https://github.com/SNEWS2/snewpy/blob/main/doc/nb/AnalyticFluence.ipynb) example notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "52c8193e", "metadata": {}, "outputs": [], "source": [ "# These numbers _almost_ reproduce the Livermore model included in the SNOwGLoBES repository.\n", "# They are obtained by calculating the total L, and from the livermore.dat\n", "# fluence file (which is modelled after a 10kpc supernova).\n", "total_energy = (5.478e+52, 5.485e+52, 4 * 5.55e+52)\n", "mean_energy = (11.5081, 15.4678, 21.0690)\n", "rms_or_pinch = \"rms\"\n", "rms_energy = (12.8788, 17.8360, 24.3913)\n", "\n", "# Make an astropy table with two times, 0s and 1s, with constant neutrino properties\n", "table = Table()\n", "table['TIME'] = np.linspace(0,1,2)\n", "table['L_NU_E'] = np.linspace(1,1,2)*total_energy[0]\n", "table['L_NU_E_BAR'] = np.linspace(1,1,2)*total_energy[1]\n", "table['L_NU_X'] = np.linspace(1,1,2)*total_energy[2]/4. #Note, L_NU_X is set to 1/4 of the total NU_X energy\n", " \n", "table['E_NU_E'] = np.linspace(1,1,2)*mean_energy[0]\n", "table['E_NU_E_BAR'] = np.linspace(1,1,2)*mean_energy[1]\n", "table['E_NU_X'] = np.linspace(1,1,2)*mean_energy[2]\n", "\n", "if rms_or_pinch == \"rms\":\n", " table['RMS_NU_E'] = np.linspace(1,1,2)*rms_energy[0]\n", " table['RMS_NU_E_BAR'] = np.linspace(1,1,2)*rms_energy[1]\n", " table['RMS_NU_X'] = np.linspace(1,1,2)*rms_energy[2]\n", " table['ALPHA_NU_E'] = (2.0 * table['E_NU_E'] ** 2 - table['RMS_NU_E'] ** 2) / (\n", " table['RMS_NU_E'] ** 2 - table['E_NU_E'] ** 2)\n", " table['ALPHA_NU_E_BAR'] = (2.0 * table['E_NU_E_BAR'] ** 2 - table['RMS_NU_E_BAR'] ** 2) / (\n", " table['RMS_NU_E_BAR'] ** 2 - table['E_NU_E_BAR'] ** 2)\n", " table['ALPHA_NU_X'] = (2.0 * table['E_NU_X'] ** 2 - table['RMS_NU_X'] ** 2) / (\n", " table['RMS_NU_X'] ** 2 - table['E_NU_X'] ** 2)\n", "elif rms_or_pinch == \"pinch\":\n", " table['ALPHA_NU_E'] = np.linspace(1,1,2)*pinch_values[0]\n", " table['ALPHA_NU_E_BAR'] = np.linspace(1,1,2)*pinch_values[1]\n", " table['ALPHA_NU_X'] = np.linspace(1,1,2)*pinch_values[2]\n", " table['RMS_NU_E'] = np.sqrt((2.0 + table['ALPHA_NU_E'])/(1.0 + table['ALPHA_NU_E'])*table['E_NU_E']**2)\n", " table['RMS_NU_E_BAR'] = np.sqrt((2.0 + table['ALPHA_NU_E_BAR'])/(1.0 + table['ALPHA_NU_E_BAR'])*table['E_NU_E_BAR']**2)\n", " table['RMS_NU_X'] = np.sqrt((2.0 + table['ALPHA_NU_X'])/(1.0 + table['ALPHA_NU_X'])*table['E_NU_X']**2 )\n", "else:\n", " print(\"incorrect second moment method: rms or pinch\")\n", "\n", "table.write(file_path,format='ascii',overwrite=True)" ] }, { "cell_type": "markdown", "id": "d967c2e5", "metadata": {}, "source": [ "## ASTERIA Simulation" ] }, { "cell_type": "code", "execution_count": null, "id": "ac96309b", "metadata": {}, "outputs": [], "source": [ "# SNEWPY model dictionary, the format must match the below example for analytic models\n", "model = {\n", " 'name': 'Analytic3Species',\n", " 'param': {\n", " 'filename': file_path\n", " }\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=10*u.s, dt=1*u.ms,\n", " mixing_scheme='AdiabaticMSW',\n", " hierarchy='normal')\n", "sim.run()" ] }, { "cell_type": "markdown", "id": "3e429b3f-6e5f-4701-9ddc-97369f7bfd68", "metadata": {}, "source": [ "### Plot Output" ] }, { "cell_type": "code", "execution_count": null, "id": "fc2640ef", "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, figsize = (6,6))\n", "dt = 0.5 * u.s\n", "scale = 1e4\n", "\n", "sim.rebin_result(dt)\n", "t, hits = sim.detector_signal(dt)\n", "bg = sim.detector.i3_bg(dt, size=hits.size) + sim.detector.dc_bg(dt, size=hits.size)\n", "\n", "\n", "ax.step(t, hits+bg, where='post', lw=2, )\n", "ax.set(xlim=(-1, 5));\n", "ax.set_xlabel(r'$t-t_\\mathrm{bounce}$ [s]', ha='right', x=1.0)\n", "ax.set_ylabel(fr'total signal [{dt} bins, $\\times 10^{{{int(np.log10(scale))}}}$]', ha='right', y=1.0)" ] }, { "cell_type": "code", "execution_count": null, "id": "106f8205", "metadata": {}, "outputs": [], "source": [] } ], "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 }