{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# SNEWPY Demo for PyHEP 2022\n", "\n", "## Getting Started\n", "\n", "### Resources\n", "\n", "* [snewpy repository on GitHub](https://github.com/SNEWS2/snewpy)\n", "* [snewpy documentation](https://snewpy.rtfd.io)\n", "* snewpy papers:\n", " * Software review: [![DOI](https://joss.theoj.org/papers/10.21105/joss.03772/status.svg)](https://doi.org/10.21105/joss.03772)\n", " * Underlying physics: [arXiv:2109.08188](https://arxiv.org/abs/2109.08188) or [DOI:10.3847/1538-4357/ac350f](https://dx.doi.org/10.3847/1538-4357/ac350f)\n", "\n", "### Installation\n", "Click the following button to open this notebook in Binder. Everything should just work.\n", "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/JostMigenda/PyHEP2022/main?labpath=snewpy_demo.ipynb)\n", "\n", "To run this notebook locally on your own machine, use `pip install snewpy==1.3b1` to install snewpy and `git clone https://github.com/SNOwGLoBES/snowglobes.git` to download [SNOwGLoBES](https://github.com/SNOwGLoBES/snowglobes), which we will need in the second half of this notebook to estimate event rates in realistic detectors.\n", "\n", "### Get Sample Data Files\n", "\n", "Let’s download data files containing neutrino fluxes predicted by a few different supernova models. Here, we will use the `Nakazato_2013` and `Bollig_2016` models, but plenty of others are available as part of snewpy." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import snewpy\n", "\n", "snewpy.get_models(models=[\"Nakazato_2013\", \"Bollig_2016\"])\n", "\n", "# You can also run `snewpy.get_models()` without arguments,\n", "# to interactively select which model files to download\n", "#snewpy.get_models()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting Different Supernova Models\n", "\n", "The Nakazato and Bollig model families use very different input file formats, different time binning, etc. Thankfully, snewpy takes care of all of that and provides us with a unified interface.\n", "\n", "Let’s take a look at one model from each model family!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from snewpy.models.ccsn import Nakazato_2013, Bollig_2016\n", "\n", "nakazato = Nakazato_2013('SNEWPY_models/Nakazato_2013/nakazato-shen-z0.004-t_rev100ms-s20.0.fits')\n", "bollig = Bollig_2016('SNEWPY_models/Bollig_2016/s27.0c') # This model has one file per flavor. Use common prefix, not full filename.\n", "\n", "nakazato" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "\n", "from snewpy.neutrino import Flavor\n", "\n", "mpl.rc('font', size=16)\n", "\n", "# Plot luminosity of both models\n", "fig, ax = plt.subplots(1, figsize=(10, 6))\n", "\n", "for flavor in Flavor:\n", " ax.plot(nakazato.time, nakazato.luminosity[flavor]/1e51, # Report luminosity in units foe/s\n", " label=flavor.to_tex() + ' (Nakazato)',\n", " color='C0' if flavor.is_electron else 'C2',\n", " ls='-' if flavor.is_neutrino else '--',\n", " lw=2)\n", "\n", "# for flavor in Flavor:\n", "# ax.plot(bollig.time, bollig.luminosity[flavor]/1e51, # Report luminosity in units foe/s\n", "# label=flavor.to_tex() + ' (Bollig)',\n", "# color='C1' if flavor.is_electron else 'C3',\n", "# ls='-' if flavor.is_neutrino else '--',\n", "# lw=1)\n", "\n", "ax.set(xlim=(-0.05, 0.5), xlabel=r'$t-t_{\\rm bounce}$ [s]', ylabel=r'luminosity [$10^{51}$ erg s$^{-1}$]')\n", "ax.grid()\n", "ax.legend(loc='upper right', ncol=2, fontsize=18);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Applying Neutrino Flavo(u)r Transformations\n", "\n", "Supernova simulations only provide us with the neutrino fluxes *emitted near the centre of the supernovae*. The fluxes *observed in a detector on Earth* can be very different.\n", "\n", "The module `snewpy.flavor_transformation` contains many different transformations that can happen along the way. Here, we will use the `AdiabaticMSW` transformation, which occurs as the neutrinos exit the star. (In brief: The electron density changes between the centre of the star and the surface, which modifies the effective mass of electron (anti-)neutrinos; whereas muon or tau (anti-)neutrinos are unaffected.)\n", "\n", "Other flavor transformations involving e.g. non-adiabatic MSW, decoherence, neutrino decay or sterile neutrinos are also available." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from snewpy.flavor_transformation import AdiabaticMSW\n", "from snewpy.neutrino import Flavor, MassHierarchy\n", "\n", "from astropy import units as u\n", "import matplotlib as mpl\n", "import numpy as np\n", "\n", "def plot_total_flux(model, xform_nmo, xform_imo):\n", " \"\"\"Plot initial and oscillated neutrino luminosities.\"\"\"\n", " \n", " energies = np.linspace(0,60,121) * u.MeV \n", " d = (10*u.kpc).to('cm').value # distance to SN\n", " \n", " times = model.get_time()\n", " burst_epoch = times <= 0.1*u.s\n", " accretion_epoch = (times > 0.1*u.s) & (times <= 0.5*u.s)\n", " cooling_epoch = (times > 0.5*u.s) & (times <= 10*u.s)\n", " \n", " ilum = {}\n", " tlum_nmo = {}\n", " tlum_imo = {}\n", " \n", " for flavor in Flavor:\n", " ilum[flavor] = np.zeros(len(times))\n", " tlum_nmo[flavor] = np.zeros(len(times))\n", " tlum_imo[flavor] = np.zeros(len(times))\n", "\n", " # Compute the transformed and untransformed flux at each time.\n", " for i, t in enumerate(times):\n", " ispec = model.get_initial_spectra(t, energies)\n", " tspec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)\n", " tspec_imo = model.get_transformed_spectra(t, energies, xform_imo)\n", "\n", " for flavor in Flavor:\n", " for j, E in enumerate(energies):\n", " ispec[flavor][j] /= (4.*np.pi*d**2)\n", " tspec_nmo[flavor][j] /= (4.*np.pi*d**2)\n", " tspec_imo[flavor][j] /= (4.*np.pi*d**2)\n", " \n", " for flavor in Flavor:\n", " ilum[flavor][i] = np.trapz(ispec[flavor].to('1/(erg*s)'), energies.to('erg')).value \n", " tlum_nmo[flavor][i] = np.trapz(tspec_nmo[flavor].to('1/(erg*s)'), energies.to('erg')).value\n", " tlum_imo[flavor][i] = np.trapz(tspec_imo[flavor].to('1/(erg*s)'), energies.to('erg')).value \n", " \n", " # make the figures \n", " fig, axes = plt.subplots(3,3, figsize=(20,12), tight_layout=True)\n", " \n", " smax = [0.,0.,0.]\n", " titles = ['Untransformed', 'Transformed (NMO)', 'Transformed (IMO)']\n", " for i, spec in enumerate([ilum, tlum_nmo, tlum_imo]):\n", " for j, phase in enumerate([burst_epoch, accretion_epoch, cooling_epoch]):\n", " ax = axes[i,j]\n", " timeunits = 'ms' if j==0 else 's'\n", " \n", " for flavor in Flavor:\n", " if i == 0:\n", " smax[j] = np.maximum(smax[j], 1.1*np.max(spec[flavor][phase]))\n", " \n", " ax.plot(times[phase].to(timeunits),\n", " spec[flavor][phase], label=flavor.to_tex(), lw=3,\n", " color='C0' if flavor.is_electron else 'C1',\n", " ls='-' if flavor.is_neutrino else ':')\n", " \n", " ax.set(xlim=(times[phase][0].to(timeunits).value, times[phase][-1].to(timeunits).value),\n", " ylim=(0, smax[j]))\n", " \n", " if j==0:\n", " ax.set(ylabel=r'flux [cm$^{-2}$ s$^{-1}$]')\n", " ax.legend(loc='upper right', ncol=1, fontsize=18)\n", " if j==1:\n", " ax.set(title=titles[i])\n", " if i < 2:\n", " ax.set(xticklabels=[])\n", " else:\n", " ax.set(xlabel='time [{}]'.format(timeunits))\n", " \n", " ax.grid(ls=':')\n", "\n", " return fig\n", "\n", "def plot_spectra(model, xform, t):\n", " \"\"\"Plot initial and oscillated neutrino spectra at a fixed time.\"\"\"\n", "\n", " energies = np.linspace(0,60,121) * u.MeV \n", " d = (10*u.kpc).to('cm').value # distance to SN\n", "\n", " #get the spectra\n", " ispec = model.get_initial_spectra(t, energies) \n", " tspec = model.get_transformed_spectra(t, energies, xform)\n", "\n", " for flavor in Flavor:\n", " for j, E in enumerate(energies):\n", " ispec[flavor][j] /= (4.*np.pi*d**2)\n", " tspec[flavor][j] /= (4.*np.pi*d**2)\n", " \n", " fig, axes = plt.subplots(1, 2, figsize=(18,7), sharey=True, tight_layout=True)\n", "\n", " for i, spec in enumerate([ispec, tspec]):\n", " axes[0].plot(energies, spec[Flavor.NU_E]/1e16, \n", " label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),\n", " color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", " axes[0].plot(energies, spec[Flavor.NU_X]/1e16, \n", " label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),\n", " color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", "\n", " axes[0].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Neutrino spectra at $t = ${:.1f}'.format(t))\n", " axes[0].grid()\n", " axes[0].legend(loc='upper right', ncol=2, fontsize=16)\n", " \n", " axes[1].plot(energies, spec[Flavor.NU_E_BAR]/1e16, \n", " label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),\n", " color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", " axes[1].plot(energies, spec[Flavor.NU_X_BAR]/1e16, \n", " label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),\n", " color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n", "\n", " axes[1].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Antineutrino spectra at $t = ${:.1f}'.format(t)) \n", " axes[1].grid()\n", " axes[1].legend(loc='upper right', ncol=2, fontsize=16) \n", " \n", " ax = axes[0]\n", " ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')\n", " \n", " return fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After defining two helper functions above, let’s use the first one to plot the neutrino flux over time! Focus on the “neutronization burst” in the left panel—this will show most clearly how some of the electron neutrinos are transformed into muon/tau neutrinos:\n", "\n", "(Note: Most supernova simulations don’t distinguish between muon & tau neutrinos, since their interactions are virtually identical; so we simply label them as NU_X here.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xform_nmo = AdiabaticMSW()\n", "xform_imo = AdiabaticMSW(mh=MassHierarchy.INVERTED)\n", "\n", "fig = plot_total_flux(nakazato, xform_nmo, xform_imo)\n", "# fig.savefig('flux_adiabaticmsw.pdf')\n", "# fig.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we’ll use the second helper function to plot the neutrino spectra at a fixed point in time. Note how the total flux of electron (anti-) neutrinos decreases due to the flavour transformation, but their mean energy increases." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig_nmo = plot_spectra(nakazato, xform_nmo, 100*u.ms)\n", "# fig_nmo.show()\n", "# fig_imo = plot_spectra(nakazato, xform_imo, 100*u.ms)\n", "# fig_imo.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate Event Rates in a Neutrino Detector\n", "\n", "Finally, we can use snewpy to calculate event rates for our models in various neutrino detectors.\n", "This requires SNOwGLoBES, a separate software package that contains cross sections for all relevant interaction channels and files describing energy smearing & efficiency for many detectors.\n", "\n", "First, we pick a detector, supernova model and time range:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from snewpy import snowglobes\n", "\n", "SNOwGLoBES_path = \"./snowglobes/\" # where SNOwGLoBES is located\n", "SNEWPY_models_base = \"./SNEWPY_models/\" # directory containing SNEWPY model files\n", "\n", "# set distance in kpc\n", "distance = 10\n", "\n", "# set SNOwGLoBES detector to use\n", "detector = \"wc100kt30prct\" # like Super- or Hyper-Kamiokande\n", "# detector = \"ar40kt\" # like DUNE\n", "\n", "# set SNEWPY model type and filename\n", "modeltype = 'Bollig_2016'\n", "model = 's27.0c'\n", "\n", "# set desired flavor transformation\n", "transformation = 'NoTransformation' #'AdiabaticMSW_IMO'\n", "\n", "# Construct file system path of model file and name of output file\n", "# The output file will be stored in the same directory as the model file.\n", "modelfile = f\"{SNEWPY_models_base}/{modeltype}/{model}\"\n", "outfile = f\"{modeltype}_{model}_{transformation}\"\n", "\n", "# There are three ways to select a time range.\n", "# Option 1 - don't specify tstart and tend, then the whole model is integrated\n", "# tstart = None\n", "# tend = None\n", "\n", "# Option 2 - specify single tstart and tend, this makes 1 fluence file integrated over the window\n", "tstart = 0.0 * u.s\n", "tend = 1.0 * u.s\n", "\n", "# Option 3 = specify sequence of time intervals, one fluence file is made for each interval\n", "# window_tstart = 0.0\n", "# window_tend = 1.0\n", "# window_bins = 10\n", "# tstart = np.linspace(window_tstart, window_tend, window_bins, endpoint=False) * u.s\n", "# tend = tstart + (window_tend - window_tstart) / window_bins * u.s" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we perform the calculation in three steps:\n", "1. `generate_fluence` integrates the flux over the specified time window(s) to get the fluence as a function of energy\n", "2. `simulate` multiplies the fluence, cross section and detector smearing & efficiency to get the number of observed events as a function of energy\n", "3. `collate` collects the results for different channels into a single dictionary" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Preparing fluences ...\")\n", "tarredfile = snowglobes.generate_fluence(modelfile, modeltype, transformation, distance, outfile, tstart, tend)\n", "\n", "print(\"Calculating events ...\")\n", "snowglobes.simulate(SNOwGLoBES_path, tarredfile, detector_input=detector)\n", "\n", "print(\"Collating results ...\")\n", "tables = snowglobes.collate(SNOwGLoBES_path, tarredfile, skip_plots=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That was quick!\n", "\n", "Finally, let’s plot the event spectra observed in this detector:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "key = f\"Collated_{outfile}_{detector}_events_smeared_weighted.dat\"\n", "# SNOwGLoBES also provides “unsmeared” events, without accounting for detector effects.\n", "# Note: If you use Option 3 above (multiple time bins), `key` also includes the number of the time bin.\n", "\n", "tables[key]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "channels = tables[key]['header'].split()\n", "energies = tables[key]['data'][0] * 1e3 # convert GeV to MeV\n", "events = {}\n", "for i, channel in enumerate(channels):\n", " if i==0: continue # skip \"Energy\" array\n", " events[channel] = tables[key]['data'][i]\n", "events['total'] = sum(tables[key]['data'][1:])\n", "\n", "for channel in ['total'] + channels[1:]:\n", " plt.plot(energies, events[channel], label=channel)\n", "\n", "plt.xlabel(\"Energy [MeV]\")\n", "plt.ylabel(\"Events per energy bin\")\n", "plt.yscale(\"log\")\n", "plt.xlim(0, 70)\n", "plt.ylim(bottom=1)\n", "plt.legend()\n", "plt.savefig(f'plot_{outfile}_{detector}.pdf')\n", "\n", "print(f\"Total Events: {sum(events['total']):.3f}\")\n", "print(f\"Mean Detected Energy: {sum(energies * events['total']) / sum(events['total']):.3f} MeV\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary and Q&A\n", "\n", "* snewpy provides …\n", " * … a unified way to access 100s of SN models\n", " * … many different flavor transformations between neutrino production and detection\n", " * … an easy way to estimate event rates in neutrino detectors\n", "* SN models and flavor transformations can be used by other software\n", "* New SN models or flavor transformations are welcome!\n", "\n", "* **Interested? Have questions?**\n", " * Message me on Slack\n", " * E-mail me at jost.migenda@kcl.ac.uk\n", " * Or [open an issue on GitHub](https://github.com/SNEWS2/snewpy/issues)\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.9.5 ('snews')", "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.9.5" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "e2528887d751495e023d57d695389d9a04f4c4d2e5866aaf6dc03a1ed45c573e" } } }, "nbformat": 4, "nbformat_minor": 2 }