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
robertdstein committed May 10, 2020
2 parents dd62ad8 + 80e43a0 commit ef573f8
Show file tree
Hide file tree
Showing 19 changed files with 1,224 additions and 70 deletions.
348 changes: 348 additions & 0 deletions examples/ipython_notebooks/analysis_with_simulated_data.ipynb
View file

Large diffs are not rendered by default.

415 changes: 415 additions & 0 deletions examples/ipython_notebooks/simulating_a_dataset.ipynb
View file

Large diffs are not rendered by default.

64 changes: 64 additions & 0 deletions flarestack/analyses/general/check_stack_bias.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
import numpy as np
from astropy.table import Table
import argparse
from flarestack.core.results import ResultsHandler
from flarestack.data.icecube import ps_v002_p01
from flarestack.shared import plot_output_dir, flux_to_k
from flarestack.analyses.ccsn.stasik_2017.ccsn_limits import limits, get_figure_limits, p_vals
from flarestack.analyses.ccsn.necker_2019.ccsn_helpers import sn_cats, updated_sn_catalogue_name, \
sn_time_pdfs, raw_output_dir, pdf_names, limit_sens
from flarestack.analyses.ccsn import get_sn_color
from flarestack.cluster import analyse
from flarestack.cluster.run_desy_cluster import wait_for_cluster
import math
import pickle
from flarestack.utils.custom_dataset import custom_dataset
import os
import logging
import time
logging.getLogger().setLevel("INFO")
injection_energy = {
"energy_pdf_name": "power_law",
"gamma": 2.0
}

injection_time = {
'time_pdf_name': 'steady'
}

llh_energy = {
"energy_pdf_name": "power_law"
}

llh_time = {
'time_pdf_name': 'steady'
}

inj_dict = {
'injection_energy_pdf': injection_energy,
'injection_sig_time_pdf': injection_time
}

llh_dict = {
"llh_name": "standard",
"llh_energy_pdf": llh_energy,
"llh_sig_time_pdf": llh_time,
"llh_bkg_time_pdf": {"time_pdf_name": "steady"}
}

mh_dict = {
"name": "examples/crosscheck_stacking",
"mh_name": 'fixed_weights',
"dataset": ps_v002_p01.get_seasons(),
"catalogue": updated_sn_catalogue_name("IIn"),
# "catalogue": ps_stack_catalogue_name(0.1, 0.3),
# "catalogue": tde_catalogue_name("jetted"),
"inj_dict": inj_dict,
"llh_dict": llh_dict,
"scale": 10.,
"n_trials": 50,
"n_steps": 5
}

analyse(mh_dict, cluster=False)
rh = ResultsHandler(mh_dict)
32 changes: 25 additions & 7 deletions flarestack/core/energy_pdf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,25 +104,37 @@ def f(energy):

def integrate_over_E(self, f, lower=None, upper=None):
"""Uses Newton's method to integrate function f over the energy
range. By default, uses 100GeV to 10PeV, unless otherwise specified.
range. By default, uses 100 GeV to 10 PeV, unless otherwise specified.
Uses 1000 logarithmically-spaced bins to calculate integral.
:param f: Function to be integrated
:param lower: Lower bound for integration
:param upper: Upper bound for integration
:return: Integral of function
"""

diff_sum, _ = self.piecewise_integrate_over_energy(f, lower, upper)
if lower is None:
lower = self.integral_e_min

int_sum = np.sum(diff_sum)
if upper is None:
upper = self.integral_e_max

return self.integrate(f, lower, upper)

@staticmethod
def integrate(f, lower, upper):
diff_sum, _ = EnergyPDF.piecewise_integrate(f, lower, upper)
int_sum = np.sum(diff_sum)
return int_sum

def piecewise_integrate_over_energy(self, f, lower=None, upper=None):
"""Uses Newton's method to integrate function f over the energy
range. By default, uses 100GeV to 10PeV, unless otherwise specified.
range. By default, uses 100 GeV to 10 PeV, unless otherwise specified.
Uses 1000 logarithmically-spaced bins to calculate integral.
:param f: Function to be integrated
:param lower: Lower bound for integration
:param upper: Upper bound for integration
:return: Integral of function bins
"""

Expand All @@ -132,14 +144,20 @@ def piecewise_integrate_over_energy(self, f, lower=None, upper=None):
if upper is None:
upper = self.integral_e_max

diff_sum, e_range = self.piecewise_integrate(f, lower, upper)

return diff_sum, e_range

@staticmethod
def piecewise_integrate(f, lower, upper):
nsteps = int(1.e3)

e_range = np.linspace(np.log10(lower), np.log10(upper), nsteps + 1)
diff_sum = []

for i, log_e in enumerate(e_range[:-1]):
e0 = np.exp(log_e)
e1 = np.exp(e_range[i + 1])
e0 = 10.**(log_e)
e1 = 10.**(e_range[i + 1])
diff_sum.append(0.5 * (e1 - e0) * (f(e0) + f(e1)))

return diff_sum, e_range
Expand Down Expand Up @@ -311,7 +329,7 @@ def fluence_integral(self, e_min=None, e_max=None):


def return_energy_parameters(self):
default = [2.]
default = [2.0]
bounds = [(gamma_range[0], gamma_range[1])]
name = ["gamma"]
return default, bounds, name
Expand Down
8 changes: 4 additions & 4 deletions flarestack/core/injector.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -458,7 +458,7 @@ def make_injection_band_mask(self):

del injection_band_mask

logging.info("Saving to {0}".format(path))
logging.info(f"Saving to {path}")

def load_band_mask(self, index):
path = self.injection_band_paths[index]
Expand Down Expand Up @@ -552,7 +552,7 @@ def inject_signal(self, scale):

random_fraction = [random.random() for _ in range(n_s)]

sim_ev["logE"] = np.log10(np.exp(convert_f(random_fraction)))
sim_ev["logE"] = convert_f(random_fraction)

# Simulates times according to Time PDF

Expand Down Expand Up @@ -618,13 +618,13 @@ def source_eff_area(log_e):
self.energy_pdf.f(log_e) * \
self.energy_proxy_mapping(log_e)

start_x = np.log(self.energy_pdf.integral_e_min)
start_x = np.log10(self.energy_pdf.integral_e_min)



x_vals = np.linspace(
start_x + 1e-7,
np.log(self.energy_pdf.integral_e_max),
np.log10(self.energy_pdf.integral_e_max),
100
)[1:]

Expand Down
16 changes: 13 additions & 3 deletions flarestack/core/time_pdf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -229,11 +229,12 @@ def effective_injection_time(self, source=None):
raise NotImplementedError

def get_livetime(self):
if self.livetime is not None:
return self.livetime
raise NotImplementedError

def get_mjd_conversion(self):
raise NotImplementedError

return self.mjd_to_livetime, self.livetime_to_mjd

@TimePDF.register_subclass('steady')
class Steady(TimePDF):
Expand Down Expand Up @@ -352,6 +353,16 @@ def __init__(self, t_pdf_dict, season):
self.pre_window -= self.offset
self.post_window += self.offset

try:
if self.t_dict["livetime"] is True:
logging.debug("Using time PDF as a detector livetime PDF.")
self.mjd_to_livetime = lambda x: x - self.sig_t0([])
self.livetime_to_mjd = lambda x: x + self.sig_t0([])
self.livetime = self.t_dict["post_window"] + self.t_dict["pre_window"]
# print(self.t1, self.t0, self.livetime, self.livetime_to_mjd(7.), self.pre_window, self.sig_t0([]))
except KeyError:
pass

def sig_t0(self, source):
"""Calculates the starting time for the window, equal to the
source reference time in MJD minus the length of the pre-reference-time
Expand Down Expand Up @@ -434,7 +445,6 @@ def effective_injection_time(self, source=None):
t0 = self.mjd_to_livetime(self.sig_t0(source))
t1 = self.mjd_to_livetime(self.sig_t1(source))
time = (t1 - t0) * 60 * 60 * 24

return max(time, 0.)

def raw_injection_time(self, source):
Expand Down
15 changes: 9 additions & 6 deletions flarestack/data/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import os
import numpy as np
import copy
import logging
from numpy.lib.recfunctions import append_fields, drop_fields
from flarestack.core.injector import MCInjector, EffectiveAreaInjector
from flarestack.utils.make_SoB_splines import make_background_spline
Expand Down Expand Up @@ -33,9 +34,9 @@ def get_current(self):
if self.current is not None:
return self.datasets[self.current]
else:
print("Warning: no file listed as current.")
logging.warning("Warning: no file listed as current.")
key = sorted(list(self.datasets.keys()))
print("Using key {0} out of {1}".format(
logging.warning("Using key {0} out of {1}".format(
key, self.datasets.keys())
)
return self.datasets[key]
Expand All @@ -55,9 +56,9 @@ def add_subseason(self, season):
def get_seasons(self, *args, **kwargs):
season_names = list(args)
if len(season_names) == 0:
return copy.copy(self)
return self.make_copy()
else:
cd = copy.copy(self)
cd = self.make_copy()
cd.seasons = dict()
cd.subseasons = dict()
for name in season_names:
Expand All @@ -82,6 +83,8 @@ def get_seasons(self, *args, **kwargs):

for season in cd.values():
season.set_subselection_fraction(subselection_fraction)


return cd

def get_single_season(self, name):
Expand All @@ -108,6 +111,8 @@ def get(self, item):
def __len__(self):
return self.seasons.__len__()

def make_copy(self):
return copy.copy(self)


class Season:
Expand Down Expand Up @@ -280,8 +285,6 @@ class SeasonWithoutMC(Season):
def __init__(self, season_name, sample_name, exp_path, pseudo_mc_path,
**kwargs):



Season.__init__(self, season_name, sample_name, exp_path, **kwargs)
self.pseudo_mc_path = pseudo_mc_path
self.all_paths.append(self.pseudo_mc_path)
Expand Down
3 changes: 0 additions & 3 deletions flarestack/data/icecube/ps_tracks/ps_v003_p01.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@

# Add in old seasons before full detector operation, and IC86_1


def old_ic_season(season):
return IceCubeSeason(
season_name=season,
Expand All @@ -26,13 +25,11 @@ def old_ic_season(season):
log_e_bins=get_ps_binning(season)[1]
)


old_seasons = ["IC40", "IC59", "IC79", "IC86_2011"]

for season in old_seasons:
ps_v003_p01.add_season(old_ic_season(season))


# Add in combined IC86 2012-2017 seasons

new_years = ["2012", "2013", "2014", "2015", "2016", "2017"]
Expand Down
6 changes: 6 additions & 0 deletions flarestack/data/public/icecube/all_sky_point_source/README.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,3 +3,9 @@ These public IceCube datasets can be found in the IceCube website at:
https://icecube.wisc.edu/science/data/PS-3years

They can be used for point source analysis.

The mapping between energy proxy and true energy was also released by IceCube at:

https://icecube.wisc.edu/science/data/TXS0506_point_source

These were released for the analysis of TXS 0506+056, but are roughly right for other cases.
3 changes: 2 additions & 1 deletion flarestack/data/public/icecube/all_sky_point_source/all_sky_3_year.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@

zip_file = src_dir + "raw_data/3year-data-release.zip"

proxy_map_file = src_dir + "raw_data/Fig_S4_tabulated.txt"

output_base_dir = public_dataset_dir + "all_sky_3_year/"
extract_dir = output_base_dir + "extracted_data"
data_dir = extract_dir + "/3year-data-release/"
Expand Down Expand Up @@ -223,7 +225,6 @@ def run_all():

run_all()


if __name__ == "__main__":
run_all()
# mc = data_loader(ps_v002_p01[0]["mc_path"])
Expand Down
32 changes: 32 additions & 0 deletions flarestack/data/public/icecube/all_sky_point_source/raw_data/Fig_S4_tabulated.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
# Each row of the table below corresponds to a row of Fig. S4.
# The first two entries for each row give the neutrino energy range in log_10(E_nu/TeV).
# Each of the remaining 22 columns of the table corresponds to a column of Fig. S4.
# The columns span the muon energy proxy range (in log_10(E_mu/TeV)) from -0.4 to 4.0 in steps of 0.2
-1.00 , -0.75 : 3.03e-01 5.04e-01 1.78e-01 1.54e-02 1.38e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
-0.75 , -0.50 : 2.94e-01 4.94e-01 1.88e-01 2.18e-02 1.10e-03 1.96e-04 0.00e+00 3.09e-05 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
-0.50 , -0.25 : 2.45e-01 4.93e-01 2.22e-01 3.60e-02 3.95e-03 3.39e-04 2.31e-05 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
-0.25 , 0.00 : 1.95e-01 4.60e-01 2.68e-01 6.48e-02 1.00e-02 1.18e-03 7.95e-05 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
0.00 , 0.25 : 1.46e-01 4.02e-01 3.13e-01 1.10e-01 2.38e-02 4.68e-03 8.07e-04 1.09e-05 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
0.25 , 0.50 : 1.04e-01 3.24e-01 3.35e-01 1.68e-01 5.46e-02 1.24e-02 2.14e-03 2.67e-04 0.00e+00 2.13e-05 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
0.50 , 0.75 : 7.18e-02 2.46e-01 3.16e-01 2.28e-01 9.90e-02 2.92e-02 7.75e-03 1.56e-03 3.31e-04 1.01e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
0.75 , 1.00 : 4.95e-02 1.83e-01 2.62e-01 2.52e-01 1.59e-01 6.89e-02 1.89e-02 4.67e-03 8.87e-04 1.18e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
1.00 , 1.25 : 3.84e-02 1.34e-01 2.00e-01 2.25e-01 2.01e-01 1.32e-01 5.19e-02 1.58e-02 2.24e-03 8.43e-04 5.78e-05 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
1.25 , 1.50 : 2.93e-02 9.66e-02 1.58e-01 1.85e-01 1.94e-01 1.64e-01 1.09e-01 4.67e-02 1.38e-02 2.55e-03 3.87e-04 1.87e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
1.50 , 1.75 : 2.38e-02 8.12e-02 1.24e-01 1.50e-01 1.55e-01 1.57e-01 1.44e-01 1.06e-01 4.48e-02 1.16e-02 3.13e-03 2.83e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
1.75 , 2.00 : 1.83e-02 6.61e-02 1.04e-01 1.19e-01 1.31e-01 1.36e-01 1.40e-01 1.26e-01 1.01e-01 4.39e-02 1.17e-02 2.38e-03 3.13e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
2.00 , 2.25 : 1.64e-02 5.74e-02 8.77e-02 9.89e-02 1.13e-01 1.14e-01 1.18e-01 1.20e-01 1.19e-01 9.58e-02 4.59e-02 1.17e-02 2.05e-03 3.99e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
2.25 , 2.50 : 1.33e-02 4.19e-02 7.28e-02 8.76e-02 9.75e-02 1.02e-01 1.07e-01 1.05e-01 1.10e-01 1.02e-01 9.84e-02 4.47e-02 1.38e-02 2.41e-03 3.76e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
2.50 , 2.75 : 1.33e-02 4.73e-02 6.33e-02 7.32e-02 8.42e-02 9.03e-02 8.96e-02 9.70e-02 9.85e-02 9.70e-02 9.11e-02 8.82e-02 5.15e-02 1.33e-02 2.22e-03 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
2.75 , 3.00 : 1.24e-02 4.26e-02 6.39e-02 6.07e-02 7.64e-02 7.85e-02 7.65e-02 8.94e-02 8.22e-02 9.07e-02 7.85e-02 9.12e-02 7.99e-02 5.68e-02 1.73e-02 2.87e-03 1.96e-04 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
3.00 , 3.25 : 9.59e-03 3.89e-02 4.99e-02 6.57e-02 7.29e-02 7.41e-02 7.49e-02 6.99e-02 7.03e-02 7.17e-02 7.88e-02 8.86e-02 7.17e-02 8.43e-02 5.36e-02 2.25e-02 1.60e-03 1.02e-03 2.07e-04 0.00e+00 0.00e+00 0.00e+00
3.25 , 3.50 : 1.13e-02 3.18e-02 4.75e-02 5.72e-02 6.86e-02 5.98e-02 6.73e-02 6.65e-02 7.17e-02 7.11e-02 7.16e-02 7.52e-02 7.05e-02 6.96e-02 7.21e-02 5.84e-02 2.54e-02 3.94e-03 5.20e-04 0.00e+00 0.00e+00 0.00e+00
3.50 , 3.75 : 9.77e-03 2.80e-02 4.52e-02 5.28e-02 5.60e-02 5.90e-02 5.98e-02 6.27e-02 6.37e-02 6.34e-02 6.44e-02 6.55e-02 6.11e-02 7.20e-02 7.43e-02 6.67e-02 5.89e-02 3.01e-02 6.12e-03 2.20e-04 2.12e-04 0.00e+00
3.75 , 4.00 : 8.02e-03 2.75e-02 4.41e-02 4.81e-02 5.10e-02 6.43e-02 5.88e-02 5.61e-02 5.59e-02 6.36e-02 6.08e-02 5.79e-02 6.21e-02 5.52e-02 6.35e-02 6.59e-02 5.98e-02 5.78e-02 3.26e-02 5.67e-03 1.20e-03 0.00e+00
4.00 , 4.25 : 4.22e-03 3.83e-02 4.79e-02 3.29e-02 5.82e-02 5.78e-02 5.93e-02 4.00e-02 5.36e-02 6.29e-02 6.10e-02 6.48e-02 4.95e-02 5.73e-02 5.16e-02 5.39e-02 5.38e-02 5.09e-02 6.21e-02 3.55e-02 4.09e-03 4.54e-04
4.25 , 4.50 : 9.10e-03 2.74e-02 4.54e-02 3.60e-02 4.90e-02 6.69e-02 5.05e-02 5.33e-02 5.63e-02 6.52e-02 5.62e-02 5.12e-02 5.17e-02 6.15e-02 3.76e-02 5.48e-02 3.52e-02 5.44e-02 4.60e-02 4.89e-02 3.28e-02 1.01e-02
4.50 , 4.75 : 3.71e-03 1.34e-02 3.86e-02 5.30e-02 3.04e-02 5.23e-02 5.68e-02 5.31e-02 5.85e-02 5.32e-02 6.57e-02 4.85e-02 4.81e-02 4.86e-02 4.11e-02 4.82e-02 3.43e-02 6.49e-02 7.48e-02 4.95e-02 2.46e-02 3.87e-02
4.75 , 5.00 : 9.87e-03 1.91e-02 3.90e-02 3.26e-02 3.34e-02 5.06e-02 5.05e-02 3.13e-02 3.39e-02 5.66e-02 6.83e-02 7.06e-02 5.89e-02 3.52e-02 5.17e-02 6.75e-02 6.13e-02 3.75e-02 5.33e-02 5.83e-02 4.41e-02 3.64e-02
5.00 , 5.25 : 1.68e-02 2.33e-02 5.65e-02 4.93e-02 3.80e-02 5.26e-02 4.99e-02 6.15e-02 5.52e-02 7.25e-02 4.74e-02 4.02e-02 5.30e-02 6.22e-02 6.67e-02 5.50e-02 3.37e-02 3.72e-02 4.54e-02 2.60e-02 3.87e-02 1.89e-02
5.25 , 5.50 : 1.07e-02 2.00e-02 2.47e-02 4.62e-02 4.26e-02 3.80e-02 2.08e-02 8.11e-02 6.54e-02 2.46e-02 5.58e-02 5.75e-02 6.57e-02 4.24e-02 5.44e-02 4.33e-02 5.75e-02 3.56e-02 7.88e-02 3.74e-02 5.25e-02 4.51e-02
5.50 , 5.75 : 9.38e-03 2.06e-02 3.21e-02 4.04e-02 3.90e-02 7.56e-02 6.42e-02 4.45e-02 4.38e-02 4.93e-02 6.74e-02 3.59e-02 6.15e-02 6.24e-02 5.03e-02 4.25e-02 6.57e-02 4.14e-02 4.50e-02 4.60e-02 2.78e-02 3.54e-02
5.75 , 6.00 : 9.02e-03 3.89e-02 4.15e-02 2.91e-02 4.17e-02 5.26e-02 4.56e-02 4.59e-02 3.84e-02 6.19e-02 5.35e-02 6.33e-02 4.64e-02 6.30e-02 4.30e-02 5.59e-02 5.58e-02 3.57e-02 4.96e-02 5.27e-02 5.27e-02 2.38e-02

0 comments on commit ef573f8

Please sign in to comment.