Skip to content

Fornax_2019 fix #417

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
jpkneller wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

Fornax_2019 fix #417

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7305d7c
Fornax_2019 fix
jpkneller Nov 10, 2025
e75e17a
Add files via upload
jpkneller Nov 10, 2025
c713c71
Update ccsn_loaders.py
jpkneller Apr 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
232 changes: 102 additions & 130 deletions python/snewpy/models/ccsn_loaders.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -293,16 +293,14 @@ def __init__(self, filename, metadata={}, cache_flux=False):
cache_flux : bool
If true, pre-compute the flux on a fixed angular grid and store the values in a FITS file.
"""
#set the parameters
self.interpolation = "linear" #"linear"/"nearest"
self.phi = 0*u.deg #Input azimuth angles
self.theta = 0*u.deg #Input zenith angles

# Open the requested filename using the model downloader.
datafile = self.request_file(filename)

# Set up model metadata.
self.filename = filename
self.filename = os.path.basename(filename)
self.metadata = metadata

self.fluxunit = 1e50 * u.erg/(u.s*u.MeV)
self.dLdE_unit = 1e50 * u.erg/(u.s*u.MeV)
self.time = None

# Conversion of flavor to key name in the model HDF5 file.
Expand All @@ -312,89 +310,81 @@ def __init__(self, filename, metadata={}, cache_flux=False):
ThreeFlavor.NU_MU_BAR: 'nu2',
ThreeFlavor.NU_TAU: 'nu2',
ThreeFlavor.NU_TAU_BAR: 'nu2'}

self.E = {}
self.dE = {}
self.dLdE = {}
self.luminosity = {}

self.is_cached = False

# Read a cached flux file in FITS format or generate one.
self.is_cached = cache_flux and 'healpy' in sys.modules
logger = logging.getLogger()
if cache_flux and not 'healpy' in sys.modules:
logger = logging.getLogger()
logger.warning("No module named 'healpy'. Cannot enable caching.")

if self.is_cached:

self.E = {}
self.dE = {}
self.dLdE = {}
self.luminosity = {}

# Check if we're initializing on a FITS file or not.
if filename.endswith('.fits'):
fitsfile = filename
else:
fitsfile = filename.replace('h5', 'fits')
# Check if we're initializing on a FITS file or not.
if filename.endswith('.fits'):
fitsfile = filename
else:
fitsfile = filename.replace('h5', 'fits')

if os.path.exists(fitsfile):
self._read_fits(fitsfile)
ntim, nene, npix = self.dLdE[Flavor.NU_E].shape
self.npix = npix
self.nside = hp.npix2nside(npix)
else:
with h5py.File(filename, 'r') as _h5file:
if self.time is None:
self.time = _h5file['nu0']['g0'].attrs['time'] * u.s

# Use a HEALPix grid with nside=4 (192 pixels) to cache the
# values of Y_lm(theta, phi).
self.nside = 4
self.npix = hp.nside2npix(self.nside)
thetac, phic = hp.pix2ang(self.nside, np.arange(self.npix))

Ylm = {}
for l in range(3):
Ylm[l] = {}
for m in range(-l, l+1):
Ylm[l][m] = self._real_sph_harm(l, m, thetac, phic)

# Store 3D tables of dL/dE for each flavor.
logger = logging.getLogger()
for flavor in ThreeFlavor:

key = self._flavorkeys[flavor]
logger.info('Caching {} for {} ({})'.format(filename, str(flavor), key))

self.E[flavor] = _h5file[key]['egroup'][()] * u.MeV
self.dE[flavor] = _h5file[key]['degroup'][()] * u.MeV

ntim, nene = self.E[flavor].shape
self.dLdE[flavor] = np.zeros((ntim, nene, self.npix), dtype=float)
# Loop over time bins.
for i in range(ntim):
# Loop over energy bins.
for j in range(nene):
dLdE_ij = 0.
# Sum over multipole moments.
for l in range(3):
for m in range(-l, l+1):
dLdE_ij += _h5file[key]['g{}'.format(j)
]['l={} m={}'.format(l, m)][i] * Ylm[l][m]
self.dLdE[flavor][i][j] = dLdE_ij

# Integrate over energy to get L(t).
factor = 1. if flavor.is_electron else 0.25
self.dLdE[flavor] = self.dLdE[flavor] * factor * self.fluxunit
self.dLdE[flavor] = self.dLdE[flavor].to('erg/(s*MeV)')

self.luminosity[flavor] = np.sum(self.dLdE[flavor] * self.dE[flavor][:, :, np.newaxis], axis=1)
# Read a cached flux file in FITS format or generate one.
if cache_flux and os.path.exists(fitsfile):
self._read_fits(fitsfile)
ntim, nene, npix = self.dLdE[Flavor.NU_E].shape
self.npix = npix
self.nside = hp.npix2nside(npix)
self.is_cached = True
else:
with h5py.File(datafile, 'r') as _h5file:
if self.time is None:
self.time = _h5file['nu0']['g0'].attrs['time'] * u.s

# Use a HEALPix grid with nside=4 (192 pixels) to cache the
# values of Y_lm(theta, phi).
self.nside = 4
self.npix = hp.nside2npix(self.nside)
thetac, phic = hp.pix2ang(self.nside, np.arange(self.npix))

Ylm = {}
for l in range(3):
Ylm[l] = {}
for m in range(-l, l+1):
Ylm[l][m] = self._real_sph_harm(l, m, thetac, phic)

# Store 3D tables of dL/dE for each flavor.
for flavor in ThreeFlavor:
key = self._flavorkeys[flavor]
logger.info('Caching {} for {} ({})'.format(filename, str(flavor), key))

self.E[flavor] = _h5file[key]['egroup'][()] * u.MeV
self.dE[flavor] = _h5file[key]['degroup'][()] * u.MeV

ntim, nene = self.E[flavor].shape
self.dLdE[flavor] = np.zeros((ntim, nene, self.npix), dtype=float)
# Loop over time bins.
for i in range(ntim):
# Loop over energy bins.
for j in range(nene):
dLdE_ij = 0.
# Sum over multipole moments.
for l in range(3):
for m in range(-l, l+1):
dLdE_ij += _h5file[key]['g{}'.format(j)]['l={} m={}'.format(l, m)][i] * Ylm[l][m]
Copy link
Copy Markdown
Contributor

@mcolomerm mcolomerm Feb 17, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see why you integrate over the multiple moments / the armonics, instead of over directions. Moreover, this should be done always to the the angle-average as the output, right? i.e., whether the file is read as fits (if, with npix from file) or if it is read as h5 (else) and the npix is set to 4. Not just in the second case. Otherwise, the next steps for expected event rates etc will fail if I am not wrong.

self.dLdE[flavor][i][j] = dLdE_ij

factor = 1. if flavor.is_electron else 0.25
Copy link
Copy Markdown
Contributor

@mcolomerm mcolomerm Feb 17, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should make sure that this convention is followed for all the models, and agrees with the original model data. I think right now it is not the case from some model comparison tests I made: some groups provided it with alread the 0.25 factor applied, while others not.

self.dLdE[flavor] = self.dLdE[flavor] * factor * self.dLdE_unit

# Write output to FITS.
self._write_fits(fitsfile, overwrite=True)
else:
# Open the requested filename using the model downloader.
datafile = self.request_file(filename)
# Open HDF5 data file.
self._h5file = h5py.File(datafile, 'r')
if cache_flux:
self._write_fits(fitsfile, overwrite=True)
self.is_cached = True

# Integrate over energy to get L(t).
for flavor in ThreeFlavor:
self.luminosity[flavor] = np.sum(self.dLdE[flavor] * self.dE[flavor][:, :, np.newaxis], axis=1)

# Get grid of model times in seconds.
self.time = self._h5file['nu0']['g0'].attrs['time'] * u.s

def _read_fits(self, filename):
"""Read cached angular data from FITS.
Expand All @@ -421,7 +411,6 @@ def _read_fits(self, filename):
self.dLdE[flavor] = hdus[ext].data * u.Unit(hdus[ext].header['BUNIT'])
self.dLdE[flavor] = self.dLdE[flavor].to('erg/(s*MeV)')

self.luminosity[flavor] = np.sum(self.dLdE[flavor] * self.dE[flavor][:, :, np.newaxis], axis=1)

def _write_fits(self, filename, overwrite=False):
"""Write angular-dependent calculated flux in FITS format.
Expand Down Expand Up @@ -451,9 +440,9 @@ def _write_fits(self, filename, overwrite=False):
hdu_dE.header['BUNIT'] = 'MeV'
hx.append(hdu_dE)

hdu_flux = fits.ImageHDU(self.dLdE[flavor].to_value(str(self.fluxunit)))
hdu_flux = fits.ImageHDU(self.dLdE[flavor].to_value(str(self.dLdE_unit)))
hdu_flux.header['EXTNAME'] = '{}_FLUX'.format(name)
hdu_flux.header['BUNIT'] = str(self.fluxunit)
hdu_flux.header['BUNIT'] = str(self.dLdE_unit)
hx.append(hdu_flux)

hx.writeto(filename, overwrite=overwrite)
Expand Down Expand Up @@ -530,44 +519,20 @@ def _get_binnedspectra(self, t, theta, phi):
# Convert input time to a time index.
t = t.to(self.time.unit)
j = (np.abs(t - self.time)).argmin()
k = hp.ang2pix(self.nside, theta.to_value('radian'), phi.to_value('radian'))

for flavor in ThreeFlavor:
# Cached data: read out the relevant time and angular rows.
if self.is_cached:
# Convert input angles to a HEALPix index.
k = hp.ang2pix(self.nside, theta.to_value('radian'), phi.to_value('radian'))
E[flavor] = self.E[flavor][j]
dE[flavor] = self.dE[flavor][j]
binspec[flavor] = self.dLdE[flavor][j, :, k]

# Read the HDF5 input file directly and extract the spectra.
else:
key = self._flavorkeys[flavor]

# Energy binning of the model for this flavor, in units of MeV.
E[flavor] = self._h5file[key]['egroup'][j] * u.MeV
dE[flavor] = self._h5file[key]['degroup'][j] * u.MeV

# Storage of differential flux per energy, angle, and time.
dLdE = np.zeros(len(E[flavor]), dtype=float)

# Loop over energy bins.
for ebin in range(len(E[flavor])):
dLdE_j = 0
# Sum over multipole moments.
for l in range(3):
for m in range(-l, l + 1):
Ylm = self._real_sph_harm(l, m, theta.to_value('radian'), phi.to_value('radian'))
dLdE_j += self._h5file[key]['g{}'.format(ebin)]['l={} m={}'.format(l, m)][j] * Ylm
dLdE[ebin] = dLdE_j

factor = 1. if flavor.is_electron else 0.25
binspec[flavor] = dLdE * factor * self.fluxunit
binspec[flavor] = binspec[flavor].to('erg/(s*MeV)')
E[flavor] = self.E[flavor][j]
dE[flavor] = self.dE[flavor][j]
binspec[flavor] = self.dLdE[flavor][j, :, k]

return E, dE, binspec

def get_initial_spectra(self, t, E, theta, phi, flavors=ThreeFlavor, interpolation='linear'):
spectra_dict = self._get_initial_spectra_dict(t, E, theta, phi, flavors, interpolation)
return Spectrum.from_dict(spectra_dict,t, E,flavor_scheme=flavors)

def _get_initial_spectra_dict(self, t, E, flavors=ThreeFlavor):
def _get_initial_spectra_dict(self, t, E, theta, phi, flavors=ThreeFlavor, interpolation='linear'):
"""Get neutrino spectra/luminosity curves before flavor transformation.

Parameters
Expand All @@ -576,17 +541,24 @@ def _get_initial_spectra_dict(self, t, E, flavors=ThreeFlavor):
Time to evaluate initial spectra.
E : astropy.Quantity or ndarray of astropy.Quantity
Energies to evaluate the initial spectra.
theta : astropy.Quantity
Zenith angle of the spectral emission.
phi : astropy.Quantity
Azimuth angle of the spectral emission.
flavors: iterable of snewpy.neutrino.Flavor
Return spectra for these flavors only (default: all)
interpolation : str
Scheme to interpolate in spectra ('nearest', 'linear').

Returns
-------
initialspectra : dict
initial_spectra : dict
Dictionary of model spectra, keyed by neutrino flavor.
"""
initialspectra = {}
initial_spectra = {}

# Extract the binned spectra for the input t, theta, phi:
_E, _dE, _spec = self._get_binnedspectra(t, self.theta, self.phi)
_E, _dE, _spec = self._get_binnedspectra(t, theta, phi)

# Avoid "division by zero" in retrieval of the spectrum.
E[E == 0] = np.finfo(float).eps * E.unit
Expand All @@ -595,21 +567,21 @@ def _get_initial_spectra_dict(self, t, E, flavors=ThreeFlavor):
for flavor in flavors:

# Linear interpolation in flux.
if self.interpolation.lower() == 'linear':
if interpolation.lower() == 'linear':
# Pad log(E) array with values where flux is fixed to zero.
_logE = np.log10(_E[flavor].to_value('MeV'))
_dlogE = np.diff(_logE)
_logEbins = np.insert(_logE, 0, np.log10(np.finfo(float).eps * E.unit/u.MeV))
_logEbins = np.append(_logEbins, _logE[-1] + _dlogE[-1])

# Pad with values where flux is fixed to zero.
_dLdE = _spec[flavor].to_value(self.fluxunit)
_dLdE = _spec[flavor].to_value(self.dLdE_unit)
_dLdE = np.insert(_dLdE, 0, 0.)
_dLdE = np.append(_dLdE, 0.)

initialspectra[flavor] = np.interp(logE, _logEbins, _dLdE) * self.fluxunit
initial_spectra[flavor] = np.interp(logE, _logEbins, _dLdE) * self.dLdE_unit / E

elif self.interpolation.lower() == 'nearest':
elif interpolation.lower() == 'nearest':
_logE = np.log10(_E[flavor].to_value('MeV'))
_dlogE = np.diff(_logE)[0]
_logEbins = _logE - _dlogE
Expand All @@ -620,14 +592,14 @@ def _get_initial_spectra_dict(self, t, E, flavors=ThreeFlavor):
select = (idx > 0) & (idx < len(_E[flavor]))

_dLdE = np.zeros(len(E))
_dLdE[np.where(select)] = np.asarray([_spec[flavor][i].to_value(self.fluxunit) for i in idx[select]])
initialspectra[flavor] = _dLdE * self.fluxunit
_dLdE[np.where(select)] = np.asarray([_spec[flavor][i].to_value(self.dLdE_unit) for i in idx[select]])

initial_spectra[flavor] = _dLdE * self.dLdE_unit / E

else:
raise ValueError('Unrecognized interpolation type "{}"'.format(self.interpolation))

return initialspectra
raise ValueError('Unrecognized interpolation type "{}"'.format(interpolation))

return initial_spectra

class Fornax_2021(SupernovaModel):
def __init__(self, filename, metadata={}):
Expand Down
Loading