Interpolation rmf rebining

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:
@@ -122,6 +123,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_itp('86_II')
+
 rmf = {year: lru_cache()(rmf_factory(year)) for year in years}
 rmf['6y'] = rmf['86_II']
 

pyicecube10/interpolation_functions.py

@@ -0,0 +1,29 @@
+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)
+