Interaction Contributions
Loop through the set of standard IceCube interactions and plot the contributions to the signal from each under various flavor transformation scenarios and hierarchies.
[1]:
from asteria.simulation import Simulation
from asteria import set_rcparams
from asteria import interactions
from snewpy.neutrino import Flavor
from scipy.interpolate import PchipInterpolator
import astropy.units as u
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import os
from tqdm import tqdm
set_rcparams(verbose=False)
/home/docs/checkouts/readthedocs.org/user_builds/asteria/envs/latest/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
In-Ice Interactions
See details in R. Abbasi+, A&A 535:A109, 2011. We’ll plot the interactions in reverse order of their relative importance.
\(^{18}\mathrm{O}\) CC interaction.
\(^{16}\mathrm{O}\) NC interaction.
\(^{16}\mathrm{O}\) CC interaction.
Electron elastic scattering.
Inverse beta decay.
[2]:
model = {'name': 'Nakazato_2013',
'param':{
'progenitor_mass': 13 * u.Msun,
'revival_time': 300 * u.ms,
'metallicity': 0.02,
'eos': 'shen'}
}
sims = []
xform = 'NoTransformation'
nmo = 'Normal'
#- Loop through individual interactions and compute the signal from each.
processes = (interactions.Oxygen18,
interactions.Oxygen16NC,
interactions.Oxygen16CC,
interactions.ElectronScatter,
interactions.InvBetaPar)
procnames = ('$^{18}\mathrm{O}$ CC',
'$^{16}\mathrm{O}$ NC',
'$^{16}\mathrm{O}$ CC',
'Electron elastic scattering',
'Inverse beta decay')
for proc in processes:
sim = Simulation(model=model,
interactions=[proc],
distance=10 * u.kpc,
Emin=0*u.MeV, Emax=100*u.MeV, dE=1*u.MeV,
tmin=-10*u.s, tmax=20*u.s, dt=1*u.ms,
mixing_scheme=xform, hierarchy=nmo)
sim.run()
sims.append(sim)
<>:19: SyntaxWarning: invalid escape sequence '\m'
<>:20: SyntaxWarning: invalid escape sequence '\m'
<>:21: SyntaxWarning: invalid escape sequence '\m'
<>:19: SyntaxWarning: invalid escape sequence '\m'
<>:20: SyntaxWarning: invalid escape sequence '\m'
<>:21: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_1538/901487725.py:19: SyntaxWarning: invalid escape sequence '\m'
procnames = ('$^{18}\mathrm{O}$ CC',
/tmp/ipykernel_1538/901487725.py:20: SyntaxWarning: invalid escape sequence '\m'
'$^{16}\mathrm{O}$ NC',
/tmp/ipykernel_1538/901487725.py:21: SyntaxWarning: invalid escape sequence '\m'
'$^{16}\mathrm{O}$ CC',
Plot the Contributions of Each Interaction
Make a stacked plot fo the total signal as a function of time.
[3]:
binsize = 500*u.ms
bg = None
signal_hits = []
for sim in sims:
print(sim.metadata['interactions'])
t, dmu = sim.detector_hits(binsize)
signal_hits.append(dmu)
if bg is None:
bg = sim.detector.i3_bg(binsize, size=dmu.size) + sim.detector.dc_bg(binsize, size=dmu.size)
hits = (dmu + bg)
print(np.sum(dmu), np.sum(bg))
Oxygen18
499.95154123590066 45638120.232824564
Oxygen16NC
3209.1714647838785 45638120.232824564
Oxygen16CC
15064.733045986883 45638120.232824564
ElectronScatter
9474.336299579074 45638120.232824564
InvBetaPar
354723.3506779512 45638120.232824564
[4]:
fig, ax = plt.subplots(figsize=(8,5), tight_layout=True)
ax.stackplot(t, np.vstack(signal_hits), step='post', labels=procnames)
ax.set(xlim=(0,20),
xlabel='time [s]',
ylim=(1,1e6),
yscale='log',
ylabel='signal hits')
ax.legend(fontsize=12);
Divide Hits into Phases of the CCSN
Plot hits showing the distinct phases of the explosion.
Consider also different oscillation scenarios and neutrino mass hierarchies.
[5]:
# Time limits for phases:
limits = [
(-0.025, 0.100), # Deleptonization, <0.1 s
(0.1, 0.6), # Accretion, 0.1 - 0.6 s
(1, 10), # Cooling, 0.6 - 10 s
]
scale = 1e4
binnings = [4e-3, 10e-3, 0.25] * u.s
No Flavor Transformations
[6]:
fig, axes = plt.subplots(1,3, figsize = (17,5))
bbox_style = {'boxstyle': 'round', 'edgecolor': 'gray', 'facecolor': 'white', 'alpha': 0.75}
for ax, binsize, xlim in zip(axes, binnings, limits):
bg = None
signal_hits = []
# Simulations iterate by mixing scheme (see cell 2)
for sim in sims:
label='x'
color='k'
# for sim, label, color in zip(sims, osc_labels, colors):
# Generate Signal hits
t, dmu = sim.detector_hits(binsize)
signal_hits.append(dmu / scale)
# ax.step(t, hits, label=label, c=color)
# Add 1-sigma band around the expected hits, assuming Poisson uncertainties.
# hits_up = ((dmu + bg) + np.sqrt(dmu + bg))/scale
# hits_lo = ((dmu + bg) - np.sqrt(dmu + bg))/scale
# ax.fill_between(t, hits_lo, hits_up, step='pre', color=color, alpha=0.25)
ax.stackplot(t, np.vstack(signal_hits), step='post', labels=procnames)
ax.set(ylim=(0, None))
# Normalized to single dom rate in Hz
# ax.step(t, bg/5160/binsize.to(u.s).value, label='Background', c='k', alpha=0.75)
ax.set(xlim=xlim)
if binsize <= 100 * u.ms:
scaled_binsize = binsize.to(u.ms)
annotation = f'{scaled_binsize.value} {scaled_binsize.unit} bins'
else:
annotation = f'{binsize.value} {binsize.unit} bins'
ax.text(0.05, 0.925, annotation, bbox=bbox_style, horizontalalignment='left',
verticalalignment='center', transform = ax.transAxes, fontsize=16)
# Plot background
axes[0].set_title('Deleptonization')
axes[0].set_ylabel(r'signal hits [$\times10^4$]', horizontalalignment='right', y=1.0)
axes[0].set_ylim(0, 0.9)
axes[1].set_title('Accretion')
axes[1].set_ylim(0, 1.5)
axes[2].set_title('PNS Cooling')
axes[2].set_xlabel(r'$t-t_\mathrm{bounce}$ [s]')
axes[2].set_ylim(0, 10)
axes[2].legend(loc='upper right', ncol=1, fontsize = 12);
Adiabatic MSW Effect + Normal Ordering
Run a new simulation and regenerate the plots.
[7]:
xform = 'AdiabaticMSW'
nmo = 'Normal'
sims = []
for proc in processes:
sim = Simulation(model=model,
interactions=[proc],
distance=10 * u.kpc,
Emin=0*u.MeV, Emax=100*u.MeV, dE=1*u.MeV,
tmin=-10*u.s, tmax=20*u.s, dt=1*u.ms,
mixing_scheme=xform, hierarchy=nmo)
sim.run()
sims.append(sim)
[8]:
fig, axes = plt.subplots(1,3, figsize = (17,5))
bbox_style = {'boxstyle': 'round', 'edgecolor': 'gray', 'facecolor': 'white', 'alpha': 0.75}
for ax, binsize, xlim in zip(axes, binnings, limits):
bg = None
signal_hits = []
# Simulations iterate by mixing scheme (see cell 2)
for sim in sims:
label='x'
color='k'
# for sim, label, color in zip(sims, osc_labels, colors):
# Generate Signal hits
t, dmu = sim.detector_hits(binsize)
signal_hits.append(dmu / scale)
# ax.step(t, hits, label=label, c=color)
# Add 1-sigma band around the expected hits, assuming Poisson uncertainties.
# hits_up = ((dmu + bg) + np.sqrt(dmu + bg))/scale
# hits_lo = ((dmu + bg) - np.sqrt(dmu + bg))/scale
# ax.fill_between(t, hits_lo, hits_up, step='pre', color=color, alpha=0.25)
ax.stackplot(t, np.vstack(signal_hits), step='post', labels=procnames)
ax.set(ylim=(0, None))
# Normalized to single dom rate in Hz
# ax.step(t, bg/5160/binsize.to(u.s).value, label='Background', c='k', alpha=0.75)
ax.set(xlim=xlim)
if binsize <= 100 * u.ms:
scaled_binsize = binsize.to(u.ms)
annotation = f'{scaled_binsize.value} {scaled_binsize.unit} bins'
else:
annotation = f'{binsize.value} {binsize.unit} bins'
ax.text(0.05, 0.925, annotation, bbox=bbox_style, horizontalalignment='left',
verticalalignment='center', transform = ax.transAxes, fontsize=16)
# Plot background
axes[0].set_title('Deleptonization')
axes[0].set_ylabel(r'signal hits [$\times10^4$]', horizontalalignment='right', y=1.0)
axes[0].set_ylim(0, 0.9)
axes[1].set_title('Accretion')
axes[1].set_ylim(0, 1.5)
axes[2].set_title('PNS Cooling')
axes[2].set_xlabel(r'$t-t_\mathrm{bounce}$ [s]')
axes[2].set_ylim(0, 10)
axes[2].legend(loc='upper right', ncol=1, fontsize = 12);
Adiabatic MSW Effect + Inverted Ordering
Run a new simulation and regenerate the plots.
[9]:
xform = 'AdiabaticMSW'
nmo = 'Inverted'
sims = []
for proc in processes:
sim = Simulation(model=model,
interactions=[proc],
distance=10 * u.kpc,
Emin=0*u.MeV, Emax=100*u.MeV, dE=1*u.MeV,
tmin=-10*u.s, tmax=20*u.s, dt=1*u.ms,
mixing_scheme=xform, hierarchy=nmo)
sim.run()
sims.append(sim)
[10]:
fig, axes = plt.subplots(1,3, figsize = (17,5))
bbox_style = {'boxstyle': 'round', 'edgecolor': 'gray', 'facecolor': 'white', 'alpha': 0.75}
for ax, binsize, xlim in zip(axes, binnings, limits):
bg = None
signal_hits = []
# Simulations iterate by mixing scheme (see cell 2)
for sim in sims:
label='x'
color='k'
# for sim, label, color in zip(sims, osc_labels, colors):
# Generate Signal hits
t, dmu = sim.detector_hits(binsize)
signal_hits.append(dmu / scale)
# ax.step(t, hits, label=label, c=color)
# Add 1-sigma band around the expected hits, assuming Poisson uncertainties.
# hits_up = ((dmu + bg) + np.sqrt(dmu + bg))/scale
# hits_lo = ((dmu + bg) - np.sqrt(dmu + bg))/scale
# ax.fill_between(t, hits_lo, hits_up, step='pre', color=color, alpha=0.25)
ax.stackplot(t, np.vstack(signal_hits), step='post', labels=procnames)
ax.set(ylim=(0, None))
# Normalized to single dom rate in Hz
# ax.step(t, bg/5160/binsize.to(u.s).value, label='Background', c='k', alpha=0.75)
ax.set(xlim=xlim)
if binsize <= 100 * u.ms:
scaled_binsize = binsize.to(u.ms)
annotation = f'{scaled_binsize.value} {scaled_binsize.unit} bins'
else:
annotation = f'{binsize.value} {binsize.unit} bins'
ax.text(0.05, 0.925, annotation, bbox=bbox_style, horizontalalignment='left',
verticalalignment='center', transform = ax.transAxes, fontsize=16)
# Plot background
axes[0].set_title('Deleptonization')
axes[0].set_ylabel(r'signal hits [$\times10^4$]', horizontalalignment='right', y=1.0)
axes[0].set_ylim(0, 0.9)
axes[1].set_title('Accretion')
axes[1].set_ylim(0, 1.5)
axes[2].set_title('PNS Cooling')
axes[2].set_xlabel(r'$t-t_\mathrm{bounce}$ [s]')
axes[2].set_ylim(0, 10)
axes[2].legend(loc='upper right', ncol=1, fontsize = 12);
