From b5a9658759840b3302a9c4f0e45ac9eb8dfafd33 Mon Sep 17 00:00:00 2001 From: Mikhail Date: Fri, 9 Feb 2024 15:37:24 +0200 Subject: [PATCH 1/4] added interpolation rebining --- pyicecube10/datareader.py | 20 +++++++++ pyicecube10/interpolation_functions.py | 61 ++++++++++++++++++++++++++ 2 files changed, 81 insertions(+) create mode 100644 pyicecube10/interpolation_functions.py diff --git a/pyicecube10/datareader.py b/pyicecube10/datareader.py index 908e106..c9a3202 100644 --- a/pyicecube10/datareader.py +++ b/pyicecube10/datareader.py @@ -6,6 +6,7 @@ from numpy.random import default_rng rng = default_rng() from functools import lru_cache +from .interpolation_functions import hist_rebin events = {} for year in years: @@ -75,6 +76,8 @@ smearing_raw[year] = smearing_raw['86_II'] +bins_Emu_min = smearing_raw['86_I']['log10(E/GeV)_min'].unique() +bins_Emu_max = smearing_raw['86_I']['log10(E/GeV)_max'].unique() def pdf_factory(iyear): if iyear in ('40', '59', '79', '86_I', '86_II'): def tmp(logEnu, dec, beta, dist): @@ -122,6 +125,23 @@ def tmp_rmf(dec, beta = 20, dist = None): return tmp_rmf else: return rmf_factory('86_II') + + +binning = np.append(bins_E_min, max(bins_E_max)) +# rebinning by interpolation +def rmf_factory_itp(jyear): + if jyear in ('40', '59', '79', '86_I', '86_II'): + + def tmp_rmf(dec, beta = 20, dist = None): + rebinned_pdfs = [] + for Enu in bins_Enu_mean: + pdf_current = pdf[jyear](Enu, dec, beta, dist) + rebinned_pdfs.append(hist_rebin(np.array(pdf_current), binning)) + return np.stack(rebinned_pdfs).transpose() + return tmp_rmf + else: + return rmf_factory('86_II') + rmf = {year: lru_cache()(rmf_factory(year)) for year in years} rmf['6y'] = rmf['86_II'] diff --git a/pyicecube10/interpolation_functions.py b/pyicecube10/interpolation_functions.py new file mode 100644 index 0000000..63ec957 --- /dev/null +++ b/pyicecube10/interpolation_functions.py @@ -0,0 +1,61 @@ +import numpy as np + + +def hist_rebin(hist, new_bins): + res = [] + amp = hist[:,2] + old_bins = np.append(hist[:,0], max(hist[:,1])) + bins = list(zip(np.insert(new_bins, 0, -np.inf), np.append(new_bins, np.inf))) + iterbins = list(enumerate(old_bins)) + flag = 1 + for ledge, redge in bins[1:-1]: + val = 0 + nledge = ledge + for i in iterbins[flag:]: + if ledge <= i[1] < redge: + val += amp[i[0] - 1] * (i[1] - nledge) / (i[1] - old_bins[i[0] - 1]) + nledge = i[1] + if redge <= old_bins[i[0] - 1]: + break + if redge <= i[1]: + val += amp[i[0] - 1] * (redge - nledge) / (i[1] - old_bins[i[0] - 1]) + flag = i[0] + break + if (type(amp[0]) is np.ndarray) and (val is 0): + res.append(np.zeros(len(amp[0]))) + else: + res.append(val) + return np.array(res) + + + +def hist_rebin_prop(amp, old_bins, new_bins): + res = [] + bins = list(zip(np.insert(new_bins, 0, -np.inf), np.append(new_bins, np.inf))) + iterbins = list(enumerate(old_bins)) + flag = 1 + for ledge, redge in bins[1:-1]: + val = 0 + prop = 0 + nledge = ledge + for i in iterbins[flag:]: + if ledge <= i[1] < redge: + val += amp[i[0] - 1] * (i[1] - nledge) / (i[1] - old_bins[i[0] - 1]) + prop += np.sum(amp[i[0] - 1]) * (i[1] - nledge) / (redge - ledge) + nledge = i[1] + if redge <= old_bins[i[0] - 1]: + break + if redge <= i[1]: + val += amp[i[0] - 1] * (redge - nledge) / (i[1] - old_bins[i[0] - 1]) + prop += np.sum(amp[i[0] - 1]) * (redge - nledge) / (redge - ledge) + flag = i[0] + break + if ((type(amp[0]) is np.ndarray) and (val is 0)) or (np.sum(val) < 1e-4): + res.append(np.zeros(len(amp[0]))) + else: + # print(np.sum(val), prop) + + res.append(val / np.sum(val) * prop) + + # print(res) + return np.array(res) From 9b7823f1bacde43da1cfc16b618dac782d77347b Mon Sep 17 00:00:00 2001 From: Mikhail Date: Fri, 9 Feb 2024 15:39:31 +0200 Subject: [PATCH 2/4] added interpolation rebining --- pyicecube10/interpolation_functions.py | 32 -------------------------- 1 file changed, 32 deletions(-) diff --git a/pyicecube10/interpolation_functions.py b/pyicecube10/interpolation_functions.py index 63ec957..9ea459f 100644 --- a/pyicecube10/interpolation_functions.py +++ b/pyicecube10/interpolation_functions.py @@ -27,35 +27,3 @@ def hist_rebin(hist, new_bins): res.append(val) return np.array(res) - - -def hist_rebin_prop(amp, old_bins, new_bins): - res = [] - bins = list(zip(np.insert(new_bins, 0, -np.inf), np.append(new_bins, np.inf))) - iterbins = list(enumerate(old_bins)) - flag = 1 - for ledge, redge in bins[1:-1]: - val = 0 - prop = 0 - nledge = ledge - for i in iterbins[flag:]: - if ledge <= i[1] < redge: - val += amp[i[0] - 1] * (i[1] - nledge) / (i[1] - old_bins[i[0] - 1]) - prop += np.sum(amp[i[0] - 1]) * (i[1] - nledge) / (redge - ledge) - nledge = i[1] - if redge <= old_bins[i[0] - 1]: - break - if redge <= i[1]: - val += amp[i[0] - 1] * (redge - nledge) / (i[1] - old_bins[i[0] - 1]) - prop += np.sum(amp[i[0] - 1]) * (redge - nledge) / (redge - ledge) - flag = i[0] - break - if ((type(amp[0]) is np.ndarray) and (val is 0)) or (np.sum(val) < 1e-4): - res.append(np.zeros(len(amp[0]))) - else: - # print(np.sum(val), prop) - - res.append(val / np.sum(val) * prop) - - # print(res) - return np.array(res) From dda103c72e68597923ba9fadf3d620acb1d2d64a Mon Sep 17 00:00:00 2001 From: Mikhail Date: Fri, 9 Feb 2024 15:40:26 +0200 Subject: [PATCH 3/4] added interpolation rebining --- pyicecube10/datareader.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyicecube10/datareader.py b/pyicecube10/datareader.py index c9a3202..4c4d87d 100644 --- a/pyicecube10/datareader.py +++ b/pyicecube10/datareader.py @@ -76,8 +76,6 @@ smearing_raw[year] = smearing_raw['86_II'] -bins_Emu_min = smearing_raw['86_I']['log10(E/GeV)_min'].unique() -bins_Emu_max = smearing_raw['86_I']['log10(E/GeV)_max'].unique() def pdf_factory(iyear): if iyear in ('40', '59', '79', '86_I', '86_II'): def tmp(logEnu, dec, beta, dist): From 567424a2dd6e4ed846fb3afcce24b34edc5288ac Mon Sep 17 00:00:00 2001 From: Mikhail Date: Tue, 13 Feb 2024 17:20:56 +0200 Subject: [PATCH 4/4] Fix rmf_factory_itp --- pyicecube10/datareader.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyicecube10/datareader.py b/pyicecube10/datareader.py index 4c4d87d..c066f17 100644 --- a/pyicecube10/datareader.py +++ b/pyicecube10/datareader.py @@ -138,7 +138,7 @@ def tmp_rmf(dec, beta = 20, dist = None): return np.stack(rebinned_pdfs).transpose() return tmp_rmf else: - return rmf_factory('86_II') + return rmf_factory_itp('86_II') rmf = {year: lru_cache()(rmf_factory(year)) for year in years} rmf['6y'] = rmf['86_II']