Detector Hits + Oscillations 2 (Obsolete)

Demonstrate the total number of hits in the IceCube detector after implementaing neutrino oscillations for a progentitor at distance 1 kpc. The following data is saved to a FITS file:

1. Luminosity corresponding to each neutrino flavor before mixing

2. Neutrino flux for each flavor in each oscillation scenario (unmixed, normal mixing, and inverted mixing)

3. Signal generated for each flavor in each oscillation scenario

The progenitor model used is the Nakazato Shen model with z = 0.02, t\(_{rev}\) = 300ms, m = 13.0M\(_\odot\)

from astropy import units as u
from astropy.table import Table

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

from asteria import config, source, detector
from asteria.interactions import Interactions
from asteria.neutrino import Flavor
from asteria.oscillation import SimpleMixing
import asteria.IO as io

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

Load Configuration

This will load the source configuration from a file.

For this to work, either the user needs to have done one of two things:

1. Run python setup.py install in the ASTERIA directory.

2. Run python setup.py develop and set the environment variable ASTERIA to point to the git source checkout.

If these were not done, the initialization will fail because the paths will not be correctly resolved.

conf = config.load_config('../../data/config/nakazato-shen-z0.02-t_rev300ms-s13.0.yaml')
ccsn = source.initialize(conf)
ic86 = detector.initialize(conf)

Prepare Iterables

Define the range of neutrino energies (E_nu) to simulate and the times (time) and distances (dist) at which to perform the simulation.

The progenitor distance is set to 1 kpc.

E_min = 0.1; E_max = 100.1; dE = 0.1;
Enu = np.arange(E_min, E_max, dE) * u.MeV

t_min = -1; t_max = 15; dt = 0.0001;
time = np.arange(t_min, t_max, dt) * u.s

ccsn.progenitor_distance = 1 * u.kpc

Redefine photonic_energy_per_volume

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).

def photonic_energy_per_vol(source, time, E, flavor, photon_spectrum, flux_list, n=1000):
        """Compute the energy deposited in a cubic meter of ice by photons
        from SN neutrino interactions.

        Parameters
        ----------

        time : float (units s)
            Time relative to core bounce.
        E : `numpy.ndarray`
            Sorted grid of neutrino energies
        flavor : :class:`asteria.neutrino.Flavor`
            Neutrino flavor.
        photon_spectrum : `numpy.ndarray` (Units vary, m**2)
            Grid of the product of lepton cross section with lepton mean energy
            and lepton path length per MeV, sorted according to parameter E
        flux_list: `numpy.ndarray`
            List containing neutrino fluxes corresponding to nu_e, nu_e_bar, nu_x, nu_x_bar
        n : int
            Maximum number of time steps to compute at once. A temporary numpy array
            of size n x time.size is created and can be very memory inefficient.

        Returns
        -------
        E_per_V
            Energy per m**3 of ice deposited  by neutrinos of requested flavor
        """
        H2O_in_ice = 3.053e28 # 1 / u.m**3

        t = time.to(u.s).value
        Enu = E.to(u.MeV)
        if Enu[0] == 0:
            Enu[0] = 1e-10 * u.MeV
        phot = photon_spectrum.to(u.m**2).value.reshape((-1,1)) # m**2

        dist = source.progenitor_distance.to(u.m).value # m**2
#         flux = source.get_flux( time, flavor ) # Unitless
        flux = flux_list
        if not flavor.is_electron:
            flux *= 2

        print('Beginning {0} simulation....'.format(flavor._name_), end='')
        # The following two lines exploit the fact that astropy quantities will
        # always return a number when numpy size is called on them, even if it is 1.
        E_per_V =  np.zeros( time.size )
        for i_part in source.parts_by_index(time, n): # Limits memory usage
            E_per_V[i_part] += np.trapz( source.energy_spectrum(t[i_part], Enu, flavor) * phot, Enu.value, axis=0)
        E_per_V *= flux * H2O_in_ice / ( 4 * np.pi * dist**2)
        print('Completed')

        return E_per_V * u.MeV / u.m**3

Compute Charged Particle Spectrum

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

The interactions are as follows:

• InvBetaTab : Tabulated inverse beta decay computation by Strumia and Vissani, Phys. Lett. B 564:42, 2003.

  • See Also: InvBetaPar : Inverse beta decay parameterization

• ElectronScatter : Elastic Neutrino-electron scattering from Marciano and Parsa, J. Phys. G 29:2969, 2003.

• Oxygen16CC : \(\nu\)-\(^{16}O\) charged current interaction, using estimates from Kolbe et al., PRD 66:013007, 2002.

• Oxygen16NC : \(\nu\)-\(^{16}O\) neutral current interaction, using estimates from Kolbe et al., PRD 66:013007, 2002.

• Oxygen18 : \(\nu\)-\(^{18}O\) charged current interaction, using estimates from Kamiokande data from Haxton and Robertson, PRC 59:515, 1999.

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\)).

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.

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.

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.

photon_spectra = np.zeros( shape=(len(Flavor), Enu.size) )

for nu, flavor in enumerate(Flavor):
    for interaction in Interactions:
        xs    = interaction.cross_section(flavor, Enu).to(u.m**2).value
        E_lep = interaction.mean_lepton_energy(flavor, Enu).value
        photon_scaling_factor = interaction.photon_scaling_factor(flavor).to( 1/u.MeV).value
        photon_spectra[nu] +=  xs * E_lep * photon_scaling_factor # u.m**2
photon_spectra *= u.m**2

Define Helper Function

def generate_signal(nu_list):
    '''Compute the number of hits generated in the detector for a given list of fluxes.
    The signal generated is by default binned to 1e-4s, i.e. the `time` binning.

    Parameter
    ---------
    nu_list: ndarray
        List of neutrino flavors, ordered by default as nu_e, nu_e_bar, nu_x, nu_x_bar
    timebin: float
        The desired time binning used to scale the deadtime fractions

    Returns
    -------
    hits: ndarray
        Number of hits in the detector
    '''
    signal = np.zeros(shape = (len(Flavor), time.size))

    for nu, (flavor, photon_spectrum) in enumerate(zip(Flavor, photon_spectra)):
        signal[nu] = photonic_energy_per_vol(ccsn, time, Enu, flavor, photon_spectrum, nu_list[nu])

    return signal

Implementing Neutrino Oscillations

Fluxes for each flavor are calculted for normal and inverted mixing using the class SimpleMixing in the neutrino.py module.

\(\theta_{12} = 33.82^{\circ}\) is used to initialize an object mix of SimpleMixing.

nu_list = np.zeros(shape = (4, 160000))
mix = SimpleMixing(33.82)

for i, flavor in enumerate(Flavor):
        nu_list[i] = ccsn.get_flux(time, flavor)

nu_norm = mix.normal_mixing(nu_list)
nu_inv = mix.inverted_mixing(nu_list)

[u_e_f, u_e_bar_f, u_x_f, u_x_bar_f] = nu_list
[n_e_f, n_e_bar_f, n_x_f, n_x_bar_f] = nu_norm
[i_e_f, i_e_bar_f, i_x_f, i_x_bar_f] = nu_inv

sig_unmixed = generate_signal(nu_list)

Beginning nu_e simulation....Completed
Beginning nu_e_bar simulation....Completed
Beginning nu_x simulation....Completed
Beginning nu_x_bar simulation....Completed

[u_e_s, u_e_bar_s, u_x_s, u_x_bar_s] = sig_unmixed

sig_norm = generate_signal(nu_norm)

Beginning nu_e simulation....Completed
Beginning nu_e_bar simulation....Completed
Beginning nu_x simulation....Completed
Beginning nu_x_bar simulation....Completed

[n_e_s, n_e_bar_s, n_x_s, n_x_bar_s] = sig_norm

sig_inv = generate_signal(nu_inv)

Beginning nu_e simulation....Completed
Beginning nu_e_bar simulation....Completed
Beginning nu_x simulation....Completed
Beginning nu_x_bar simulation....Completed

[i_e_s, i_e_bar_s, i_x_s, i_x_bar_s] = sig_inv

lum = np.zeros(shape = (len(Flavor), len(time)))
for nu, flavor in enumerate(Flavor):
    lum[nu] = ccsn.get_luminosity(time, flavor)
[e_lum, e_bar_lum, x_lum, x_bar_lum] = lum

# Create astropy table containing all relevant data

t = Table([time,
          e_lum, e_bar_lum, x_lum, x_bar_lum,
          u_e_f, u_e_bar_f, u_x_f, u_x_bar_f,
          n_e_f, n_e_bar_f, n_x_f, n_x_bar_f,
          i_e_f, i_e_bar_f, i_x_f, i_x_bar_f,
          u_e_s, u_e_bar_s, u_x_s, u_x_bar_s,
          n_e_s, n_e_bar_s, n_x_s, n_x_bar_s,
          i_e_s, i_e_bar_s, i_x_s, i_x_bar_s],
         names = ('time',
                 'nu_e luminosity', 'nu_e_bar luminosity','nu_x luminosity','nu_x_bar luminosity',
                 'unmixed nu_e flux', 'unmixed nu_e_bar flux', 'unmixed nu_x flux', 'unmixed nu_x_bar flux',
                 'normal nu_e flux', 'normal nu_e_bar flux', 'normal nu_x flux', 'normal nu_x_bar flux',
                 'inverted nu_e flux', 'inverted nu_e_bar flux', 'inverted nu_x flux', 'inverted nu_x_bar flux',
                 'unmixed nu_e sig', 'unmixed nu_e_bar sig', 'unmixed nu_x sig', 'unmixed nu_x_bar sig',
                 'normal nu_e sig', 'normal nu_e_bar sig', 'normal nu_x sig', 'normal nu_x_bar sig',
                 'inverted nu_e sig', 'inverted nu_e_bar sig', 'inverted nu_x sig', 'inverted nu_x_bar sig'
                 ))

t.write('nakazato-shen-z0.02-t_rev300ms-s13.0_e_per_v.fits', format = 'fits', overwrite = True)