PUSH model support

doc/scripts/PUSH_demo.py

@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+from snewpy.rate_calculator import RateCalculator
+from snewpy.models import ccsn, ccsn_loaders
+from snewpy.flavor_transformation import AdiabaticMSW
+from snewpy.neutrino import MixingParameters
+
+import astropy.units as u
+import numpy as np
+
+path = "~/.astropy/cache/snewpy/models/PUSH/"
+
+model = ccsn_loaders.PUSH("luminosity.d","mutau_luminosity.d") # SN model
+transformation = AdiabaticMSW(MixingParameters('NORMAL')) # Desired flavor transformation
+
+times    = model.get_time()
+energies = np.linspace(0,100,501)<<u.MeV
+distance = 10*u.kpc
+
+flux = model.get_flux(t=times, E=energies, distance=distance, flavor_xform=transformation)
+fluence = flux.integrate('time')
+
+detector = "wc100kt30prct"
+rc = RateCalculator()
+
+events = rc.run(fluence, detector, detector_effects=False)        
+events_smeared = rc.run(fluence, detector, detector_effects=True)
+
+# Compute number of events in all interaction channels
+total_events  = sum([chan.integrate_or_sum('energy').array.squeeze().value for chan in events.values()])        
+total_events_smeared  = sum([chan.integrate_or_sum('energy').array.squeeze().value for chan in events_smeared.values()])
+
+print("Total events in Super-K-like detector (with smearing):" , 0.32*total_events_smeared)
+

python/snewpy/flux.py

@@ -390,7 +390,26 @@ def _load_quantity(name):
             return cls(data=array,
                        **{name:_load_quantity(name) for name in ['time','energy','flavor']},
                        integrable_axes=f['_integrable_axes'])
-
+
+    def __add__(self,other:'Container'):
+        # Overload the + operator. 
+        # Don't compare the flavors, only that they have the same number
+        if self.__class__==other.__class__ and \
+            self.unit == other.unit and \
+            self.flavor_scheme==other.flavor_scheme and \
+            len(self.flavor)==len(other.flavor) and \
+            all([np.allclose(self.axes[ax], other.axes[ax]) for ax in list(Axes)[1:]]):
+                array = self.array+other.array
+                axes = list(self.axes)
+                return Container(array,*axes)
+        else:
+            return NotImplemented
+
+    def __radd__(self,other):
+        if other == 0:
+            return self
+        return self.__add__(other)    
+
     def __eq__(self, other:'Container')->bool:
         "Check if two Containers are equal"
         result = self.__class__==other.__class__ and \

python/snewpy/models/ccsn.py

@@ -34,14 +34,52 @@ class method to get a list of all valid combinations and filter it:
 from astropy import units as u
 from astropy.table import Table
 
+from snewpy.models import base
 from snewpy.models import ccsn_loaders as loaders
-from .base import PinchedModel
 
 from snewpy.models.registry_model import RegistryModel, Parameter
 from snewpy.models.registry_model import all_models
 from textwrap import dedent
 
+@RegistryModel(
+    progenitor_mass= [11.2, 27.] * u.Msun,
+    eos = ['SFHo'],
+    callibration = ['calI', 'calII'] 
+)
+class Ebinger_2018(loaders.PUSHArchiveModel):
+    """Model from the PUSH collaboration described in Ebinger et al.
+    """
+    def __init__(self, progenitor_mass:u.Quantity, eos:str='SFHo', callibration:str='calI'):
+        filename = f's{progenitor_mass.value:2.1f}_{eos}_{callibration}_Ebinger_luminosity.h5'
+        return super().__init__(filename=filename, metadata=self.metadata)
 
+@RegistryModel(
+    progenitor_mass= [11.2, 27.] * u.Msun,
+    eos = ['SFHo'],
+    callibration = ['calI', 'calII'] 
+)
+class Curtis_2019(loaders.PUSHArchiveModel):
+    """Model from the PUSH collaboration described in Curtis et al.
+    """
+    def __init__(self, progenitor_mass:u.Quantity, eos:str='SFHo', callibration:str='calI'):
+        filename = f's{progenitor_mass.value:2.1f}_{eos}_{callibration}_Curtis_luminosity.h5'
+        return super().__init__(filename=filename, metadata=self.metadata)
+
+
+@RegistryModel(
+    #progenitor_mass = np.concat( (np.arange(10.8,28.2+0.01,0.2),np.arange(29,40+0.01,1)) ) * u.Msun,
+    progenitor_mass = [ 10.8, 27.6, 28.2, 29, 40] * u.Msun,
+    eos = ['SFHo', 'SFHx', 'DD2', 'BHB', 'TM1', 'NL3'],
+    callibration = ['calI'] 
+)
+class Wolfe_2023(loaders.PUSHArchiveModel):
+    """Model from the PUSH collaboration described in Wolfe et al.
+    """
+    def __init__(self, progenitor_mass:u.Quantity, eos:str='SFHo', callibration:str='calI'):
+        filename = f's{progenitor_mass.value:2.1f}_{eos}_{callibration}_Wolfe_luminosity.h5'
+        return super().__init__(filename=filename, metadata=self.metadata)
+
+
 @RegistryModel()
 class Fischer_2020(loaders.Fischer_2020):
     """Model based on simulations from `Fischer et al. (2020) <https://arxiv.org/abs/1804.10890>`
@@ -52,12 +90,12 @@ def __init__(self):
         filename='Fischer_2020.tar.gz'
         return super().__init__(filename, metadata=self.metadata)
 
+
 @RegistryModel(
     progenitor_mass = [13, 20, 30, 50] * u.Msun,
     revival_time = [0, 100, 200, 300] * u.ms,
     metallicity = [0.02, 0.004],
     eos = ['LS220', 'shen', 'togashi'],
-
     _param_validator = lambda p: (p['revival_time'] == 0 * u.ms and p['progenitor_mass'] == 30 * u.Msun
                                   and p['metallicity'] == 0.004) or \
                                  (p['revival_time'] != 0 * u.ms and p['eos'] == 'shen'
@@ -190,7 +228,7 @@ def __init__(self, progenitor_mass:u.Quantity):
 @RegistryModel(
     progenitor_mass = Parameter(values=(list(range(16, 27)) + [19.89, 22.39, 30, 33]) * u.Msun,
                                 desc_values = '[16..26, 19.89, 22.39, 30, 33] solMass'
-                               ),
+                               )
 )
 class Zha_2021(loaders.Zha_2021):
     """Model based on the hadron-quark phse transition models from `Zha et al. 2021 <https://arxiv.org/abs/2103.02268>`_.
@@ -212,7 +250,7 @@ def __init__(self, *, progenitor_mass:u.Quantity):
                               name='turbmixing_param',
                               label='Turb. mixing param.',
                               description='Turbulent mixing parameter alpha_lambda',
-                              ),
+                              )
 )
 class Warren_2020(loaders.Warren_2020):
     """Model based on simulations from Warren et al., ApJ 898:139, 2020.
@@ -535,9 +573,10 @@ def get_fluence(self, t):
 
         return fluence
 
-class Analytic3Species(PinchedModel):
-    """An analytical model calculating spectra given total luminosity,
-    average energy, and rms or pinch, for each species.
+class Analytic3Species(base.PinchedModel):
+    """This is the basically the loader version of base.PinchedModel i.e. it reads the data 
+    for the PinchedModel from a file. The format of the file is that made by 
+    the `doc/scripts/Analytic.py` script
     """
 
     param = "There are no input files available for this class. Use `doc/scripts/Analytic.py` in the SNEWPY GitHub repo to create a custom input file."

python/snewpy/models/ccsn_loaders.py

@@ -9,6 +9,7 @@
 import re
 import sys
 import tarfile
+from pathlib import Path
 
 from astropy import units as u
 from astropy.table import Table, join
@@ -19,14 +20,14 @@
 import numpy as np
 from scipy.special import gamma, lpmv
 
+from snewpy.models import base
 from snewpy.flux import Spectrum
-from snewpy.models.base import PinchedModel, SupernovaModel
 from snewpy.flavor import ThreeFlavor
 from snewpy import _model_downloader
 
 import multiprocessing
 
-class GarchingArchiveModel(PinchedModel):
+class GarchingArchiveModel(base.PinchedModel):
     """Subclass that reads models in the format used in the
     `Garching Supernova Archive <https://wwwmpa.mpa-garching.mpg.de/ccsnarchive/>`_."""
     def __init__(self, filename, eos='LS220', metadata={}):
@@ -88,7 +89,63 @@ def __init__(self, filename, eos='LS220', metadata={}):
             }
         super().__init__(simtab, metadata)
 
-class Nakazato_2013(PinchedModel):
+
+class PUSHArchiveModel(base.PinchedModel):
+    """Subclass that reads models in the format used
+    by the PUSH collaboration
+    """
+
+    def __init__(self, filename, metadata={}):
+        """
+        Parameters
+        ----------
+        filename : str
+            Absolute or relative path to model data
+        """
+        datafile = self.request_file(filename) 
+        f = h5py.File(datafile, 'r')
+
+        simtab = Table()
+
+        tbounce = f['metadata'].attrs['bounce_time'] * u.s
+
+        data = np.array(f['data'])
+
+        # Keep row only if all elements are >= 0
+        columns = data[:, 1:]
+        mask = np.any(columns<=0,axis=1)
+        data = data[~mask]
+
+        simtab['TIME'] = data[:,0] * u.s - tbounce
+
+        simtab['L_NU_E'] = data[:,3] << u.erg/u.s
+        simtab['L_NU_E_BAR'] = data[:,4] << u.erg/u.s
+        simtab['L_NU_X'] = data[:,6] << u.erg/u.s
+
+        simtab['E_NU_E'] = data[:,3] / data[:,1] << u.erg
+        simtab['E_NU_E_BAR'] = data[:,4] / data[:,2] << u.erg
+        simtab['E_NU_X'] = data[:,6] / data[:,5] << u.erg
+
+        simtab['ALPHA_NU_E'] = np.full(len(simtab['TIME']),3)
+        simtab['ALPHA_NU_E_BAR'] = simtab['ALPHA_NU_E']
+        simtab['ALPHA_NU_X'] = simtab['ALPHA_NU_E']
+
+        # prevent negative luminosities
+        simtab['L_NU_E'][simtab['L_NU_E'] < 0] = 1
+        simtab['L_NU_E_BAR'][simtab['L_NU_E_BAR'] < 0] = 1
+        simtab['L_NU_X'][simtab['L_NU_X'] < 0] = 1
+
+        metadata = dict(f['metadata'].attrs)
+
+        k = f['metadata/compactness_forpush'].attrs['columns']        
+        compacness = np.transpose( np.array( f['metadata/compactness_forpush'] ) )
+
+        metadata = metadata | dict(zip(k,compacness))
+
+        super().__init__(simtab, metadata)
+
+
+class Nakazato_2013(base.PinchedModel):
     def __init__(self, filename, metadata={}):
         """Model initialization.
 
@@ -131,7 +188,7 @@ class Walk_2019(GarchingArchiveModel):
     pass
 
 
-class OConnor_2013(PinchedModel):
+class OConnor_2013(base.PinchedModel):
     """Model based on the black hole formation simulation in `O'Connor & Ott (2013) <https://arxiv.org/abs/1207.1100>`_.
     """
 
@@ -163,7 +220,7 @@ def __init__(self, filename, metadata={}):
         super().__init__(simtab, metadata)
 
 
-class OConnor_2015(PinchedModel):
+class OConnor_2015(base.PinchedModel):
     """Model based on the black hole formation simulation in `O'Connor (2015) <https://arxiv.org/abs/1411.7058>`_.
     """
 
@@ -209,7 +266,7 @@ def __init__(self, filename, metadata={}):
 class Zha_2021(OConnor_2015):
     pass
 
-class Warren_2020(PinchedModel):
+class Warren_2020(base.PinchedModel):
     def __init__(self, filename, metadata={}):
         """
         Parameters
@@ -255,7 +312,7 @@ def __init__(self, filename, metadata={}):
         super().__init__(simtab, metadata)
 
 
-class Kuroda_2020(PinchedModel):
+class Kuroda_2020(base.PinchedModel):
     def __init__(self, filename, metadata={}):
         """
         Parameters
@@ -282,7 +339,7 @@ def __init__(self, filename, metadata={}):
 
         super().__init__(simtab, metadata)
 
-class Fornax_2019(SupernovaModel):
+class Fornax_2019(base.SupernovaModel):
     def __init__(self, filename, metadata={}, cache_flux=False):
         """
         Parameters
@@ -657,7 +714,7 @@ def _get_initial_spectra_dict(self, t, E, theta, phi, flavors=ThreeFlavor, inter
 
         return initial_spectra
 
-class Fornax_2021(SupernovaModel):
+class Fornax_2021(base.SupernovaModel):
     def __init__(self, filename, metadata={}):
         """
         Parameters
@@ -877,7 +934,7 @@ def __init__(self, filename, metadata={}):
             self.luminosity[flavor] = np.sum(dLdE*dE, axis=1) * factor * 1e50 * u.erg/u.s
 
 
-class Mori_2023(PinchedModel):
+class Mori_2023(base.PinchedModel):
     def __init__(self, filename, metadata={}):
         """
         Parameters
@@ -925,7 +982,7 @@ def __init__(self, filename, metadata={}):
         super().__init__(simtab, metadata)
 
 
-class Takata_2025(PinchedModel):
+class Takata_2025(base.PinchedModel):
     def __init__(self, filename, metadata={}):
         """
         Parameters
@@ -971,7 +1028,7 @@ def __init__(self, filename, metadata={}):
         super().__init__(simtab, metadata)
 
 
-class Bugli_2021(PinchedModel):
+class Bugli_2021(base.PinchedModel):
     """Model based on `Buggli (2021) <https://arxiv.org/abs/2105.00665>`_.
     """
 
@@ -1002,7 +1059,7 @@ def __init__(self, filename, metadata={}):
         super().__init__(simtab, metadata)
 
 
-class Fischer_2020(PinchedModel):
+class Fischer_2020(base.PinchedModel):
     def __init__(self, filename, metadata={}):
         """
         Parameters
@@ -1069,3 +1126,6 @@ def __init__(self, filename, metadata={}):
         tf.close()
 
         super().__init__(simtab, metadata)
+
+
+