{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Detector Hits + Oscillations 2 (Obsolete)\n", "\n", "Demonstrate the total number of hits in the IceCube detector after implementaing neutrino oscillations for a progentitor at distance 1 kpc.\n", "The following data is saved to a FITS file:\n", "1. Luminosity corresponding to each neutrino flavor before mixing\n", "2. Neutrino flux for each flavor in each oscillation scenario (unmixed, normal mixing, and inverted mixing)\n", "3. Signal generated for each flavor in each oscillation scenario\n", "\n", "\n", "The progenitor model used is the Nakazato Shen model with z = 0.02, t$_{rev}$ = 300ms, m = 13.0M$_\\odot$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from astropy import units as u\n", "from astropy.table import Table\n", "\n", "import numpy as np\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "\n", "from asteria import config, source, detector\n", "from asteria.interactions import Interactions\n", "from asteria.neutrino import Flavor\n", "from asteria.oscillation import SimpleMixing\n", "import asteria.IO as io\n", "\n", "mpl.rc('font', size=16)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load Configuration\n", "\n", "This will load the source configuration from a file.\n", "\n", "For this to work, either the user needs to have done one of two things:\n", "1. Run `python setup.py install` in the ASTERIA directory.\n", "2. Run `python setup.py develop` and set the environment variable `ASTERIA` to point to the git source checkout.\n", "\n", "If these were not done, the initialization will fail because the paths will not be correctly resolved." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "conf = config.load_config('../../data/config/nakazato-shen-z0.02-t_rev300ms-s13.0.yaml')\n", "ccsn = source.initialize(conf)\n", "ic86 = detector.initialize(conf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Prepare Iterables\n", "\n", "Define the range of neutrino energies (`E_nu`) to simulate and the times (`time`) and distances (`dist`) at which to perform the simulation. \n", "\n", "The progenitor distance is set to 1 kpc." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "E_min = 0.1; E_max = 100.1; dE = 0.1;\n", "Enu = np.arange(E_min, E_max, dE) * u.MeV\n", "\n", "t_min = -1; t_max = 15; dt = 0.0001;\n", "time = np.arange(t_min, t_max, dt) * u.s\n", "\n", "ccsn.progenitor_distance = 1 * u.kpc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Redefine `photonic_energy_per_volume`\n", "\n", "`photonic_energy_per_volume` is originally defined in the `source.py` module. The function signature has been modified for the purposes of this program to read in the neutrino flux for each flavor as an ndarray to account for flavor mixing, as opposed to calculcing it using `source.get_flux(time, flavor)`." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def photonic_energy_per_vol(source, time, E, flavor, photon_spectrum, flux_list, n=1000):\n", " \"\"\"Compute the energy deposited in a cubic meter of ice by photons\n", " from SN neutrino interactions.\n", "\n", " Parameters\n", " ----------\n", "\n", " time : float (units s)\n", " Time relative to core bounce. \n", " E : `numpy.ndarray`\n", " Sorted grid of neutrino energies\n", " flavor : :class:`asteria.neutrino.Flavor`\n", " Neutrino flavor.\n", " photon_spectrum : `numpy.ndarray` (Units vary, m**2)\n", " Grid of the product of lepton cross section with lepton mean energy\n", " and lepton path length per MeV, sorted according to parameter E\n", " flux_list: `numpy.ndarray`\n", " List containing neutrino fluxes corresponding to nu_e, nu_e_bar, nu_x, nu_x_bar\n", " n : int\n", " Maximum number of time steps to compute at once. A temporary numpy array\n", " of size n x time.size is created and can be very memory inefficient.\n", "\n", " Returns\n", " -------\n", " E_per_V\n", " Energy per m**3 of ice deposited by neutrinos of requested flavor\n", " \"\"\"\n", " H2O_in_ice = 3.053e28 # 1 / u.m**3\n", " \n", " t = time.to(u.s).value\n", " Enu = E.to(u.MeV)\n", " if Enu[0] == 0:\n", " Enu[0] = 1e-10 * u.MeV\n", " phot = photon_spectrum.to(u.m**2).value.reshape((-1,1)) # m**2\n", " \n", " dist = source.progenitor_distance.to(u.m).value # m**2\n", "# flux = source.get_flux( time, flavor ) # Unitless\n", " flux = flux_list\n", " if not flavor.is_electron:\n", " flux *= 2\n", " \n", " print('Beginning {0} simulation....'.format(flavor._name_), end='')\n", " # The following two lines exploit the fact that astropy quantities will\n", " # always return a number when numpy size is called on them, even if it is 1.\n", " E_per_V = np.zeros( time.size ) \n", " for i_part in source.parts_by_index(time, n): # Limits memory usage\n", " E_per_V[i_part] += np.trapz( source.energy_spectrum(t[i_part], Enu, flavor) * phot, Enu.value, axis=0)\n", " E_per_V *= flux * H2O_in_ice / ( 4 * np.pi * dist**2)\n", " print('Completed')\n", " \n", " return E_per_V * u.MeV / u.m**3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute Charged Particle Spectrum\n", "\n", "Compute the number of photons produced by $\\nu$ component interactions with charged particles given neutrino flavor and energy. `Interactions` contains a list of the interactions that are simulated. This list may be changed to turn 'off/on' specific interactions\n", "\n", "The interactions are as follows:\n", "- `InvBetaTab` : Tabulated inverse beta decay computation by Strumia and Vissani, Phys. Lett. B 564:42, 2003.\n", " - See Also: `InvBetaPar` : Inverse beta decay parameterization\n", "- `ElectronScatter` : Elastic Neutrino-electron scattering from Marciano and Parsa, J. Phys. G 29:2969, 2003.\n", "- `Oxygen16CC` : $\\nu$-$^{16}O$ charged current interaction, using estimates from Kolbe et al., PRD 66:013007, 2002.\n", "- `Oxygen16NC` : $\\nu$-$^{16}O$ neutral current interaction, using estimates from Kolbe et al., PRD 66:013007, 2002.\n", "- `Oxygen18` : $\\nu$-$^{18}O$ charged current interaction, using estimates from Kamiokande data from Haxton and Robertson, PRC 59:515, 1999.\n", "\n", "These `Interaction` objects may be used to compute the neutrino cross sections and mean energy of the produced lepton, both as a function of neutrino energy. The final state lepton energy has been integrated out. This cross section with a component of $H_2O$ is then scaled as appropriate for a $H_2O$ molecule (IE Electron scattering cross section is scaled by 10, as there are 10 electrons in $H_2O$).\n", "\n", "`photon_scaling_factor` is the number of photones per MeV of lepton energy. It is computed by taking product of the data members `photons_per_lepton_MeV` and `p2e_path_ratio` which are respectively, the number of photons emitted per unit lepton path length, and the ratio of positron path length to electron path length in ice.\n", "\n", "`photons_per_lepton_MeV` is computed by finding number of photon emitted per unit lepton path length and multiplying it by the lepton path length per MeV. This is done using the Frank-Tamm formula and index of refraction from Price and Bergstrom, AO 36:004181, 1997.\n", "\n", "This result estimates the number of photons as a function of neutrino energy. It will have units $m^2$ at the end of this cell but is later scaled by the $r^2$ where $r$ is the progenitor distance, accounting for losses." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "photon_spectra = np.zeros( shape=(len(Flavor), Enu.size) )\n", "\n", "for nu, flavor in enumerate(Flavor):\n", " for interaction in Interactions: \n", " xs = interaction.cross_section(flavor, Enu).to(u.m**2).value\n", " E_lep = interaction.mean_lepton_energy(flavor, Enu).value\n", " photon_scaling_factor = interaction.photon_scaling_factor(flavor).to( 1/u.MeV).value\n", " photon_spectra[nu] += xs * E_lep * photon_scaling_factor # u.m**2 \n", "photon_spectra *= u.m**2 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define Helper Function" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def generate_signal(nu_list):\n", " '''Compute the number of hits generated in the detector for a given list of fluxes.\n", " The signal generated is by default binned to 1e-4s, i.e. the `time` binning.\n", " \n", " Parameter\n", " ---------\n", " nu_list: ndarray\n", " List of neutrino flavors, ordered by default as nu_e, nu_e_bar, nu_x, nu_x_bar\n", " timebin: float\n", " The desired time binning used to scale the deadtime fractions\n", " \n", " Returns\n", " -------\n", " hits: ndarray\n", " Number of hits in the detector\n", " '''\n", " signal = np.zeros(shape = (len(Flavor), time.size))\n", " \n", " for nu, (flavor, photon_spectrum) in enumerate(zip(Flavor, photon_spectra)):\n", " signal[nu] = photonic_energy_per_vol(ccsn, time, Enu, flavor, photon_spectrum, nu_list[nu])\n", "\n", " return signal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementing Neutrino Oscillations\n", "\n", "Fluxes for each flavor are calculted for normal and inverted mixing using the class `SimpleMixing` in the `neutrino.py` module.\n", "\n", "$\\theta_{12} = 33.82^{\\circ}$ is used to initialize an object `mix` of `SimpleMixing`." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "nu_list = np.zeros(shape = (4, 160000))\n", "mix = SimpleMixing(33.82)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "for i, flavor in enumerate(Flavor):\n", " nu_list[i] = ccsn.get_flux(time, flavor)\n", " \n", "nu_norm = mix.normal_mixing(nu_list)\n", "nu_inv = mix.inverted_mixing(nu_list)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "[u_e_f, u_e_bar_f, u_x_f, u_x_bar_f] = nu_list\n", "[n_e_f, n_e_bar_f, n_x_f, n_x_bar_f] = nu_norm\n", "[i_e_f, i_e_bar_f, i_x_f, i_x_bar_f] = nu_inv" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Beginning nu_e simulation....Completed\n", "Beginning nu_e_bar simulation....Completed\n", "Beginning nu_x simulation....Completed\n", "Beginning nu_x_bar simulation....Completed\n" ] } ], "source": [ "sig_unmixed = generate_signal(nu_list)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "[u_e_s, u_e_bar_s, u_x_s, u_x_bar_s] = sig_unmixed" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Beginning nu_e simulation....Completed\n", "Beginning nu_e_bar simulation....Completed\n", "Beginning nu_x simulation....Completed\n", "Beginning nu_x_bar simulation....Completed\n" ] } ], "source": [ "sig_norm = generate_signal(nu_norm)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "[n_e_s, n_e_bar_s, n_x_s, n_x_bar_s] = sig_norm" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Beginning nu_e simulation....Completed\n", "Beginning nu_e_bar simulation....Completed\n", "Beginning nu_x simulation....Completed\n", "Beginning nu_x_bar simulation....Completed\n" ] } ], "source": [ "sig_inv = generate_signal(nu_inv)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "[i_e_s, i_e_bar_s, i_x_s, i_x_bar_s] = sig_inv" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "lum = np.zeros(shape = (len(Flavor), len(time)))\n", "for nu, flavor in enumerate(Flavor):\n", " lum[nu] = ccsn.get_luminosity(time, flavor)\n", "[e_lum, e_bar_lum, x_lum, x_bar_lum] = lum" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# Create astropy table containing all relevant data\n", "\n", "t = Table([time, \n", " e_lum, e_bar_lum, x_lum, x_bar_lum,\n", " u_e_f, u_e_bar_f, u_x_f, u_x_bar_f,\n", " n_e_f, n_e_bar_f, n_x_f, n_x_bar_f,\n", " i_e_f, i_e_bar_f, i_x_f, i_x_bar_f,\n", " u_e_s, u_e_bar_s, u_x_s, u_x_bar_s,\n", " n_e_s, n_e_bar_s, n_x_s, n_x_bar_s,\n", " i_e_s, i_e_bar_s, i_x_s, i_x_bar_s],\n", " names = ('time',\n", " 'nu_e luminosity', 'nu_e_bar luminosity','nu_x luminosity','nu_x_bar luminosity',\n", " 'unmixed nu_e flux', 'unmixed nu_e_bar flux', 'unmixed nu_x flux', 'unmixed nu_x_bar flux',\n", " 'normal nu_e flux', 'normal nu_e_bar flux', 'normal nu_x flux', 'normal nu_x_bar flux',\n", " 'inverted nu_e flux', 'inverted nu_e_bar flux', 'inverted nu_x flux', 'inverted nu_x_bar flux',\n", " 'unmixed nu_e sig', 'unmixed nu_e_bar sig', 'unmixed nu_x sig', 'unmixed nu_x_bar sig',\n", " 'normal nu_e sig', 'normal nu_e_bar sig', 'normal nu_x sig', 'normal nu_x_bar sig',\n", " 'inverted nu_e sig', 'inverted nu_e_bar sig', 'inverted nu_x sig', 'inverted nu_x_bar sig'\n", " ))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "t.write('nakazato-shen-z0.02-t_rev300ms-s13.0_e_per_v.fits', format = 'fits', overwrite = True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }