diff --git a/doc/scripts/PUSH_demo.py b/doc/scripts/PUSH_demo.py new file mode 100644 index 000000000..aa08a17d0 --- /dev/null +++ b/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)<,\n", + " 'eos': ['SFHo', 'SFHx', 'DD2', 'BHB', 'TM1', 'NL3'],\n", + " 'callibration': ['calI']}" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Wolfe_2023.param" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "54dbcde6-4aac-4f96-9dbe-ced31225919a", + "metadata": {}, + "outputs": [], + "source": [ + "# prepare the neutrino flux from the model\n", + "model = Wolfe_2023(progenitor_mass=27.6*u.solMass,eos='SFHo') # SN model\n", + "transformation = AdiabaticMSW(MixingParameters('NORMAL')) # Desired flavor transformation\n", + "\n", + "times = model.get_time()\n", + "energies = np.linspace(0,60,601)<" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "t = 50*u.ms\n", + "\n", + "ispec = model.get_initial_spectra(t, energies)\n", + "ospec_nmo = model.get_transformed_spectra(t, energies, transformation)\n", + "\n", + "fig, axes = plt.subplots(1,2, figsize=(12,5), sharex=True, sharey=True, tight_layout=True)\n", + "\n", + "for i, spec in enumerate([ispec, ospec_nmo]):\n", + " ax = axes[i]\n", + " plt.sca(ax)\n", + " spec.plot('energy')\n", + " \n", + " ax.set(title='Initial Spectra: $t = ${:.1f}'.format(t) if i==0 else 'Oscillated Spectra: $t = ${:.1f}'.format(t))\n", + " ax.grid()\n", + " ax.legend(loc='upper right', ncol=2, fontsize=16)\n", + "\n", + "ax = axes[0]\n", + "ax.set(ylabel=r'flux, MeV')\n", + "\n", + "fig.tight_layout();" + ] + }, + { + "cell_type": "markdown", + "id": "abc382c2-f06a-4a52-a829-07947991ec62", + "metadata": {}, + "source": [ + "## Using a detector config from SNOwGLoBES" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "0e7889dd-f442-4841-819b-ce7e31cd9a85", + "metadata": {}, + "outputs": [], + "source": [ + "from snewpy.rate_calculator import RateCalculator\n", + "\n", + "#load the RateCalculator object\n", + "rc = RateCalculator()" + ] + }, + { + "cell_type": "markdown", + "id": "6990252e-ab7f-4f6f-ad80-e86dd3c7794b", + "metadata": {}, + "source": [ + "### List available detectors" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "312c67d9-2556-4932-88b9-20ae7d042405", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['wc100kt30prct',\n", + " 'wc100kt15prct',\n", + " 'ar40kt',\n", + " 'scint20kt',\n", + " 'halo1',\n", + " 'halo2',\n", + " 'novaND',\n", + " 'novaFD',\n", + " 'wc100kt30prct_he',\n", + " 'ar40kt_he',\n", + " 'icecube',\n", + " 'km3net',\n", + " 'ds20',\n", + " 'argo',\n", + " 'lz',\n", + " 'xent',\n", + " 'pandax']" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#list available detectors\n", + "list(rc.detectors)" + ] + }, + { + "cell_type": "markdown", + "id": "389865fa-25e9-430a-9736-365a299fc7f9", + "metadata": {}, + "source": [ + "### Read the detector you need" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "c565f89f-e855-47b1-a591-da1561dd0bf5", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Detector(name=\"ar40kt\", mass=40.0 kt, channels=['nue_e', 'nuebar_e', 'numu_e', 'numubar_e', 'nutau_e', 'nutaubar_e', 'nue_Ar40', 'nuebar_Ar40', 'nc_nue_Ar40', 'nc_numu_Ar40', 'nc_nutau_Ar40', 'nc_nuebar_Ar40', 'nc_numubar_Ar40', 'nc_nutaubar_Ar40'])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#read the detector\n", + "detector = rc.read_detector('ar40kt')\n", + "detector" + ] + }, + { + "cell_type": "markdown", + "id": "8bfeb5a5-946b-4155-a72d-4adffdbf9b78", + "metadata": {}, + "source": [ + "### Inspecting the detector" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "a88a95fd-0e24-43e3-8205-41ac030a8632", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "{'nue_e': DetectionChannel (flavor=NU_E, smearing=True, weight=0.45),\n", + " 'nuebar_e': DetectionChannel (flavor=NU_E_BAR, smearing=True, weight=0.45),\n", + " 'numu_e': DetectionChannel (flavor=NU_MU, smearing=True, weight=0.45),\n", + " 'numubar_e': DetectionChannel (flavor=NU_MU_BAR, smearing=True, weight=0.45),\n", + " 'nutau_e': DetectionChannel (flavor=NU_TAU, smearing=True, weight=0.45),\n", + " 'nutaubar_e': DetectionChannel (flavor=NU_TAU_BAR, smearing=True, weight=0.45),\n", + " 'nue_Ar40': DetectionChannel (flavor=NU_E, smearing=True, weight=0.025),\n", + " 'nuebar_Ar40': DetectionChannel (flavor=NU_E_BAR, smearing=True, weight=0.025),\n", + " 'nc_nue_Ar40': DetectionChannel (flavor=NU_E, smearing=True, weight=0.025),\n", + " 'nc_numu_Ar40': DetectionChannel (flavor=NU_MU, smearing=True, weight=0.025),\n", + " 'nc_nutau_Ar40': DetectionChannel (flavor=NU_TAU, smearing=True, weight=0.025),\n", + " 'nc_nuebar_Ar40': DetectionChannel (flavor=NU_E_BAR, smearing=True, weight=0.025),\n", + " 'nc_numubar_Ar40': DetectionChannel (flavor=NU_MU_BAR, smearing=True, weight=0.025),\n", + " 'nc_nutaubar_Ar40': DetectionChannel (flavor=NU_TAU_BAR, smearing=True, weight=0.025)}" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#list all the channels\n", + "detector.channels" + ] + }, + { + "cell_type": "markdown", + "id": "11bfb80a-cef7-446f-83b7-eaa97b20c0d6", + "metadata": {}, + "source": [ + "### Running the rate calculation" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "948ca075-66eb-449b-b310-a9e599d1e2e9", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\jpknelle\\AppData\\Roaming\\Python\\Python312\\site-packages\\snewpy\\rate_calculator.py:353: RuntimeWarning: divide by zero encountered in log\n", + " return np.interp(np.log(E)/np.log(10), xp, yp, left=0, right=0)*E*1e-38 <" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "events = detector.run(fluence, detector_effects=True)\n", + "\n", + "plot_events(sum_events(list(events.values())), axis='energy', label='Total', lw=2, color='k')\n", + "for chan,numbers in events.items():\n", + " plot_events(numbers, axis='energy', label=chan)\n", + "\n", + "plt.legend(ncols=3)\n", + "plt.ylabel(f'Events per {numbers.energy.diff()[0]<bool: "Check if two Containers are equal" result = self.__class__==other.__class__ and \ diff --git a/python/snewpy/models/ccsn.py b/python/snewpy/models/ccsn.py index 253bf1b5a..bb3930067 100644 --- a/python/snewpy/models/ccsn.py +++ b/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) ` @@ -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 `_. @@ -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." diff --git a/python/snewpy/models/ccsn_loaders.py b/python/snewpy/models/ccsn_loaders.py index f52ee9eec..e2e383f76 100644 --- a/python/snewpy/models/ccsn_loaders.py +++ b/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 `_.""" 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) `_. """ @@ -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) `_. """ @@ -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) `_. """ @@ -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) + + + diff --git a/python/snewpy/models/model_files.yml b/python/snewpy/models/model_files.yml index a7126dce6..92ec5ac02 100644 --- a/python/snewpy/models/model_files.yml +++ b/python/snewpy/models/model_files.yml @@ -12,6 +12,12 @@ models: ccsn: Bollig_2016: repository: *ccsn_repository + + Bugli_2021: + repository: *ccsn_repository + + Fischer_2020: + repository: *ccsn_repository Fornax_2019: repository: *ccsn_repository @@ -27,6 +33,9 @@ models: Kuroda_2020: repository: *ccsn_repository + + Mori_2023: + repository: *ccsn_repository Nakazato_2013: repository: *ccsn_repository @@ -39,6 +48,9 @@ models: Sukhbold_2015: repository: *ccsn_repository + + Takata_2025: + repository: *ccsn_repository Tamborra_2014: repository: *ccsn_repository @@ -51,22 +63,14 @@ models: Warren_2020: repository: *ccsn_repository - - Zha_2021: - repository: *ccsn_repository - - Mori_2023: - repository: *ccsn_repository - Takata_2025: + Wolfe_2023: repository: *ccsn_repository - - Bugli_2021: - repository: *ccsn_repository - - Fischer_2020: + + Zha_2021: repository: *ccsn_repository + presn: Odrzywolek_2010: