From 6bbb81c3ec78998ab95a9615cfecf62475a2819f Mon Sep 17 00:00:00 2001
From: Johannes Schumann <johannes.schumann@fau.de>
Date: Mon, 17 Jan 2022 13:47:19 +0100
Subject: [PATCH] Refine flux determination method
---
km3buu/output.py | 7 +++++--
km3buu/tests/test_output.py | 3 +++
2 files changed, 8 insertions(+), 2 deletions(-)
diff --git a/km3buu/output.py b/km3buu/output.py
index 1c0802e..31766f0 100644
--- a/km3buu/output.py
+++ b/km3buu/output.py
@@ -583,10 +583,13 @@ class GiBUUOutput:
def fluxfunc(x, a, b):
return a * x**b
- mask = ~np.isclose(self.flux_data["events"], 0)
+ lower_limit = np.exp(np.log(np.max(self.flux_data["flux"])) * 0.2)
+ upper_limit = np.exp(np.log(np.max(self.flux_data["flux"])) * 0.8)
+ mask = (self.flux_data["flux"] > lower_limit) & (self.flux_data["flux"]
+ < upper_limit)
popt, pcov = curve_fit(fluxfunc,
self.flux_data["energy"][mask],
- self.flux_data["events"][mask],
+ self.flux_data["flux"][mask],
p0=[1, -1])
self._flux_index = popt[1]
diff --git a/km3buu/tests/test_output.py b/km3buu/tests/test_output.py
index cb8de46..e6ff083 100644
--- a/km3buu/tests/test_output.py
+++ b/km3buu/tests/test_output.py
@@ -73,6 +73,9 @@ class TestGiBUUOutput(unittest.TestCase):
assert self.output.Z == 8
assert self.output.A == 16
+ def test_flux_index(self):
+ assert np.isclose(self.output.flux_index, -0.904, rtol=1e-3)
+
def test_w2weights(self):
w2 = self.output.w2weights(123.0, 2.6e28, 4 * np.pi)
np.testing.assert_array_almost_equal(
--
GitLab