SNEWPY Demo for PyHEP 2022¶

import snewpy

snewpy.get_models(models=["Nakazato_2013", "Bollig_2016"])

# You can also run `snewpy.get_models()` without arguments,
# to interactively select which model files to download
#snewpy.get_models()

from snewpy.models.ccsn import Nakazato_2013, Bollig_2016

nakazato = Nakazato_2013('SNEWPY_models/Nakazato_2013/nakazato-shen-z0.004-t_rev100ms-s20.0.fits')
bollig = Bollig_2016('SNEWPY_models/Bollig_2016/s27.0c')  # This model has one file per flavor. Use common prefix, not full filename.

nakazato

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt

from snewpy.neutrino import Flavor

mpl.rc('font', size=16)

# Plot luminosity of both models
fig, ax = plt.subplots(1, figsize=(10, 6))

for flavor in Flavor:
    ax.plot(nakazato.time, nakazato.luminosity[flavor]/1e51,  # Report luminosity in units foe/s
            label=flavor.to_tex() + ' (Nakazato)',
            color='C0' if flavor.is_electron else 'C2',
            ls='-' if flavor.is_neutrino else '--',
            lw=2)

# for flavor in Flavor:
#     ax.plot(bollig.time, bollig.luminosity[flavor]/1e51,  # Report luminosity in units foe/s
#             label=flavor.to_tex() + ' (Bollig)',
#             color='C1' if flavor.is_electron else 'C3',
#             ls='-' if flavor.is_neutrino else '--',
#             lw=1)

ax.set(xlim=(-0.05, 0.5), xlabel=r'$t-t_{\rm bounce}$ [s]', ylabel=r'luminosity [$10^{51}$ erg s$^{-1}$]')
ax.grid()
ax.legend(loc='upper right', ncol=2, fontsize=18);

from snewpy.flavor_transformation import AdiabaticMSW
from snewpy.neutrino import Flavor, MassHierarchy

from astropy import units as u
import matplotlib as mpl
import numpy as np

def plot_total_flux(model, xform_nmo, xform_imo):
    """Plot initial and oscillated neutrino luminosities."""
    
    energies = np.linspace(0,60,121) * u.MeV    
    d = (10*u.kpc).to('cm').value # distance to SN
        
    times = model.get_time()
    burst_epoch = times <= 0.1*u.s
    accretion_epoch = (times > 0.1*u.s) & (times <= 0.5*u.s)
    cooling_epoch = (times > 0.5*u.s) & (times <= 10*u.s)
    
    ilum = {}
    tlum_nmo = {}
    tlum_imo = {}
    
    for flavor in Flavor:
        ilum[flavor] = np.zeros(len(times))
        tlum_nmo[flavor] = np.zeros(len(times))
        tlum_imo[flavor] = np.zeros(len(times))

    # Compute the transformed and untransformed flux at each time.
    for i, t in enumerate(times):
        ispec = model.get_initial_spectra(t, energies)
        tspec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)
        tspec_imo = model.get_transformed_spectra(t, energies, xform_imo)

        for flavor in Flavor:
            for j, E in enumerate(energies):
                ispec[flavor][j] /= (4.*np.pi*d**2)
                tspec_nmo[flavor][j] /= (4.*np.pi*d**2)
                tspec_imo[flavor][j] /= (4.*np.pi*d**2)
        
        for flavor in Flavor:
            ilum[flavor][i] = np.trapz(ispec[flavor].to('1/(erg*s)'), energies.to('erg')).value 
            tlum_nmo[flavor][i] = np.trapz(tspec_nmo[flavor].to('1/(erg*s)'), energies.to('erg')).value
            tlum_imo[flavor][i] = np.trapz(tspec_imo[flavor].to('1/(erg*s)'), energies.to('erg')).value            
            
    # make the figures 
    fig, axes = plt.subplots(3,3, figsize=(20,12), tight_layout=True)
    
    smax = [0.,0.,0.]
    titles = ['Untransformed', 'Transformed (NMO)', 'Transformed (IMO)']
    for i, spec in enumerate([ilum, tlum_nmo, tlum_imo]):
        for j, phase in enumerate([burst_epoch, accretion_epoch, cooling_epoch]):
            ax = axes[i,j]
            timeunits = 'ms' if j==0 else 's'
                
            for flavor in Flavor:
                if i == 0:
                    smax[j] = np.maximum(smax[j], 1.1*np.max(spec[flavor][phase]))
                    
                ax.plot(times[phase].to(timeunits),
                        spec[flavor][phase], label=flavor.to_tex(), lw=3,
                        color='C0' if flavor.is_electron else 'C1',
                        ls='-' if flavor.is_neutrino else ':')
            
            ax.set(xlim=(times[phase][0].to(timeunits).value, times[phase][-1].to(timeunits).value),
                   ylim=(0, smax[j]))
            
            if j==0:
               ax.set(ylabel=r'flux [cm$^{-2}$ s$^{-1}$]')
               ax.legend(loc='upper right', ncol=1, fontsize=18)
            if j==1:
                ax.set(title=titles[i])
            if i < 2:
                ax.set(xticklabels=[])
            else:
                ax.set(xlabel='time [{}]'.format(timeunits))
            
            ax.grid(ls=':')

    return fig

def plot_spectra(model, xform, t):
    """Plot initial and oscillated neutrino spectra at a fixed time."""

    energies = np.linspace(0,60,121) * u.MeV   
    d = (10*u.kpc).to('cm').value # distance to SN

    #get the spectra
    ispec = model.get_initial_spectra(t, energies) 
    tspec = model.get_transformed_spectra(t, energies, xform)

    for flavor in Flavor:
        for j, E in enumerate(energies):
            ispec[flavor][j] /= (4.*np.pi*d**2)
            tspec[flavor][j] /= (4.*np.pi*d**2)
    
    fig, axes = plt.subplots(1, 2, figsize=(18,7), sharey=True, tight_layout=True)

    for i, spec in enumerate([ispec, tspec]):
        axes[0].plot(energies, spec[Flavor.NU_E]/1e16, 
                    label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),
                    color='C0', ls='-' if i==0 else ':', lw=2,  alpha=0.7)
        axes[0].plot(energies, spec[Flavor.NU_X]/1e16, 
                    label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),
                    color='C1', ls='-' if i==0 else ':', lw=2,  alpha=0.7)

        axes[0].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Neutrino spectra at $t = ${:.1f}'.format(t))
        axes[0].grid()
        axes[0].legend(loc='upper right', ncol=2, fontsize=16)
    
        axes[1].plot(energies, spec[Flavor.NU_E_BAR]/1e16, 
                    label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),
                    color='C0', ls='-' if i==0 else ':', lw=2,  alpha=0.7)
        axes[1].plot(energies, spec[Flavor.NU_X_BAR]/1e16, 
                    label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),
                    color='C1', ls='-' if i==0 else ':', lw=2,  alpha=0.7)

        axes[1].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Antineutrino spectra at $t = ${:.1f}'.format(t)) 
        axes[1].grid()
        axes[1].legend(loc='upper right', ncol=2, fontsize=16)    
    
    ax = axes[0]
    ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')
    
    return fig

xform_nmo = AdiabaticMSW()
xform_imo = AdiabaticMSW(mh=MassHierarchy.INVERTED)

fig = plot_total_flux(nakazato, xform_nmo, xform_imo)
# fig.savefig('flux_adiabaticmsw.pdf')
# fig.show()

fig_nmo = plot_spectra(nakazato, xform_nmo, 100*u.ms)
# fig_nmo.show()
# fig_imo = plot_spectra(nakazato, xform_imo, 100*u.ms)
# fig_imo.show()

from snewpy import snowglobes

SNOwGLoBES_path = "./snowglobes/"  # where SNOwGLoBES is located
SNEWPY_models_base = "./SNEWPY_models/"  # directory containing SNEWPY model files

# set distance in kpc
distance = 10

# set SNOwGLoBES detector to use
detector = "wc100kt30prct"  # like Super- or Hyper-Kamiokande
# detector = "ar40kt"  # like DUNE

# set SNEWPY model type and filename
modeltype = 'Bollig_2016'
model = 's27.0c'

# set desired flavor transformation
transformation = 'NoTransformation' #'AdiabaticMSW_IMO'

# Construct file system path of model file and name of output file
# The output file will be stored in the same directory as the model file.
modelfile = f"{SNEWPY_models_base}/{modeltype}/{model}"
outfile = f"{modeltype}_{model}_{transformation}"

# There are three ways to select a time range.
# Option 1 - don't specify tstart and tend, then the whole model is integrated
# tstart = None
# tend = None

# Option 2 - specify single tstart and tend, this makes 1 fluence file integrated over the window
tstart = 0.0 * u.s
tend = 1.0 * u.s

# Option 3 = specify sequence of time intervals, one fluence file is made for each interval
# window_tstart = 0.0
# window_tend = 1.0
# window_bins = 10
# tstart = np.linspace(window_tstart, window_tend, window_bins, endpoint=False) * u.s
# tend = tstart + (window_tend - window_tstart) / window_bins * u.s

print("Preparing fluences ...")
tarredfile = snowglobes.generate_fluence(modelfile, modeltype, transformation, distance, outfile, tstart, tend)

print("Calculating events ...")
snowglobes.simulate(SNOwGLoBES_path, tarredfile, detector_input=detector)

print("Collating results ...")
tables = snowglobes.collate(SNOwGLoBES_path, tarredfile, skip_plots=True)

key = f"Collated_{outfile}_{detector}_events_smeared_weighted.dat"
# SNOwGLoBES also provides “unsmeared” events, without accounting for detector effects.
# Note: If you use Option 3 above (multiple time bins), `key` also includes the number of the time bin.

tables[key]

channels = tables[key]['header'].split()
energies = tables[key]['data'][0] * 1e3  # convert GeV to MeV
events = {}
for i, channel in enumerate(channels):
    if i==0: continue  # skip "Energy" array
    events[channel] = tables[key]['data'][i]
events['total'] = sum(tables[key]['data'][1:])

for channel in ['total'] + channels[1:]:
    plt.plot(energies, events[channel], label=channel)

plt.xlabel("Energy [MeV]")
plt.ylabel("Events per energy bin")
plt.yscale("log")
plt.xlim(0, 70)
plt.ylim(bottom=1)
plt.legend()
plt.savefig(f'plot_{outfile}_{detector}.pdf')

print(f"Total Events: {sum(events['total']):.3f}")
print(f"Mean Detected Energy: {sum(energies * events['total']) / sum(events['total']):.3f} MeV")