<jupyter_start><jupyter_text>Introduction to $\texttt{flarestack}$ 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 $\texttt{flarestack}$ 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)<jupyter_code>import logging
logging.basicConfig(level='INFO')<jupyter_output><empty_output><jupyter_text>1. Installation Execute```bashpip install --prefix=/path/to/lib flarestack```to install $\texttt{flarestack}$ 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//` Set up external directories To specify your **dataset directory** do```bashexport 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```bashexport FLARESTACK_SCRATCH_DIR=/path/to/scratch```If not given, defaults to home directory. For now we can use `/scratch//flarestack_scratch`<jupyter_code>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}')
from flarestack.shared import fs_scratch_dir
print(f'Scratch directory is {fs_scratch_dir}')<jupyter_output>Scratch directory is /afs/ifh.de/user/n/neckerja/scratch/fs_dev_branch/flarestack__data/<jupyter_text>2. Using Flarestack ClassesClasses used in $\texttt{flarestack}$'s core functionality (e.g. `flarestack.core.energy_pdf.EnergyPDF`, `flarestack.core.minimisation.MinimisationHandler`, etc) have a class attribute `.subclasses`.  This is a dictionary with the structure `{: }`.<jupyter_code>from flarestack.core.minimisation import MinimisationHandler
MinimisationHandler.subclasses<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>from flarestack.cluster.submitter import Submitter
Submitter.submitter_dict<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>name = "analyses/10yr_ps_sens_one_declination"<jupyter_output><empty_output><jupyter_text>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:<jupyter_code>from flarestack.shared import plot_output_dir
plot_output_dir(name)<jupyter_output><empty_output><jupyter_text>Many dataset implementations are available in `flarestack.data`. We will use the PS Tracks v3.2<jupyter_code>from flarestack.data.icecube import ps_v003_p02<jupyter_output><empty_output><jupyter_text>We want to inject a steady neutrino signal with a power law spectrum with $\gamma=2.5$. For other Energy or Time PDFs check `flarestack.core.energy_pdf` and `flarestack.core.time_pdf`.   \This is as straight forward as:<jupyter_code>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
}<jupyter_output><empty_output><jupyter_text>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`<jupyter_code>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
}<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>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)<jupyter_output>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<jupyter_text>Now we make a guess for our sensitivity.   Note: $\texttt{flarestack}$ is using its own scale factor $k$<jupyter_code>from flarestack.shared import flux_to_k
flux_to_k(1)<jupyter_output><empty_output><jupyter_text>Here we know where the sensitivity should be. Because the analysis has been done before.<jupyter_code>from flarestack.icecube_utils.reference_sensitivity import reference_sensitivity
scale = flux_to_k(reference_sensitivity(sindec)) * 3<jupyter_output><empty_output><jupyter_text>Now we just have to put all the info into one dictionary to pass to the `MinimisationHanddler`<jupyter_code>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)
}<jupyter_output><empty_output><jupyter_text>To execute the analysis we defined above we create a submitter instance<jupyter_code>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)<jupyter_output>----- Submitter for analyses/10yr_ps_sens_one_declination -----
not using cluster 
using 1 CPUs locally
job-id: None 
no scale estimation<jupyter_text>Energise ......<jupyter_code>submitter.analyse()<jupyter_output>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
W[...]<jupyter_text>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.<jupyter_code>from flarestack.core.results import ResultsHandler
results_handler = ResultsHandler(submitter.mh_dict)
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}')<jupyter_output>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<jupyter_text>4. Example: Upper LimitsThere are some really useful functions in `flarestack.cosmo`!  Take the implementations of the IceCube diffuse flux measurements for example:<jupyter_code>from flarestack.cosmo.icecube_diffuse_flux import contours, get_diffuse_flux_contour
contours.keys()
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()<jupyter_output>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<jupyter_text>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.<jupyter_code>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()<jupyter_output>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)<jupyter_text>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:<jupyter_code>from flarestack.cosmo.neutrino_cosmology import calculate_transient_cosmology<jupyter_output><empty_output><jupyter_text>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 $\gamma=2$<jupyter_code>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)<jupyter_output>INFO:flarestack.core.energy_pdf:Minimum Energy is 1e+03 GeV.
INFO:flarestack.core.energy_pdf:Maximum Energy is 1e+06 GeV.<jupyter_text>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:<jupyter_code>fit = "joint_15"
integrated_nu_flux_1_gev = calculate_transient_cosmology(e_pdf_dict, frb_rate, "frb_limit", zmax=8.0, diffuse_fit=fit)<jupyter_output>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)[...]<jupyter_text>Let's plot our very interesting findings so we can publish in a prestegious journal!<jupyter_code>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()<jupyter_output>INFO:flarestack.cosmo.icecube_diffuse_flux:Loaded contour 'joint_15' from https://arxiv.org/abs/1507.03991
{'fraction_non_alphanumeric': 0.10454406831351022, 'fraction_numerical': 0.030314120158584933, 'mean_word_length': 5.417221135029354, 'pattern_counts': {'":': 27, '<': 79, '<?xml version=': 0, '>': 79, 'https://': 11, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '20859', 'license': 'MIT', 'n_tokens_mistral': 6062, 'n_tokens_neox': 5542, 'n_words': 1431, 'path': 'examples/ipython_notebooks/flarestack_llh_workshop.ipynb', 'repo_name': 'robertdstein/flarestack', 'size': '16395'}