Skip to content

Commit

Permalink
Merge pull request #330 from SNEWS2/Sheshuk/add_Yoshida_model
Browse files Browse the repository at this point in the history 
Add presn.Yoshida_2016 model
  • Loading branch information
@JostMigenda
JostMigenda committed May 31, 2024
2 parents f5281c2 + e57875a commit a97833c
Show file tree
Hide file tree
Showing 5 changed files with 337 additions and 2 deletions.
287 changes: 287 additions & 0 deletions doc/nb/presn/Yoshida_2016.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,287 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "83d2d170-a5b3-44d9-ad05-762448b08c80",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from astropy import units as u\n",
"\n",
"#import the snewpy modules\n",
"from snewpy.models import presn\n",
"from snewpy.flavor_transformation import NoTransformation, AdiabaticMSW, MassHierarchy\n",
"from snewpy.neutrino import Flavor"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f2ad86ec-ff7e-4ef6-818e-f244ac8112f4",
"metadata": {},
"outputs": [],
"source": [
"# set the plot parameters\n",
"plt.rc('grid', ls=':')\n",
"plt.rc('axes', grid=True)\n",
"plt.rc('legend', fontsize=12, loc='upper right')\n",
"\n",
"# define drawing styles for each flavor\n",
"styles = {f: dict(color='C0' if f.is_electron else 'C1',\n",
" ls='-' if f.is_neutrino else ':',\n",
" label=f.to_tex()) for f in Flavor}"
]
},
{
"cell_type": "markdown",
"id": "4f1b6439-6702-417a-a9b6-0ca0d5798f47",
"metadata": {},
"source": [
"## Initialize the model and calculate the flux"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ffff2732-a87a-401c-b85d-a0f0d5598657",
"metadata": {},
"outputs": [],
"source": [
"# See what progenitors are available …\n",
"print(f\"Available progenitors:\\n{presn.Yoshida_2016.param}\")\n",
"\n",
"# … and select one of them to initialise the model\n",
"model = presn.Yoshida_2016(progenitor_mass=15*u.Msun)\n",
"model"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4d8e7b15-0a7e-49f3-b6bd-9463e8b371a5",
"metadata": {},
"outputs": [],
"source": [
"# Energy array and time to compute spectra.\n",
"# Note that any convenient units can be used and the calculation will remain internally consistent.\n",
"E = np.linspace(0, 25, 201) * u.MeV\n",
"t = np.geomspace(-2*u.day, -1*u.s, 101)\n",
"distance = 1*u.kpc"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "efb466d7-1b4c-4949-a983-5134fba3f35b",
"metadata": {},
"outputs": [],
"source": [
"flux = model.get_flux(t, E,\n",
" distance=distance,\n",
" flavor_xform=NoTransformation())\n",
"flux"
]
},
{
"cell_type": "markdown",
"id": "f06fc5ed-e39f-4f91-a8c6-09a671a7f6d0",
"metadata": {},
"source": [
"## Plotting the integral neutrino fluence and rates"
]
},
{
"cell_type": "markdown",
"id": "34cd5637-201f-493c-9aa8-10a653d4158b",
"metadata": {},
"source": [
"### Integral neutrino fluence vs. Energy for the last hour before collapse"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4f8f5b53-257d-43b7-b933-415738e8b4a1",
"metadata": {},
"outputs": [],
"source": [
"#integrate the flux over the last hour before the collapse\n",
"fluence = flux.integrate('time', limits=[-1, 0]<<u.hour)\n",
"#get the relevant arrays for plotting\n",
"x = fluence.energy\n",
"y = fluence.array.squeeze()\n",
"#plot each flavor\n",
"for flv in fluence.flavor:\n",
" plt.plot(x,y[flv], **styles[flv])\n",
"#add the legend\n",
"plt.legend(ncol=2, loc='lower center')\n",
"#adjust the scales\n",
"plt.autoscale(tight=True,axis='x')\n",
"plt.yscale('log')\n",
"plt.ylim(1)\n",
"#define the labels\n",
"xlabel = 'Neutrino energy'\n",
"ylabel = 'Neutrino fluence'\n",
"plt.xlabel(f'{xlabel} [{x.unit._repr_latex_()}]')\n",
"plt.ylabel(f'{ylabel} [{y.unit._repr_latex_()}]')\n",
"#add the plot title\n",
"plt.title('Neutrinos in the last 1 hour before collapse')\n",
"#add the model parameters\n",
"plt.annotate(str(model) + f'\\ndistance : {distance}',\n",
" xy=(0.98,0.98),\n",
" xycoords='axes fraction',\n",
" va='top', ha='right'\n",
" )\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "8347961e",
"metadata": {},
"source": [
"### Integral neutrino rate vs. time"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dba0da41-6061-4c7a-b1fc-4b73f29ea960",
"metadata": {},
"outputs": [],
"source": [
"#integrate the flux over all energy bins\n",
"rate = flux.integrate('energy')\n",
"#get the relevant arrays for plotting\n",
"x = rate.time\n",
"y = rate.array.squeeze()\n",
"#we can also convert to a more convenient unit\n",
"y = y.to('1/(hour*m**2)')\n",
"#plot each flavor\n",
"for flv in rate.flavor:\n",
" plt.plot(x,y[flv], **styles[flv])\n",
"#add the legend\n",
"plt.legend(ncol=2, loc='lower center')\n",
"#adjust the scales\n",
"plt.autoscale(tight=True,axis='x')\n",
"plt.yscale('log')\n",
"#plt.ylim(1)\n",
"#define the labels\n",
"xlabel = 'Time before collapse'\n",
"ylabel = 'Neutrino rate'\n",
"plt.xlabel(f'{xlabel} [{x.unit._repr_latex_()}]')\n",
"plt.ylabel(f'{ylabel} [{y.unit._repr_latex_()}]')\n",
"#add the plot title\n",
"plt.title('Integral neutrino rate')\n",
"#add the model parameters\n",
"plt.annotate(str(model) + f'\\ndistance : {distance}',\n",
" xy=(0.02,0.98),\n",
" xycoords='axes fraction',\n",
" va='top', ha='left'\n",
" )\n",
"plt.xscale('symlog')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "2dfec9e5-2819-4993-b352-819224b65279",
"metadata": {},
"source": [
"## Plot energy spectrum with neutrino oscillations"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "94e4605b-fcee-4d3a-8344-caa4100e3ee1",
"metadata": {},
"outputs": [],
"source": [
"# Energy array and time to compute spectra.\n",
"# Note that any convenient units can be used and the calculation will remain internally consistent.\n",
"E = np.linspace(0,25,201) * u.MeV\n",
"t = -1*u.s\n",
"\n",
"xforms = [NoTransformation(), AdiabaticMSW(mh=MassHierarchy.NORMAL), AdiabaticMSW(mh=MassHierarchy.INVERTED)]\n",
"labels = ['No transformation', 'AdiabaticMSW NMO', 'AdiabaticMSW IMO']\n",
"xform_styles = ['-','--',':']\n",
"flavor_colors = ['C0','C1','C2','C3']\n",
"\n",
"#calculate the flux for different transformations\n",
"fluxes = {xform: model.get_flux(t, E, distance, flavor_xform=xform) for xform in xforms}\n",
"#make a big figure\n",
"fig, axes = plt.subplots(1,4,figsize=(16,4), sharex=True, sharey=True)\n",
"\n",
"for flv, ax in zip(Flavor,axes):\n",
" plt.sca(ax)\n",
" for xform, line_style, label in zip(xforms,xform_styles, labels):\n",
" flux = fluxes[xform]\n",
" #get the relevant arrays for plotting\n",
" x = flux.energy\n",
" y = flux.array[flv].squeeze()\n",
" #plot each flavor\n",
" plt.plot(x,y, \n",
" ls=line_style, \n",
" color=flavor_colors[flv], \n",
" label=f'{Flavor(flv).to_tex()} {label}'\n",
" )\n",
" \n",
" #add the legend\n",
" plt.legend(ncol=1, loc='lower left', fontsize=10)\n",
" \n",
" #define the labels\n",
" xlabel = 'Neutrino energy'\n",
" ylabel = 'Neutrino flux'\n",
" plt.xlabel(f'{xlabel} [{x.unit._repr_latex_()}]')\n",
" #add the plot title\n",
"\n",
" #add the model parameters\n",
" plt.annotate(str(model) + f'\\ndistance : {distance}',\n",
" xy=(0.98,0.98),\n",
" xycoords='axes fraction',\n",
" va='top', ha='right'\n",
" )\n",
" axes[0].set_ylabel(f'{ylabel} [{y.unit._repr_latex_()}]')\n",
" #adjust the scales\n",
" plt.autoscale(tight=True,axis='x')\n",
" plt.autoscale(tight=False,axis='y')\n",
" plt.yscale('log')\n",
" plt.ylim(1)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d12483ed-e5f7-48a0-99a7-9f7f1bfdd289",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
4 changes: 3 additions & 1 deletion python/snewpy/models/model_files.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
config:
- &snewpy "https://github.com/SNEWS2/snewpy/raw/v{snewpy_version}/models/{model}/{filename}"
- &ccsn_repository "https://github.com/SNEWS2/snewpy-models-ccsn/raw/v0.2/models/{model}/{filename}"
- &presn_repository "https://github.com/SNEWS2/snewpy-models-presn/raw/v0.1/models/{model}/{filename}"
- &presn_repository "https://github.com/SNEWS2/snewpy-models-presn/raw/v0.2/models/{model}/{filename}"

models:
ccsn:
Expand Down Expand Up @@ -65,3 +65,5 @@ models:
Kato_2017:
repository: *presn_repository

Yoshida_2016:
repository: *presn_repository
14 changes: 14 additions & 0 deletions python/snewpy/models/presn.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,3 +45,17 @@ class Kato_2017(loaders.Kato_2017):
def __init__(self, progenitor_mass:u.Quantity):
path=f"pre_collapse/m{progenitor_mass.to_value('Msun'):.0f}"
super().__init__(path)


@RegistryModel(
progenitor_mass = [12, 15, 20]<<u.Msun
)
class Yoshida_2016(loaders.Yoshida_2016):
"""Presupernova model based on
`Yoshida et al. (2016), PRD 93, 123012 <https://doi.org/10.1103/PhysRevD.93.123012>`_
Dataset available on `Zenodo <https://zenodo.org/records/3778014>`__
"""
def __init__(self, progenitor_mass:u.Quantity):
path=f"t_spc_m{progenitor_mass.to_value('Msun'):.0f}_1.2.txt"
super().__init__(path)
32 changes: 32 additions & 0 deletions python/snewpy/models/presn_loaders.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,3 +139,35 @@ def get_initial_spectra(self, t, E, flavors=Flavor):
flux = self.interpolated(t, E) / (u.MeV * u.s)
return {f: flux[f] for f in flavors}

class Yoshida_2016(SupernovaModel):
"""Set up a presupernova model based on
[Yoshida et al. (2016), PRD 93, 123012.]
"""
def __init__(self, filename):
with open(self.request_file(filename)) as f:
data = []
T = []
while (line := f.readline()):
if not line: break
T += [float(line.split()[1])]
data += [[np.loadtxt(f, max_rows=100).flatten() for i in range(4)]]
times = np.array(T)
energies = np.concatenate([
np.linspace(0,10,1001)[1:],
np.linspace(10,20,501)[1:]
])
dNdEdT = np.stack(data, axis=1)
#rearrange flavors from NU_E, NU_E_BAR,_NU_X, NU_X_BAR to current
indices = np.argsort([Flavor.NU_E,Flavor.NU_E_BAR, Flavor.NU_X, Flavor.NU_X_BAR])
dNdEdT = dNdEdT.take(indices, axis=0)

self.interpolated = _interp_TE(
times, energies, dNdEdT, ax_t=1, ax_e=2
)
super().__init__(times << u.s, self.metadata)

def get_initial_spectra(self, t, E, flavors=Flavor):
t = np.array(-t.to_value("s"), ndmin=1)
E = np.array(E.to_value("MeV"), ndmin=1)
flux = self.interpolated(t, E) / (u.MeV * u.s)
return {f: flux[f] for f in flavors}
2 changes: 1 addition & 1 deletion python/snewpy/test/test_presn_rates.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
T = np.geomspace(-1*u.hour, -1*u.min,1000)
E = np.linspace(0,20,100)*u.MeV

@pytest.mark.parametrize('model_class',[presn.Odrzywolek_2010, presn.Kato_2017, presn.Patton_2017])
@pytest.mark.parametrize('model_class',[presn.Odrzywolek_2010, presn.Kato_2017, presn.Patton_2017, presn.Yoshida_2016])
@pytest.mark.parametrize('transformation',[AdiabaticMSW(mh=mh) for mh in MassHierarchy])
@pytest.mark.parametrize('detector', ["wc100kt30prct"])
def test_presn_rate(model_class, transformation, detector):
Expand Down

0 comments on commit a97833c

Please sign in to comment.