Introduction to 
Links
Slides: https://events.icecube.wisc.edu/event/140/contributions/7689/attachments/6114/7347/IceCube_tech_workshop_flarestack.pdf
GitHub: https://github.com/icecube/flarestack
More examples: https://mybinder.org/v2/gh/icecube/flarestack/master
Logging
uses logging. So to get an idea of what is going on, let’s turn it on!
(We will use the INFO level. If you run into a problem and really want to know what’s going on use DEBUG)
import logging
logging.basicConfig(level='INFO')
1. Installation
Execute
pip install --prefix=/path/to/lib flarestack
to install into /path/to/lib space prior to executing this notebook.
This will also install all necessary dependencies.
For the purpose of this workshop: path/to/lib=/scratch/<your_username>/
Set up external directories
To specify your dataset directory do
export FLARESTACK_DATASET_DIR=/path/to/datasets
If you are running at DESY or WIPAC this is not necesarry.
The code products such as data, plots, caches, etc will be stored in a separate scratch directory.
To specify run
export FLARESTACK_SCRATCH_DIR=/path/to/scratch
If not given, defaults to home directory. For now we can use /scratch/<your_username>/flarestack_scratch
from flarestack.shared import host_server
from flarestack.data.icecube.ic_season import icecube_dataset_dir
print(f'Running at {host_server}, data directory is {icecube_dataset_dir}')
INFO:flarestack.shared:Scratch Directory is: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/storage/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/cluster/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/pull_corrections/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/cluster/logs/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/catalogues/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/acceptance_functions/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/energy_pdf_splines/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/storage/pickles/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/plots/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/SoB_splines/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/analysis/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/plots/illustrations/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/catalogues/transients/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/bkg_splines/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/dataset_plots/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/unblinding_results/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/limits/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/pull_corrections/pulls/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/pull_corrections/floors/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/storage/cache/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/storage/cache/catalogue_cache/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/public_datasets/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/energy_proxy_weighting/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/dataset_plots/effective_area_plots/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/median_angular_resolution/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/dataset_plots/angular_resolution_plots/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/dataset_plots/energy_proxy_map/
INFO:flarestack.shared:Found Directory: /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/sim_datasets/
INFO:numexpr.utils:Note: NumExpr detected 64 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
INFO:numexpr.utils:NumExpr defaulting to 8 threads.
INFO:flarestack.data.icecube.ic_season:Loading datasets from /lustre/fs22/group/icecube/data_mirror/ (DESY)
Running at DESY, data directory is /lustre/fs22/group/icecube/data_mirror/
from flarestack.shared import fs_scratch_dir
print(f'Scratch directory is {fs_scratch_dir}')
Scratch directory is /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/
2. Using Flarestack Classes
Classes used in ’s core functionality (e.g. flarestack.core.energy_pdf.EnergyPDF, flarestack.core.minimisation.MinimisationHandler, etc) have a class attribute <class>.subclasses.
This is a dictionary with the structure {<subclass name>: <subclass>}.
from flarestack.core.minimisation import MinimisationHandler
MinimisationHandler.subclasses
{'fixed_weights': flarestack.core.minimisation.FixedWeightMinimisationHandler,
'large_catalogue': flarestack.core.minimisation.LargeCatalogueMinimisationHandler,
'fit_weights': flarestack.core.minimisation.FitWeightMinimisationHandler,
'flare': flarestack.core.minimisation.FlareMinimisationHandler}
For analyses we only have to pass a dictionary of the subclass names and corresponding parameters.
To execute use flarestack.cluster.submitter.Submitter. This always works locally. For using the cluster, again, if you are running at DESY or WIPAC, you do not have to worry. We got you covered.
from flarestack.cluster.submitter import Submitter
Submitter.submitter_dict
{'local': flarestack.cluster.submitter.LocalSubmitter,
'DESY': flarestack.cluster.submitter.DESYSubmitter,
'WIPAC': flarestack.cluster.submitter.WIPACSubmitter}
3. Example: Point Source Sensitivity
Let’s try to calculate the 10-year point source sensitivity for one declination.
First we have to specify a name for the analysis.
name = "analyses/10yr_ps_sens_one_declination"
The input directory (with the analysis dictionaries), the output directory (plots, p-values, etc) and the cache directory (saved trials, etc) will be created accordingly.
For example our plot output directory will be:
from flarestack.shared import plot_output_dir
plot_output_dir(name)
'/afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/plots/analyses/10yr_ps_sens_one_declination'
Many dataset implementations are available in flarestack.data. We will use the PS Tracks v3.2
from flarestack.data.icecube import ps_v003_p02
We want to inject a steady neutrino signal with a power law spectrum with 
. For other Energy or Time PDFs check flarestack.core.energy_pdf and flarestack.core.time_pdf.
This is as straight forward as:
injection_energy = {
"energy_pdf_name": "power_law",
"gamma": 2.5
}
injection_time = {
"time_pdf_name": "steady"
}
inj_kwargs = {
"injection_energy_pdf": injection_energy,
"injection_sig_time_pdf": injection_time
}
We are looking for a steady signal with a power law spectrum.
We assume the background to be constant in time.
We want to use the “standard” point source likelihood. More likelihood implementations in flarestack.core.llh
llh_time = {
"time_pdf_name": "steady"
}
llh_energy = {
"energy_pdf_name": "power_law",
}
llh_time_bkg = {
"time_pdf_name": "steady"
}
llh_kwargs = {
"llh_name": "standard",
"llh_energy_pdf": llh_energy,
"llh_sig_time_pdf": llh_time,
"llh_bkg_time_pdf": llh_time_bkg
}
We need a source catalogue. This catalogue will be a numpy array stored as a .npy file and we only pass the filename.
For point sources the is a uitility function to generate dummy sources.
from flarestack.utils.prepare_catalogue import ps_catalogue_name
import numpy as np
sindec = 0.5
catalogue_path = ps_catalogue_name(sindec)
print(f'your catalogue is located at {catalogue_path}')
np.load(catalogue_path)
your catalogue is located at /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/input/catalogues/single_source/sindec_0.50.npy
array([(3.14159265, 0.52359878, 1., 1., 55800.4164699, 55750.4164699, 55900.4164699, 1., b'PS_dec=0.5')],
dtype=[('ra_rad', '<f8'), ('dec_rad', '<f8'), ('base_weight', '<f8'), ('injection_weight_modifier', '<f8'), ('ref_time_mjd', '<f8'), ('start_time_mjd', '<f8'), ('end_time_mjd', '<f8'), ('distance_mpc', '<f8'), ('source_name', 'S30')])
Now we make a guess for our sensitivity.
Note: is using its own scale factor 
from flarestack.shared import flux_to_k
flux_to_k(1)
999999999.9999999
Here we know where the sensitivity should be. Because the analysis has been done before.
from flarestack.icecube_utils.reference_sensitivity import reference_sensitivity
scale = flux_to_k(reference_sensitivity(sindec)) * 3
Now we just have to put all the info into one dictionary to pass to the MinimisationHanddler
mh_dict = {
"name": name, # unique name for the analysis
"mh_name": "fixed_weights", # name of the MinimisationHandler subcalss
"dataset": ps_v003_p02, # the neutrino dataset
"catalogue": catalogue_path, # path to the .npy catalogue file
"inj_dict": inj_kwargs, # info for the Injector
"llh_dict": llh_kwargs, # info for the LLH
"scale": scale, # a guess for the sensitivity scale
"n_trials": 10, # number of trials to run (background trials will be run ten times this number!)
"n_steps": 10, # number of steps when injecting signal
"allow_extrapolated_sensitivity": True # allow extrapolation in the sensitivity calculation (here we do because we only run very few trials)
}
To execute the analysis we defined above we create a submitter instance
submitter = Submitter.get_submitter(
mh_dict=mh_dict, # the analysis info
use_cluster=False, # run it on the cluster if True
n_cpu=1, # number of LOCAL CPUs to use, NOTE: the number of cluster CPUs has to be specified in the cluster_kwargs!
do_sensitivity_scale_estimation=False, # make a guess of the sensitivity scale, for options check flarestack.cluster.submitter
remove_old_results=True, # if you are running the analysis again and something changed, maybe you want to remove old trials?
# **cluster_kwargs # keyword arguments used when running the cluster, This depends on the cluster obviously
)
print(submitter)
----- Submitter for analyses/10yr_ps_sens_one_declination -----
not using cluster
using 1 CPUs locally
job-id: None
no scale estimation
Energise …
submitter.analyse()
INFO:flarestack.core.minimisation:Using 'standard' LLH class
INFO:flarestack.core.injector:Initialising Injector for IC40
WARNING:flarestack.icecube_utils.dataset_loader:No field called 'good_i3' found in GoodRunList. Cannot check if all runs in GoodRunList are actually good.
INFO:flarestack.core.injector:Initialising Injector for IC59
WARNING:flarestack.icecube_utils.dataset_loader:No field called 'good_i3' found in GoodRunList. Cannot check if all runs in GoodRunList are actually good.
INFO:flarestack.core.injector:Initialising Injector for IC79
WARNING:flarestack.icecube_utils.dataset_loader:No field called 'good_i3' found in GoodRunList. Cannot check if all runs in GoodRunList are actually good.
INFO:flarestack.core.injector:Initialising Injector for IC86_2011
WARNING:flarestack.icecube_utils.dataset_loader:No field called 'good_i3' found in GoodRunList. Cannot check if all runs in GoodRunList are actually good.
INFO:flarestack.core.injector:Initialising Injector for IC86_2012_17
WARNING:flarestack.icecube_utils.dataset_loader:No field called 'good_i3' found in GoodRunList. Cannot check if all runs in GoodRunList are actually good.
ERROR:flarestack.data.icecube.ic_season:The IceCube GoodRunList was not produced correctly.
ERROR:flarestack.data.icecube.ic_season:Some runs in the GoodRunList start immediately after the preceding run ends.
ERROR:flarestack.data.icecube.ic_season:There should be gaps between every run due to detector downtime, but some are missing here.
ERROR:flarestack.data.icecube.ic_season:The first missing gap is between runs 120700 and 120701.
ERROR:flarestack.data.icecube.ic_season:Any livetime estimates using this GoodRunList will not be accurate.
ERROR:flarestack.data.icecube.ic_season:This is a known problem affecting older IceCube GoodRunLists.
ERROR:flarestack.data.icecube.ic_season:You should use a newer, corrected GoodRunList.
ERROR:flarestack.data.icecube.ic_season:Flarestack will attempt to stitch these runs together.
ERROR:flarestack.data.icecube.ic_season:However, livetime estimates may be off by several percentage points, or even more for very short timescales.
ERROR:flarestack.data.icecube.ic_season:You have been warned!
ERROR:flarestack.data.icecube.ic_season:The IceCube GoodRunList was not produced correctly.
ERROR:flarestack.data.icecube.ic_season:Some runs in the GoodRunList start before the preceding run has ended.
ERROR:flarestack.data.icecube.ic_season:Under no circumstances should runs overlap.
ERROR:flarestack.data.icecube.ic_season:The first overlap is between runs 126953 and 126954.
ERROR:flarestack.data.icecube.ic_season:Any livetime estimates using this GoodRunList will not be accurate.
ERROR:flarestack.data.icecube.ic_season:This is a known problem affecting older IceCube GoodRunLists.
ERROR:flarestack.data.icecube.ic_season:You should use a newer, corrected GoodRunList.
ERROR:flarestack.data.icecube.ic_season:Flarestack will attempt to stitch these runs together.
ERROR:flarestack.data.icecube.ic_season:However, livetime estimates may be off by several percentage points, or even more for very short timescales.
ERROR:flarestack.data.icecube.ic_season:You have been warned!
INFO:flarestack.core.multiprocess_wrapper:Added 190 trials to queue. Now processing.
INFO:flarestack.core.multiprocess_wrapper:190 tasks remaining.
INFO:flarestack.core.minimisation:Using 'standard' LLH class
INFO:flarestack.core.multiprocess_wrapper:146 tasks remaining.
INFO:flarestack.core.multiprocess_wrapper:93 tasks remaining.
INFO:flarestack.core.multiprocess_wrapper:43 tasks remaining.
INFO:flarestack.core.multiprocess_wrapper:Finished processing 190 tasks.
To get the results we use the ResultsHandler. This will also create some plots like the sensitivity fit, bias plots, etc. in the plot directory.
from flarestack.core.results import ResultsHandler
results_handler = ResultsHandler(submitter.mh_dict)
INFO:flarestack.core.results:Saving bias plot to /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/plots/analyses/10yr_ps_sens_one_declination/bias_n_s.pdf
INFO:flarestack.core.results:Saving bias plot to /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/output/plots/analyses/10yr_ps_sens_one_declination/bias_gamma.pdf
INFO:flarestack.core.results:Fraction of overfluctuations is 0.50 above 0.04 (N_trials=100) (Scale=0)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.20 above 0.04 (N_trials=10) (Scale=0.1924)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.60 above 0.04 (N_trials=10) (Scale=0.3847)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.40 above 0.04 (N_trials=10) (Scale=0.5771)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.80 above 0.04 (N_trials=10) (Scale=0.7695)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.70 above 0.04 (N_trials=10) (Scale=0.9618)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.50 above 0.04 (N_trials=10) (Scale=1.154)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.60 above 0.04 (N_trials=10) (Scale=1.347)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.30 above 0.04 (N_trials=10) (Scale=1.539)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.50 above 0.04 (N_trials=10) (Scale=1.731)
WARNING:flarestack.core.results:The sensitivity is beyond the range of the tested scales.The number is probably not good.
INFO:flarestack.core.results:EXTRAPOLATED Sensitivity is 4.22e-09
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=100) (Scale=0)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=100) (Scale=0)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=0.1924)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=0.1924)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=0.3847)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=0.3847)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=0.5771)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=0.5771)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=0.7695)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=0.7695)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=0.9618)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=0.9618)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=1.154)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=1.154)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=1.347)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=1.347)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=1.539)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=1.539)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 26.63 (N_trials=10) (Scale=1.731)
INFO:flarestack.core.results:Fraction of overfluctuations is 0.00 above 25 (N_trials=10) (Scale=1.731)
WARNING:flarestack.core.results:RuntimeError for discovery potential:
Optimal parameters not found: Number of calls to function has reached maxfev = 800.
print(fr'sensitivity flux: {results_handler.sensitivity:.2e} +{results_handler.sensitivity_err[1]} -{results_handler.sensitivity_err[0]}')
print(f'reference: {reference_sensitivity(sindec)[0]}')
print(fr'sensitivity n_s: {results_handler.sensitivity * results_handler.flux_to_ns:.2e} +{results_handler.sensitivity_err[1] * results_handler.flux_to_ns} -{results_handler.sensitivity_err[0] * results_handler.flux_to_ns}')
sensitivity flux: 4.22e-09 +3.73093684037954e-09 -1.348702961817186e-09
reference: 5.770894709671131e-10
sensitivity n_s: 5.39e-01 +0.4759829693036423 -0.1720639260162485
4. Example: Upper Limits
There are some really useful functions in flarestack.cosmo!
Take the implementations of the IceCube diffuse flux measurements for example:
from flarestack.cosmo.icecube_diffuse_flux import contours, get_diffuse_flux_contour
contours.keys()
dict_keys(['joint_15', 'northern_tracks_16', 'northern_tracks_17', 'northern_tracks_19'])
import matplotlib.pyplot as plt
plt.figure()
ax = plt.subplot(111)
for fit in contours.keys():
best_fit, upper_butterfly, lower_butterfly, e_range = get_diffuse_flux_contour(fit)
plt.plot(e_range, best_fit(e_range) * e_range**2, label=fit)
plt.fill_between(e_range, upper_butterfly(e_range)* e_range**2, lower_butterfly(e_range)* e_range**2, alpha=0.3)
plt.yscale("log")
plt.xscale("log")
plt.xlabel(r"$E_{\nu}$")
plt.ylabel(r"$E_{\nu}^{2} \frac{dN}{dE}$")
ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.4), ncol=2, fancybox=True, shadow=True, fontsize=12)
plt.tight_layout()
plt.show()
plt.close()
INFO:flarestack.cosmo.icecube_diffuse_flux:Loaded contour 'joint_15' from https://arxiv.org/abs/1507.03991
INFO:flarestack.cosmo.icecube_diffuse_flux:Loaded contour 'northern_tracks_16' from https://arxiv.org/abs/1607.08006
INFO:flarestack.cosmo.icecube_diffuse_flux:Loaded contour 'northern_tracks_17' from https://doi.org/10.22323/1.301.1005
INFO:flarestack.cosmo.icecube_diffuse_flux:Loaded contour 'northern_tracks_19' from https://arxiv.org/abs/1908.09551
If you are interested in transients you also are at the right place!
With your favourte transients population’s rate, the flux normaisation at 1 GeV and the corresponding spectral index you can easily get the rate in a redshift shell, the neutrino flux per source at a certain redshift, the neutrino flux per redshift and the cumulatice neutrino flux.
from flarestack.cosmo.neutrino_cosmology import define_cosmology_functions
from flarestack.cosmo.rates import get_rate
from astropy import units as u
frb_rate = get_rate('FRB')
frb_dummy_flux = 2e+46 / u.GeV
frb_example_gamma = 2
rate_per_z, nu_flux_per_z, nu_flux_per_source, cumulative_nu_flux = define_cosmology_functions(frb_rate, frb_dummy_flux, frb_example_gamma)
redshift = np.linspace(0.1, 4, 1000)
fig, axs = plt.subplots(4, sharex='all', figsize=[5, 12])
axs[0].plot(redshift, rate_per_z(redshift).to('yr-1').value)
axs[0].set_ylabel('rate [yr$^{-1}$]')
axs[1].plot(redshift, nu_flux_per_z(redshift).to('GeV-1 cm-2 s-1 sr-1').value)
axs[1].set_ylabel(r'$\nu$ flux per redshift [GeV$^{-1}$ cm$^{-2}$]')
axs[2].plot(redshift, nu_flux_per_source(redshift).to('GeV-1 cm-2').value)
axs[2].set_ylabel(r'$\nu$ flux per source [GeV$^{-1}$ cm$^{-2}$]')
axs[3].plot(redshift[1:-1], [i.to('1 / (cm2 GeV s sr)').value for i in cumulative_nu_flux(redshift)])
axs[3].set_ylabel(r'cumulative $\nu$ flux [GeV$^{-1}$ cm$^{-2}$]')
plt.show()
plt.close()
INFO:flarestack.cosmo.rates:Loading source class 'frb'
INFO:root:No evolution specified. Assuming default evolution.
INFO:root:Loaded evolution 'madau_14' (http://arxiv.org/abs/1403.0007v3)
INFO:root:No rate specified. Assuming default rate.
INFO:root:Loaded rate 'bochenek_20' (https://arxiv.org/abs/2005.10828)
NOTE: The result of cumulative_nu_flux() is already your result for the contribution of the popultion to the diffuse flux!
All the above is packed into one convenience function:
from flarestack.cosmo.neutrino_cosmology import calculate_transient_cosmology
Let’s use this to get some actual super interesting, timely results and revisit the FRB asscociated with the galactic Magnetar SGR 1935+2154 (https://arxiv.org/abs/2005.10828).
IceCube performed a search for neutrinos and found upper limits:
IceCube Limit is E^2 dN/dE = 5.2 × 10−2 GeV cm^-2 @ 1 GeV
(http://www.astronomerstelegram.org/?read=13689)
The Magnetar is 16 kpc away (conservative). Let’s assume a spectrum with 
from flarestack.core.energy_pdf import EnergyPDF
dist = 16 * u.kpc
atel_flux_norm_lim = 5.2 * 10**-2. * (u. GeV / u.cm**2) / u.GeV**2.
e_pdf_dict = {
"energy_pdf_name": "power_law",
"gamma": 2.0,
"e_min_gev": 10.**3,
"e_max_gev": 10.**6,
"nu_flux_at_1_gev": atel_flux_norm_lim * 4 * np.pi * dist**2.
}
epdf = EnergyPDF.create(e_pdf_dict)
INFO:flarestack.core.energy_pdf:Minimum Energy is 1e+03 GeV.
INFO:flarestack.core.energy_pdf:Maximum Energy is 1e+06 GeV.
With these information it is now super straight forward to get the upper limits from a population of FRB’s, that share SGR 1935+2154’s properties, assuming they are all standard candels:
fit = "joint_15"
integrated_nu_flux_1_gev = calculate_transient_cosmology(e_pdf_dict, frb_rate, "frb_limit", zmax=8.0, diffuse_fit=fit)
INFO:flarestack.cosmo.icecube_diffuse_flux:Loaded contour 'joint_15' from https://arxiv.org/abs/1507.03991
INFO:flarestack.cosmo.neutrino_cosmology:Using the joint_15 best fit values of the diffuse flux.
INFO:flarestack.cosmo.neutrino_cosmology:Diffuse Flux at 1 GeV: 7.0624201077093805e-06 1 / (cm2 GeV s sr)
INFO:flarestack.cosmo.neutrino_cosmology:Diffuse Spectral Index is 2.5
INFO:flarestack.cosmo.neutrino_cosmology:Neutrino Flux at 1 GeV is 1.6e+45 1 / GeV
INFO:flarestack.cosmo.neutrino_cosmology:Local rate is 0.072 1 / (Mpc3 yr)
INFO:flarestack.cosmo.neutrino_cosmology:Cumulative sources at z=8.0: 2.9e+11
INFO:flarestack.cosmo.neutrino_cosmology:Cumulative flux at z=8.0 (1 GeV): 2.4e-10 1 / (cm2 GeV s sr)
INFO:flarestack.cosmo.neutrino_cosmology:Cumulative annual flux at z=8.0 (1 GeV): 0.0075 1 / (cm2 GeV sr)
INFO:flarestack.cosmo.neutrino_cosmology:Fraction of diffuse flux at 1GeV: 3.4e-05
INFO:flarestack.cosmo.neutrino_cosmology:Cumulative neutrino flux 2.4e-10 1 / (cm2 GeV s sr)
INFO:flarestack.cosmo.neutrino_cosmology:Fraction of flux from nearby (z<0.1) sources: 0.057
INFO:flarestack.cosmo.neutrino_cosmology:Fraction of flux from nearby (z<0.3) sources: 0.17
Let’s plot our very interesting findings so we can publish in a prestegious journal!
best_fit, upper_butterfly, lower_butterfly, e_range = get_diffuse_flux_contour(fit=fit)
plt.figure()
plt.plot(e_range, best_fit(e_range) * e_range**2, label="IceCube Diffuse Flux")
plt.fill_between(e_range, upper_butterfly(e_range)* e_range**2, lower_butterfly(e_range)* e_range**2, alpha=0.3)
x = [epdf.e_min, np.exp(0.5*(np.log(epdf.e_min) + np.log(epdf.e_max))), epdf.e_max]
y = np.array([integrated_nu_flux_1_gev.value for _ in range(3)])
plt.errorbar(x, y, yerr=0.25*y, uplims=True, label='limit')
plt.yscale("log")
plt.xscale("log")
plt.xlabel(r"$E_{\nu} [GeV] $")
plt.ylabel(r"$E_{\nu}^{2} \frac{dN}{dE}$ [GeV cm$^{-2}$ s$^{-1}$ sr$^{-1}$]")
plt.legend()
plt.show()
plt.close()
INFO:flarestack.cosmo.icecube_diffuse_flux:Loaded contour 'joint_15' from https://arxiv.org/abs/1507.03991
