Skip to content

Commit

Permalink
Progress for simulate.
Browse files Browse the repository at this point in the history
  • Loading branch information
@robertdstein
robertdstein committed May 9, 2020
1 parent 8d77181 commit 77cf54f
Show file tree
Hide file tree
Showing 7 changed files with 249 additions and 20 deletions.
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
21 changes: 9 additions & 12 deletions flarestack/data/simulate/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,8 @@ def add_subseason(self, season):
def add_sim_season(self, name, sim_season_f):
self.sim_seasons[name] = sim_season_f

def set_sim_params(self, name, bkg_time_pdf_dict, bkg_flux_norm,
bkg_e_pdf_dict, **kwargs):
def set_sim_params(self, name, bkg_flux_model,
**kwargs):

if name not in self.sim_seasons.keys():
raise KeyError("SimSeasonClass {0} not found. The following "
Expand All @@ -45,35 +45,31 @@ def set_sim_params(self, name, bkg_time_pdf_dict, bkg_flux_norm,
"function.".format(name, self.sim_seasons.keys()))


self.seasons[name] = self.sim_seasons[name](
bkg_time_pdf_dict, bkg_flux_norm, bkg_e_pdf_dict, **kwargs)
self.seasons[name] = self.sim_seasons[name](bkg_flux_model, **kwargs)


class SimSeason(SeasonWithoutMC):

def __init__(self, season_name, sample_name, pseudo_mc_path,
event_dtype, effective_area_f, load_angular_resolution,
bkg_t_pdf_dict, bkg_flux_norm, bkg_e_pdf_dict,
bkg_flux_model,
energy_proxy_map, **kwargs):

self.bkg_t_pdf_dict = bkg_t_pdf_dict
self.event_dtype = event_dtype

self.bkg_flux_norm = bkg_flux_norm
self.bkg_e_pdf_dict = bkg_e_pdf_dict
self.bkg_flux_model = bkg_flux_model

self.load_angular_resolution = load_angular_resolution
self.angular_res_f = self.load_angular_resolution()

self.base_dataset_path = sim_dataset_dir_path(
sample_name, season_name, bkg_flux_norm, bkg_e_pdf_dict)
sample_name, season_name, bkg_flux_model)

exp_path = self.dataset_path()

SeasonWithoutMC.__init__(
self, season_name, sample_name, exp_path, pseudo_mc_path, **kwargs
)
self.bkg_energy_pdf = EnergyPDF.create(bkg_e_pdf_dict)

self.energy_proxy_map = energy_proxy_map

Expand All @@ -85,7 +81,7 @@ def load_energy_proxy_mapping(self):
return self.energy_proxy_map

def build_time_pdf_dict(self):
return self.bkg_t_pdf_dict
return self.bkg_flux_model.build_time_pdf_dict()

def check_sim(self, resimulate=False, **kwargs):

Expand Down Expand Up @@ -124,7 +120,8 @@ def generate_sim_data(self, fluence):
return

def get_time_integrated_flux(self):

return k_to_flux(
self.bkg_flux_norm * self.get_time_pdf().effective_injection_time())
self.bkg_flux_model.get_norm() * self.get_time_pdf().effective_injection_time())


195 changes: 195 additions & 0 deletions flarestack/data/simulate/potemkin.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,195 @@
import numpy as np
import random
from scipy.interpolate import interp1d
import logging
from flarestack.data.public import icecube_ps_3_year
from flarestack.core.energy_pdf import EnergyPDF
from flarestack.core.time_pdf import TimePDF
from flarestack.data.simulate import SimSeason, SimDataset
from flarestack.shared import flux_to_k

class BackgroundFluxModel:


def __init__(self, flux_norm, bkg_time_pdf_dict):
self.flux_norm = flux_to_k(flux_norm)
self.bkg_time_pdf_dict = bkg_time_pdf_dict

def get_norm(self):
return self.flux_norm

def flux_model_f(self, e, sindec):
return NotImplementedError

def flux_range(self):
return NotImplementedError

def unique_name(self):
return self.flux_norm

def build_time_pdf_dict(self):
return self.bkg_time_pdf_dict

class IdealBackgroundFluxModel(BackgroundFluxModel):

def __init__(self, flux_norm, bkg_time_pdf_dict):
BackgroundFluxModel.__init__(self, flux_norm, bkg_time_pdf_dict)
self.atmo = self.atmo_flux()

@staticmethod
def atmo_flux():
bkg_e_pdf_dict = {
"energy_pdf_name": "power_law",
"gamma": 3.7,
"e_min_gev": 100.,
"e_max_gev": 10.**7
}
return EnergyPDF.create(bkg_e_pdf_dict)

def flux_model_f(self, e, sindec):
return self.atmo.f(e)

def flux_range(self):
return self.atmo.e_min, self.atmo.e_max


class PotemkinSeason(SimSeason):

def __init__(self, season_name, sample_name, pseudo_mc_path,
event_dtype, effective_area_f, load_angular_resolution,
bkg_flux_model,
energy_proxy_map, sin_dec_bins, log_e_bins, **kwargs):

self.log_e_bins = log_e_bins
self.sin_dec_bins = sin_dec_bins

SimSeason.__init__(
self, season_name, sample_name, pseudo_mc_path,
event_dtype, effective_area_f,
load_angular_resolution, bkg_flux_model,
energy_proxy_map, **kwargs
)

def generate_sim_data(self, fluence):
logging.info("Simulating events:")
sim_events = np.empty((0,),
dtype=self.event_dtype)

for i, lower_sin_dec in enumerate(self.sin_dec_bins[:-1]):
upper_sin_dec = self.sin_dec_bins[i + 1]
new_events = self.simulate_dec_range(
fluence,lower_sin_dec, upper_sin_dec)

logging.info("Simulated {0} events between sin(dec)={1} and "
"sin(dec)={2}".format(
len(new_events), lower_sin_dec, upper_sin_dec))

# Joins the new events to the signal events
sim_events = np.concatenate((sim_events, new_events))

sim_events = np.array(sim_events)

logging.info("Simulated {0} events in total".format(len(sim_events)))

return sim_events

def simulate_dec_range(self, fluence, lower_sin_dec, upper_sin_dec):
mean_sin_dec = 0.5 * (lower_sin_dec + upper_sin_dec)
solid_angle = 2 * np.pi * (upper_sin_dec - lower_sin_dec)
sim_fluence = fluence * solid_angle # GeV^-1 cm^-2

def source_eff_area(e):
return self.effective_area_f(
np.log10(e), mean_sin_dec) * self.bkg_flux_model.flux_model_f(e, mean_sin_dec)

lower, upper = self.bkg_flux_model.flux_range()

logging.debug(f"Simulating between {lower:.2g} GeV and {upper:.2g} GeV")

int_eff_a = EnergyPDF.integrate(source_eff_area, lower, upper)

# Effective areas are given in m2, but flux is in per cm2

int_eff_a *= 10**4

n_exp = sim_fluence * int_eff_a

n_sim = np.random.poisson(n_exp)

new_events = np.empty((n_sim,), dtype=self.event_dtype)

new_events["ra"] = [random.random() * 2 * np.pi for _ in range(n_sim)]

new_events["sinDec"] = [
lower_sin_dec + (random.random() * (upper_sin_dec - lower_sin_dec))
for _ in range(n_sim)
]
new_events["dec"] = np.arcsin(new_events["sinDec"])
new_events["time"] = self.get_time_pdf().simulate_times([], n_sim)

fluence_ints, log_e_range = \
EnergyPDF.piecewise_integrate(
source_eff_area, lower, upper
)
fluence_ints = fluence_ints

fluence_ints = np.array(fluence_ints)
fluence_ints /= np.sum(fluence_ints)

fluence_cumulative = [np.sum(fluence_ints[:i])
for i, _ in enumerate(fluence_ints)]

fluence_cumulative = [0.] + fluence_cumulative + [1.]

log_e_range = list(log_e_range) + [log_e_range[-1]]

sim_true_e = interp1d(fluence_cumulative, log_e_range)

true_e_vals = np.array(
[10**sim_true_e(random.random()) for _ in range(n_sim)])

new_events["logE"] = self.energy_proxy_map(true_e_vals)

new_events["sigma"] = self.angular_res_f(new_events["logE"]).copy()
new_events["raw_sigma"] = new_events["sigma"].copy()

return new_events

potemkin_dataset = SimDataset()

for (name, season) in icecube_ps_3_year.get_seasons().items():

def ideal_energy_proxy(e):
return np.log10(e)

def wrapper_f(bkg_flux_model, energy_proxy_map=None, sim_name=None, **kwargs):

if np.logical_and(energy_proxy_map is None, sim_name is None):
energy_proxy_map = ideal_energy_proxy
sim_name = "default"

if np.logical_and(energy_proxy_map != ideal_energy_proxy,
sim_name is None):
raise ValueError("Non-default energy proxy mapping was used, "
"but no unique sim_name was provided. Please "
"provide a unique 'sim_name' to describe this "
"simulation.")

sin_dec_bins = season.sin_dec_bins[season.sin_dec_bins > 0.]

sim_season = PotemkinSeason(
season_name=name,
sample_name="SimCube_{0}".format(sim_name),
pseudo_mc_path=season.pseudo_mc_path,
effective_area_f=season.load_effective_area(),
event_dtype=season.get_background_dtype(),
load_angular_resolution=season.load_angular_resolution,
bkg_flux_model=bkg_flux_model,
energy_proxy_map=energy_proxy_map,
sin_dec_bins=sin_dec_bins,
log_e_bins=season.log_e_bins,
**kwargs
)
return sim_season

potemkin_dataset.add_sim_season(name, wrapper_f)
5 changes: 2 additions & 3 deletions flarestack/data/simulate/simcube.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ def source_eff_area(e):
def ideal_energy_proxy(e):
return np.log10(e)

def wrapper_f(bkg_time_pdf_dict, bkg_flux_norm, bkg_e_pdf_dict,
def wrapper_f(bkg_time_pdf_dict, bkg_flux_model,
energy_proxy_map=None, sim_name=None, **kwargs):

if np.logical_and(energy_proxy_map is None, sim_name is None):
Expand All @@ -208,8 +208,7 @@ def wrapper_f(bkg_time_pdf_dict, bkg_flux_norm, bkg_e_pdf_dict,
event_dtype=season.get_background_dtype(),
load_angular_resolution=season.load_angular_resolution,
bkg_time_pdf_dict=bkg_time_pdf_dict,
bkg_flux_norm=bkg_flux_norm,
bkg_e_pdf_dict=bkg_e_pdf_dict,
bkg_flux_model=bkg_flux_model,
energy_proxy_map=energy_proxy_map,
sin_dec_bins=sin_dec_bins,
log_e_bins=season.log_e_bins,
Expand Down
7 changes: 3 additions & 4 deletions flarestack/shared.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,10 +195,9 @@ def limit_output_path(name):
return path


def sim_dataset_dir_path(sample_name, season_name, flux, e_pdf_dict):
return sim_dataset_dir + sample_name + "/" + \
season_name + '/' + str(deterministic_hash(e_pdf_dict)) + \
"/" + str(flux)
def sim_dataset_dir_path(sample_name, season_name, bkg_flux_model):
return f"{sim_dataset_dir}{sample_name}/" \
f"{season_name}/{bkg_flux_model.unique_name()}/"

def get_base_sob_plot_dir(season):
return dataset_plot_dir + "Signal_over_background/" + \
Expand Down

0 comments on commit 77cf54f

Please sign in to comment.