Skip to content

Commit

Permalink
Merge branch 'master' of github.com:IceCubeOpenSource/flarestack
Browse files Browse the repository at this point in the history
  • Loading branch information
robertdstein committed Aug 5, 2020
2 parents 1e94b6c + 0a6644d commit cb65458
Show file tree
Hide file tree
Showing 9 changed files with 173 additions and 68 deletions.
63 changes: 30 additions & 33 deletions flarestack/cosmo/icecube_diffuse_flux/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -107,51 +107,48 @@ def get_diffuse_flux_contour(fit="joint_15", contour_name="68"):

return best_f, upper_f, lower_f, e_range

def plot_diffuse_flux(label, contour_68, contour_95, e_range):
def plot_diffuse_flux(fit="joint_15"):
plt.figure()

best_fit_flux, best_fit_gamma, contour_68, contour_95, e_range = load_fit(fit)

plt.figure()
for i, contour in enumerate([68., 95.]):
all_c = [contour_68, contour_95][i]

# Plot 68% contour
best_f, upper_f, lower_f, e_range = get_diffuse_flux_contour(fit=fit, contour_name=contour)

plt.plot(e_range,
e_range ** 2 * best_f(e_range),
alpha=1.0, color="k")

for (index, norm) in contour_68:
plt.plot(e_range,
e_range ** 2 * flux_f(e_range, norm, index),
alpha=0.6)
plt.plot(
e_range,
e_range ** 2 * upper_contour(e_range, contour_68),
color="k", label="68% contour", linestyle="-"
)
plt.plot(
e_range,
e_range ** 2 * lower_contour(e_range, contour_68),
color="k", linestyle="-",
)

# Plot 95% contour

for (index, norm) in contour_95:
e_range ** 2 * lower_f(e_range),
color=f"C{i}",
label=f"{contour} %",
)

plt.plot(e_range,
e_range ** 2 * flux_f(e_range, norm, index),
alpha=0.3)
plt.plot(
e_range,
e_range ** 2 * upper_contour(e_range, contour_95),
color="k", linestyle="--", label="95% contour"
)
plt.plot(
e_range,
e_range ** 2 * lower_contour(e_range, contour_95),
color="k", linestyle="--"
)
e_range ** 2 * upper_f(e_range),
color=f"C{i}",
)

for (index, norm) in all_c:

plt.plot(e_range,
e_range**2 * flux_f(e_range, norm, index),
alpha=0.15
)

plt.yscale("log")
plt.xscale("log")
plt.legend()
plt.title(r"Diffuse Flux Global Best Fit ($\nu_{\mu} + \bar{\nu}_{\mu})$")
plt.ylabel(r"$E^{2}\frac{dN}{dE}$[GeV cm$^{-2}$ s$^{-1}$ sr$^{-1}$]")
plt.xlabel(r"$E_{\nu}$ [GeV]")
savepath = illustration_dir + "diffuse_flux_global_fit.pdf"
plt.tight_layout()
savepath = f"{illustration_dir}diffuse_flux_global_fit_{fit}.pdf"

logger.info(f"Saving to {savepath}")

plt.savefig(savepath)
plt.close()
2 changes: 1 addition & 1 deletion flarestack/cosmo/icecube_diffuse_flux/nt_16.py
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@
best_fit_flux,
best_fit_gamma,
contour_68,
contour_68,
contour_95,
e_range,
"https://arxiv.org/abs/1607.08006"
)
Expand Down
2 changes: 1 addition & 1 deletion flarestack/cosmo/icecube_diffuse_flux/nt_17.py
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@
best_fit_flux,
best_fit_gamma,
contour_68,
contour_68,
contour_95,
e_range,
"https://doi.org/10.22323/1.301.1005"
)
Expand Down
2 changes: 1 addition & 1 deletion flarestack/cosmo/icecube_diffuse_flux/nt_19.py
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@
best_fit_flux,
best_fit_gamma,
contour_68,
contour_68,
contour_95,
e_range,
"https://arxiv.org/abs/1908.09551"
)
Expand Down
9 changes: 6 additions & 3 deletions flarestack/cosmo/rates/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,21 +3,24 @@
from flarestack.cosmo.rates.ccsn_rates import get_ccsn_rate
from flarestack.cosmo.rates.tde_rates import get_tde_rate
from flarestack.cosmo.rates.grb_rates import get_grb_rate
from flarestack.cosmo.rates.fbot_rates import get_fbot_rate

logger = logging.getLogger(__name__)

source_maps = {
"tde": ["TDE", "tidal_disruption_event"],
"sfr": ["SFR", "star_formation_rate"],
"ccsn": ["CCSN", "sn", "supernova", "core_collapse_supernova"],
"grb": ["GRB", "gamma_ray_burst"]
"grb": ["GRB", "gamma_ray_burst"],
"fbot": ["FBOT", "fast_blue_optical_transient"]
}

sources = {
"tde": get_tde_rate,
"sfr": get_sfr_rate,
"ccsn": get_ccsn_rate,
"grb": get_grb_rate
"grb": get_grb_rate,
"fbot": get_fbot_rate
}

def get_rate(source_name, evolution_name=None, rate_name=None, fraction=1.0, **kwargs):
Expand Down Expand Up @@ -49,6 +52,6 @@ def get_rate(source_name, evolution_name=None, rate_name=None, fraction=1.0, **k
f = sources[source_name](evolution_name, rate_name, **kwargs)

if fraction != 1.0:
logger.info(f"Assuming a modified rate that is {100.*fraction:.2f} of that total.")
logger.info(f"Assuming a modified rate that is {100.*fraction:.2f}% of that total.")

return lambda z: f(z) * fraction
64 changes: 46 additions & 18 deletions flarestack/cosmo/rates/ccsn_rates.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,16 +6,40 @@

# Taken from https://arxiv.org/pdf/1509.06574.pdf

sn_types = {
"IIn": 0.064,
"IIP": 0.52,
"Ib": 0.069,
"Ic": 0.176,
"Ibc": 0.069 + 0.176,
"all": 1.0
sn_subclass_rates = {
"li_11": ({
"IIn": 0.064,
"IIP": 0.52,
"IIL": 0.073,
"II": 0.064 + 0.52 + 0.073,
"Ib": 0.069,
"Ic": 0.176,
"Ibc": 0.069 + 0.176,
}, "https://arxiv.org/abs/1006.4612")
}

def get_sn_fraction(sn_subclass=None):
def get_sn_subfraction(sn_subclass_fractions_name=None):
"""Return a value of kcc (SN per unit star formation)
:param sn_subclass_fractions_name: Name of kcc to be used
:return: Value of kcc
"""

if sn_subclass_fractions_name is None:
logging.info("No specified sn_subclass_fractions_name. Assuming default.")
sn_subclass_fractions_name = "li_11"

if sn_subclass_fractions_name not in sn_subclass_rates.keys():
raise Exception(f"Subclass name '{sn_subclass_fractions_name}' not recognised. "
f"The following kcc values are available: {sn_subclass_rates.keys()}")
else:
sn_rates, ref = sn_subclass_rates[sn_subclass_fractions_name]
logging.info(f"Loaded SN subclass fractions '{sn_subclass_fractions_name}' ({ref})")

return sn_rates


def get_sn_fraction(sn_subclass=None, sn_subclass_fractions_name=None):
"""Return SN rates for specific types. These are taken from
https://arxiv.org/pdf/1509.06574.pdf, and are assumed to be fixed
fractions of the overall SN rate. Acceptable types are:
Expand All @@ -26,18 +50,24 @@ def get_sn_fraction(sn_subclass=None):
SNIbc (equal to Ib + Ic)
:param sn_subclass: Type of SN
:param sn_subclass_fractions_name: Name of estimates to use for relative rates of each subclass
:return: fraction represented by that subtype
"""
if sn_subclass is None:
logging.info("No specified subclass of supernova. Assuming overall CCSN rate.")
sn_subclass = "all"
return 1.0

if sn_subclass in sn_types.keys():
logging.info(f"Subclass '{sn_subclass}' is equal to {100.*sn_types[sn_subclass]:.2f}% of the CCSN rate.")
return sn_types[sn_subclass]
else:
raise Exception(f"Supernova type '{sn_subclass}' not recognised. "
f"The following types are available: {sn_types.keys()}")

sn_types = get_sn_subfraction(sn_subclass_fractions_name)

if sn_subclass in sn_types.keys():
logging.info(f"Subclass '{sn_subclass}' is equal to "
f"{100.*sn_types[sn_subclass]:.2f}% of the CCSN rate.")
return sn_types[sn_subclass]
else:
raise Exception(f"Supernova type '{sn_subclass}' not recognised. "
f"The following types are available: {sn_types.keys()}")

kcc_rates = {
"madau_14": (0.0068 / u.solMass, "http://arxiv.org/abs/1403.0007v3"),
Expand Down Expand Up @@ -65,13 +95,12 @@ def get_kcc_rate(kcc_name=None):
return kcc


def get_local_ccsn_rate(rate_name=None, kcc_name=None, sn_subclass=None, fraction=1.0):
def get_local_ccsn_rate(rate_name=None, kcc_name=None, sn_subclass=None):
"""Returns a local rate of core-collapse supernovae (CCSNe).
:param rate_name: Name of local Star Formation Rate (sfr) to be used
:param kcc_name: Name of kcc (sn per unit star formation) to be used
:param sn_subclass: Name of subclass of CCSNe to use
:param fraction: Fraction of specified rate to include
:return: Local rate
"""

Expand All @@ -85,7 +114,7 @@ def get_local_ccsn_rate(rate_name=None, kcc_name=None, sn_subclass=None, fractio

subclass_fraction = get_sn_fraction(sn_subclass)

return sfr_rate * kcc * subclass_fraction * fraction
return sfr_rate * kcc * subclass_fraction

def get_ccsn_rate(evolution_name=None, rate_name=None, kcc_name=None, sn_subclass=None, **kwargs):
"""Returns a local rate of core-collapse supernovae (CCSNe) as a function of redshift.
Expand All @@ -94,7 +123,6 @@ def get_ccsn_rate(evolution_name=None, rate_name=None, kcc_name=None, sn_subclas
:param rate_name: Name of local Star Formation Rate (sfr) to be used
:param kcc_name: Name of kcc (sn per unit star formation) to be used
:param sn_subclass: Name of subclass of CCSNe to use
:param fraction: Fraction of specified rate to include
:return: Rate as a function of redshift
"""

Expand Down
43 changes: 43 additions & 0 deletions flarestack/cosmo/rates/fbot_rates.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
import logging
import numpy as np
from astropy.cosmology import Planck15 as cosmo
from astropy import units as u
from flarestack.cosmo.rates.sfr_rates import get_sfr_evolution

local_fbot_rates = {
"ho_20_high": (7. * 10**-7. / (u.Mpc**3 * u.yr), "https://arxiv.org/abs/2003.01222"),
"ho_20_low": (4. * 10**-7. / (u.Mpc**3 * u.yr), "https://arxiv.org/abs/2003.01222"),
}

def get_local_fbot_rate(rate_name=None):
"""Returns a local rate of Fast Blue Optical Transients (FBOTs).
:param rate_name: Name of local FBOT rate to be used
:return: Local rate
"""

if rate_name is None:
logging.info("No rate specified. Assuming default rate.")
rate_name = "ho_20_high"

if rate_name not in local_fbot_rates.keys():
raise Exception(f"Rate name '{rate_name}' not recognised. "
f"The following source evolutions are available: {local_fbot_rates.keys()}")
else:
local_rate, ref = local_fbot_rates[rate_name]
logging.info(f"Loaded rate '{rate_name}' ({ref})")

return local_rate.to("Mpc-3 yr-1")

def get_fbot_rate(evolution_name=None, rate_name=None, **kwargs):
"""Returns a local rate of core-collapse supernovae (CCSNe) as a function of redshift.
:param evolution_name: Name of Star Formation evolution to use
:param rate_name: Name of local FBOT rate to be used
:return: Rate as a function of redshift
"""

normed_evolution = get_sfr_evolution(evolution_name, **kwargs)
local_rate = get_local_fbot_rate(rate_name)

return lambda z: normed_evolution(z) * local_rate
13 changes: 10 additions & 3 deletions tests/test_cosmo.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,8 @@
"""
import logging
import unittest
from flarestack.cosmo import get_diffuse_flux_at_1GeV, get_diffuse_flux_at_100TeV, calculate_transient_cosmology, get_rate
from flarestack.cosmo import get_diffuse_flux_at_1GeV, get_diffuse_flux_at_100TeV, \
get_diffuse_flux, calculate_transient_cosmology, get_rate
from flarestack.cosmo.icecube_diffuse_flux import contours, plot_diffuse_flux
from astropy import units as u

Expand Down Expand Up @@ -37,6 +38,9 @@ def test_get_diffuse_flux(self):

self.assertEqual(expected_1_gev, res_1_gev[0])

ratio = (res_1_gev[0]/get_diffuse_flux(1., fit=fit)[0]).to("").value
self.assertAlmostEqual(ratio, 1.0, places=3)

self.assertEqual(res_100_tev[1], res_1_gev[1])

fits = ["joint_15", "northern_tracks_17", "northern_tracks_19"]
Expand All @@ -46,6 +50,9 @@ def test_get_diffuse_flux(self):

self.assertEqual(res_100_tev, default_flux_100TeV[i])

ratio = (res_100_tev[0]/get_diffuse_flux(10.**5, fit=fit)[0]).to("").value
self.assertAlmostEqual(ratio, 1.0, places=3)

logging.info("Calculated values {0}".format(res_100_tev))
logging.info("Reference values {0}".format(default_flux_100TeV[i]))

Expand Down Expand Up @@ -75,8 +82,8 @@ def test_neutrino_cosmology(self):

def test_plotting(self):

for label, (_, _, contour_68, contour_95, e_range, _) in contours.items():
plot_diffuse_flux(label, contour_68, contour_95, e_range)
for label in contours.keys():
plot_diffuse_flux(label)

if __name__ == '__main__':
unittest.main()

0 comments on commit cb65458

Please sign in to comment.