{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Neutrino Cross Sections\n", "\n", "Plot the cross sections for the key in-ice interactions for ~10 MeV neutrinos in IceCube." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from asteria import interactions\n", "from snewpy.neutrino import Flavor\n", "\n", "import astropy.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 importlib.resources import files\n", "\n", "from itertools import cycle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup styles for Plotting" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "axes_style = { 'grid' : 'True',\n", " 'labelsize' : '24',\n", " 'labelpad' : '8.0' }\n", "\n", "xtick_style = { 'direction' : 'out',\n", " 'labelsize' : '20.',\n", " 'major.size' : '5.', \n", " 'major.width' : '1.',\n", " 'minor.visible' : 'True',\n", " 'minor.size' : '2.5',\n", " 'minor.width' : '1.' }\n", "\n", "ytick_style = { 'direction' : 'out',\n", " 'labelsize' : '20.',\n", " 'major.size' : '5', \n", " 'major.width' : '1.',\n", " 'minor.visible' : 'True',\n", " 'minor.size' : '2.5',\n", " 'minor.width' : '1.' }\n", "\n", "grid_style = { 'alpha' : '0.75' }\n", "legend_style = { 'fontsize' : '16' }\n", "font_syle = { 'size' : '20'}\n", "text_style = { 'usetex' : 'True' }\n", "figure_style = { 'subplot.hspace' : '0.05' }\n", "\n", "mpl.rc( 'font', **font_syle )\n", "mpl.rc( 'text', **text_style )\n", "mpl.rc( 'axes', **axes_style )\n", "mpl.rc( 'xtick', **xtick_style )\n", "mpl.rc( 'ytick', **ytick_style )\n", "mpl.rc( 'grid', **grid_style )\n", "mpl.rc( 'legend', **legend_style )\n", "mpl.rc( 'figure', **figure_style )\n", "\n", "mpl.rcParams['text.usetex'] = True \n", "# mpl.rcParams['text.latex.preamble'] = [r'\\usepackage[cm]{sfmath}']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Neutrino Interactions\n", "\n", "There are several important neutrino interactions to consider:\n", "\n", "1. Inverse beta decay: $\\bar{\\nu}_e+p \\rightarrow e^{+} + n$.\n", "2. Elastic neutrino scattering: $\\nu_e + e^{-} \\rightarrow \\nu_e + e^{-}$ (plus antineutrino, plus $\\mu$ and $\\tau$).\n", "3. Oxygen-16 charged-current interaction: $\\nu_e + ^{16}\\mathrm{O}\\rightarrow e^{-} + \\mathrm{X}$ (plus antineutrino).\n", "4. Oxygen-16 neutral-current interaction: $\\nu_\\mathrm{all} + ^{16}\\mathrm{O} \\rightarrow \\nu_\\mathrm{all} + \\mathrm{X}$.\n", "5. Oxygen-18 interactions: $\\nu_e + ^{17/18}\\mathrm{O} / ^{2}_{1}\\mathrm{H} \\rightarrow e^{-} + \\mathrm{X}$.\n", "\n", "For details, see [R. Abbasi et al., A&A 535:A109, 2011](http://dx.doi.org/10.1051/0004-6361/201117810).\n", "\n", "__NOTE (3/13/19)__: Inverse Beta Decay has two implentations in ASTERIA, `InvBetaTab()` and `InvBetaPar()`. Below the latter is shown, but comparison plots are generated for both." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Interactions = [ interactions.InvBetaPar(),\n", " interactions.ElectronScatter(),\n", " interactions.Oxygen16CC(),\n", " interactions.Oxygen16NC(),\n", " interactions.Oxygen18() ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Neutrino Interactions\n", "\n", "The differential cross section as a function of neutrino energy is plotted versus neutrino energy for every flavor. Only non-zero cross sections are shown (IE Only the `Oxygen18()` cross section for $\\nu_e$ is shown as all other flavors' `Oxygen18()` cross sections are zero.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Enu = np.arange(0, 100, 0.1) * u.MeV\n", "lines = [\"-\", \"--\", \"-.\", \":\"]\n", "\n", "fig, ax = plt.subplots(1,1, figsize=(9,6))\n", "\n", "for interaction in Interactions:\n", " color = None\n", " line = cycle(lines)\n", " for flavor in Flavor:\n", " xs = interaction.cross_section(flavor, Enu).to(u.m**2)\n", " if xs.value.any():\n", " label='{}: {}'.format(interaction.__class__.__name__, flavor.to_tex())\n", " if color is None:\n", " p = ax.plot(Enu, xs, next(line), label=label)\n", " \n", " color = p[0].get_color()\n", " else:\n", " ax.plot(Enu, xs, next(line), label=label, color=color)\n", "\n", "\n", "\n", "ax.set(xlim=[0,100], ylim=[1e-49, 1e-42], yscale='log', title='Primary SN Interactions' )\n", "\n", "ax.set_xlabel(r'$E_\\nu$ [MeV]', horizontalalignment='right', x=1.0)\n", "ax.set_ylabel(r'$\\sigma(E_\\nu)$ [m$^2$]', horizontalalignment='right', y=1.0)\n", "\n", "ax.legend( bbox_to_anchor=(1.05,1))\n", "\n", "fig.subplots_adjust(left=0.075, right=0.8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define Helper Functions\n", "\n", "Define Functions for plotting, retrieving information from Data Files\n", "\n", "- `drawComparison` plots USSR's reported values against ASTERIA's reported values and computes the absolute difference and percent difference of ASTERIA's reported value relative to USSR's reported value. The average percent difference is denoted $\\bar{\\Delta}_\\%$. This average was taken over the entire curve excluding any points where USSR returned a value of 0.\n", "\n", " - `drawComparison` assumes that the two quantities being compared are plotted on the same domain, which is the case in this notebook.\n", " \n", " - In other notebooks, `drawComparison` plots the raw difference, here it plots the absolute difference for the cross sections, so it may be plotted with the cross section on a semi-log plot.\n", " \n", "\n", "- `getUSSRdata` retrieves the specified file from `data\\USSR\\`. The array size is specified in the first line of the file's header. \n", " - `getUSSRdata` assumes that the data is sorted into columns such that the first is time in seconds, the next four columns are USSR's reported cross sections for each flavor in the order $\\nu_e$, $\\bar{\\nu}_e$, $\\nu_x$, $\\bar{\\nu}_x$, and the next four are USSR's reported mean lepton energies (in the same order). The files in `data\\USSR\\` were written accordingly.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def drawComparison(t, ussr_y, astr_y, flavor, label='', units='', xs=False): \n", " # Compute difference and percent difference relative to USSR\n", " if xs:\n", " diff = abs(ussr_y - astr_y) \n", " else:\n", " diff = ussr_y - astr_y\n", " pct_diff = 100*np.divide( abs(diff) , ussr_y, \n", " where=ussr_y>0, \n", " out=np.zeros_like(ussr_y) )\n", " \n", " # Compute average percent difference excluding where USSR Reported 0\n", " avg_pct_diff = np.mean( pct_diff[pct_diff>0] )\n", " \n", " fig, (ax1, ax2) = plt.subplots(2,1, figsize=(9,6), \n", " gridspec_kw = {'height_ratios':[5, 1]}, \n", " sharex=True)\n", " \n", " # Plot USSR and ASTERIA against Each other.\n", " # Plot the abs difference between the two on the same figure.\n", " # Plot the percent difference on a subplot with a shared x-axis\n", " ax1.step( t, ussr_y, label='USSR')\n", " ax1.step( t, astr_y, label='ASTERIA')\n", " ax1.step( t, diff, 'k', label = 'Residual: USSR-ASTERIA', alpha=.50)\n", " if avg_pct_diff < 0.01:\n", " tmp_label = r'\\% Diff. $\\bar{\\Delta}_\\%$'+r'={0:6.1E}\\%'.format(avg_pct_diff)\n", " else:\n", " tmp_label = r'\\% Diff. $\\bar{\\Delta}_\\%$'+r'={0:4.2f}\\%'.format(avg_pct_diff)\n", " ax2.step( t, pct_diff, 'r', label = tmp_label)\n", " \n", " ttl = ax1.set_title( label+' Comparison: '+ flavor.to_tex() )\n", " # If this is a cross section plot, limit the y-axis and make it log-scaled.\n", " # - This was done to make the cross sections easier to read, even though\n", " # the graph of the residual is often cut off.\n", " if xs: \n", " ymin = np.min( np.append(ussr_y[ussr_y>0], astr_y[astr_y>0]) )\n", " ymax = np.max( np.append(ussr_y, astr_y) ) \n", " ax1.set( ylim=(ymin/2, ymax*2), yscale='log' )\n", " ttl.set_position([.5, 1.025])\n", " \n", " ax1.set_ylabel( label+' '+units, horizontalalignment='right', y = 1)\n", " ax2.set( xlim=(-1, 101), yscale='log')\n", " ax2.set_xlabel(r'$E_\\nu$ [MeV]', horizontalalignment='right', x = 1)\n", " \n", " \n", " # Group the labels of all lines in one legend object and display on axis 1\n", " handles, labels = ax1.get_legend_handles_labels()\n", " handles.append(ax2.get_legend_handles_labels()[0][0])\n", " labels.append(ax2.get_legend_handles_labels()[1][0])\n", " \n", " # Lock the legend to the lower right, In this notebook, this interferes with the \n", " # graphs the least, despite differing with loc='best' for the legend.\n", " legend = ax1.legend( handles, labels, loc='lower right' )\n", " \n", " ax1.set_title( label+' Comparison: '+ flavor.to_tex() )\n", " \n", " # Return figure and axis handles for additional manipulation.\n", " return fig, (ax1, ax2)\n", " \n", "def getUSSRdata(filename):\n", " fullpath = files('asteria.data').joinpath(f'USSR/{filename}')\n", " if not fullpath.exists():\n", " raise FileNotFoundError(f'{str(fullpath)} does not exist.')\n", "\n", " ussr = Table.read(fullpath,\n", " format='ascii',\n", " names=['Enu', 'xs_NU_E', 'xs_NU_E_BAR', 'xs_NU_X', 'xs_NU_X_BAR', \n", " 'Elep_NU_E', 'Elep_NU_E_BAR', 'Elep_NU_X', 'Elep_NU_X_BAR'])\n", "\n", " return ussr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Comparison of USSR's IBD with ASTERIA's Tabulated IBD\n", "Only $\\bar{\\nu}_e$ interact via IBD. This will compare the differential cross section and mean lepton energy reported by USSR with those reported by ASTERIA's `InvBetaTab()`. This implementation differs with USSR more than `InvBetaPar()` as a table is interpolated rather than using a parameterization. \n", "\n", "__NOTE (03/13/19)__: It might be worth checking `InvBetaTab()` and comparing with Vissani Table 1 to make sure everything is correct." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ussrdata = getUSSRdata('InvBeta.txt')\n", "Enu = ussrdata['Enu']\n", "interaction = interactions.InvBetaTab()\n", "\n", "for nu, flavor in enumerate(Flavor):\n", " xsname = f'xs_{flavor.name}'\n", " ussrxs = ussrdata[xsname]\n", " astrxs = interaction.cross_section(flavor=flavor, e_nu=Enu*u.MeV).to_value('m**2')\n", "\n", " if astrxs.any() and ussrxs.any():\n", " drawComparison(Enu, ussrxs, astrxs, flavor, r'$\\sigma(E_\\nu)$', '[m$^2$]', xs=True)\n", "\n", " mEname = f'Elep_{flavor.name}'\n", " ussrmE = ussrdata[mEname]\n", " astrmE = interaction.mean_lepton_energy(flavor=flavor, e_nu=Enu*u.MeV).to_value('MeV')\n", "\n", " if astrmE.any() and ussrmE.any():\n", " drawComparison(Enu, ussrmE, astrmE, flavor, r'$\\langle E \\rangle_\\mathrm{lepton}$', '[MeV]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Comparison of USSR's IBD with ASTERIA's Parameterized IBD\n", "Only $\\bar{\\nu}_e$ interact via IBD. This will compare the differential cross section and mean lepton energy reported by USSR with those reported by ASTERIA's `InvBetaPar()`. This implementation is the same as USSR's implementation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ussrdata = getUSSRdata('InvBeta.txt')\n", "Enu = ussrdata['Enu']\n", "interaction = interactions.InvBetaPar()\n", "\n", "for nu, flavor in enumerate(Flavor):\n", " xsname = f'xs_{flavor.name}'\n", " ussrxs = ussrdata[xsname]\n", " astrxs = interaction.cross_section(flavor=flavor, e_nu=Enu*u.MeV).to_value('m**2')\n", "\n", " if astrxs.any() and ussrxs.any():\n", " drawComparison(Enu, ussrxs, astrxs, flavor, r'$\\sigma(E_\\nu)$', '[m$^2$]', xs=True)\n", "\n", " mEname = f'Elep_{flavor.name}'\n", " ussrmE = ussrdata[mEname]\n", " astrmE = interaction.mean_lepton_energy(flavor=flavor, e_nu=Enu*u.MeV).to_value('MeV')\n", "\n", " if astrmE.any() and ussrmE.any():\n", " drawComparison(Enu, ussrmE, astrmE, flavor, r'$\\langle E \\rangle_\\mathrm{lepton}$', '[MeV]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Electron Scattering Comparisons for Each Flavor\n", "All Flavors interact via Electron Scattering. This will compare each flavor's differential cross section and mean energy as it is reported by ASTERIA and USSR. In this case, the mean lepton energy that is reported is the product of the differential cross section with the lepton energy, integrated w.r.t lepton energy, It has units m$^2$ MeV.\n", "\n", "__NOTE (03/13/19)__: I have thought about changing the implementation of `ElectronScatter()`'s mean energy to return a quantity with units MeV. That is, take the integrated product we have now, and divide it by the differential cross section. This is NOT how USSR implements it, and at the time of writing ASTERIA does NOT do this.\n", "\n", "__NOTE (04/22/19)__: I have implemented the change made in my previous comment. Taking a slight hit to performance, the `ElectronScatter.mean_lepton_energy` method now returns units $MeV$. This is a achieved by dividing the previous result by the cross section." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "astrmE = interactions.ElectronScatter().mean_lepton_energy(flavor=Flavor.NU_E, e_nu=Enu*u.MeV).to_value('MeV')\n", "ussrdata" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ussrdata = getUSSRdata('ElectronScatter.txt')\n", "Enu = ussrdata['Enu']\n", "interaction = interactions.ElectronScatter()\n", "\n", "for nu, flavor in enumerate(Flavor):\n", " xsname = f'xs_{flavor.name}'\n", " ussrxs = ussrdata[xsname]\n", " astrxs = interaction.cross_section(flavor=flavor, e_nu=Enu*u.MeV).to_value('m**2')\n", "\n", " if astrxs.any() and ussrxs.any():\n", " drawComparison(Enu, ussrxs, astrxs, flavor, r'$\\sigma(E_\\nu)$', '[m$^2$]', xs=True)\n", "\n", " mEname = f'Elep_{flavor.name}'\n", " ussrmE = ussrdata[mEname]\n", " astrmE = interaction.mean_lepton_energy(flavor=flavor, e_nu=Enu*u.MeV).to_value('MeV') * astrxs\n", "\n", " if astrmE.any() and ussrmE.any():\n", " drawComparison(Enu, ussrmE, astrmE, flavor, r'$\\sigma \\times \\langle E \\rangle_\\mathrm{lepton}$', '[m$^2$ MeV]', xs=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Oxygen-16 Charged Current Comparisons for $\\nu_e$ and $\\bar{\\nu}_e$\n", "Only $\\nu_e$ and $\\bar{\\nu}_e$ interact via charged currents with Oxygen-16. This will compare each flavor's differential cross section and mean energy as it is reported by ASTERIA and USSR. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ussrdata = getUSSRdata('Oxygen16CC.txt')\n", "Enu = ussrdata['Enu']\n", "interaction = interactions.Oxygen16CC()\n", "\n", "for nu, flavor in enumerate(Flavor):\n", " xsname = f'xs_{flavor.name}'\n", " ussrxs = ussrdata[xsname]\n", " astrxs = interaction.cross_section(flavor=flavor, e_nu=Enu*u.MeV).to_value('m**2')\n", "\n", " if astrxs.any() and ussrxs.any():\n", " drawComparison(Enu, ussrxs, astrxs, flavor, r'$\\sigma(E_\\nu)$', '[m$^2$]', xs=True)\n", "\n", " mEname = f'Elep_{flavor.name}'\n", " ussrmE = ussrdata[mEname]\n", " astrmE = interaction.mean_lepton_energy(flavor=flavor, e_nu=Enu*u.MeV).to_value('MeV')\n", "\n", " if astrmE.any() and ussrmE.any():\n", " drawComparison(Enu, ussrmE, astrmE, flavor, r'$\\langle E \\rangle_\\mathrm{lepton}$', '[MeV]', xs=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Oxygen-16 Neutral Current Comparisons for Each Flavor\n", "All Flavors interact via neutral currents with Oxygen-16. This will compare each flavor's differential cross section and mean energy as it is reported by ASTERIA and USSR. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ussrdata = getUSSRdata('Oxygen16NC.txt')\n", "Enu = ussrdata['Enu']\n", "interaction = interactions.Oxygen16NC()\n", "\n", "for nu, flavor in enumerate(Flavor):\n", " xsname = f'xs_{flavor.name}'\n", " ussrxs = ussrdata[xsname]\n", " astrxs = interaction.cross_section(flavor=flavor, e_nu=Enu*u.MeV).to_value('m**2')\n", "\n", " if astrxs.any() and ussrxs.any():\n", " drawComparison(Enu, ussrxs, astrxs, flavor, r'$\\sigma(E_\\nu)$', '[m$^2$]', xs=True)\n", "\n", " mEname = f'Elep_{flavor.name}'\n", " ussrmE = ussrdata[mEname]\n", " astrmE = interaction.mean_lepton_energy(flavor=flavor, e_nu=Enu*u.MeV).to_value('MeV')\n", "\n", " if astrmE.any() and ussrmE.any():\n", " drawComparison(Enu, ussrmE, astrmE, flavor, r'$\\langle E \\rangle_\\mathrm{lepton}$', '[MeV]', xs=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Oxygen-18 Charged Current Comparisons for $\\nu_e$.\n", "Only $\\nu_e$ interact via charged currents with Oxygen-18. This will compare each flavor's differential cross section and mean energy as it is reported by ASTERIA and USSR. \n", "\n", "__NOTE (03/13/19)__: This parameterization was obtained using a quadratic fit to cross section estimated from Kamiokande data from Haxton and Robertson, PRC 59:515, 1999. *See also*, the page on the Mainz Wiki, Neutrino cross sections on natural oxygen. The cross section includes one additional scaling factor of one over the abundance of O$^{18}$, and it is unclear why this is present. It would make more sense if the cross sectin was scaled by the oxygen abundance. This scaling factor currently features both in USSR and ASTERIA." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ussrdata = getUSSRdata('Oxygen18.txt')\n", "Enu = ussrdata['Enu']\n", "interaction = interactions.Oxygen18()\n", "\n", "for nu, flavor in enumerate(Flavor):\n", " xsname = f'xs_{flavor.name}'\n", " ussrxs = ussrdata[xsname]\n", " astrxs = interaction.cross_section(flavor=flavor, e_nu=Enu*u.MeV).to_value('m**2')\n", "\n", " if astrxs.any() and ussrxs.any():\n", " drawComparison(Enu, ussrxs, astrxs, flavor, r'$\\sigma(E_\\nu)$', '[m$^2$]', xs=True)\n", "\n", " mEname = f'Elep_{flavor.name}'\n", " ussrmE = ussrdata[mEname]\n", " astrmE = interaction.mean_lepton_energy(flavor=flavor, e_nu=Enu*u.MeV).to_value('MeV')\n", "\n", " if astrmE.any() and ussrmE.any():\n", " drawComparison(Enu, ussrmE, astrmE, flavor, r'$\\langle E \\rangle_\\mathrm{lepton}$', '[MeV]', xs=True)" ] } ], "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": 4 }