From d9c2d4d67872ddc611cf9c54b0f0939315adbb63 Mon Sep 17 00:00:00 2001
From: Johannes Schumann <johannes.schumann@fau.de>
Date: Thu, 26 Oct 2023 00:05:15 +0000
Subject: [PATCH] Flux Index Fit
---
CHANGELOG.rst | 1 +
km3buu/output.py | 32 +++++++++++++++++++-------------
km3buu/tests/test_output.py | 2 +-
3 files changed, 21 insertions(+), 14 deletions(-)
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index d6f1c0a..d32846f 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,5 +1,6 @@
Unreleased changes
------------------
+* Improved curve for flux index
v1.1.2
----------------------------
diff --git a/km3buu/output.py b/km3buu/output.py
index 973b9ab..9f3747f 100644
--- a/km3buu/output.py
+++ b/km3buu/output.py
@@ -16,7 +16,7 @@ import re
import numpy as np
from io import StringIO
from os import listdir, environ
-from os.path import isfile, join, abspath
+from os.path import isfile, join, abspath, exists
from tempfile import TemporaryDirectory
import awkward as ak
import uproot
@@ -237,10 +237,9 @@ class GiBUUOutput:
self._generated_events = -1
self._flux_index = np.nan
- try:
- self._read_flux_file()
+ if self._read_flux_file():
self._determine_flux_index()
- except OSError:
+ else:
self._read_single_energy()
def _read_root_output(self):
@@ -285,6 +284,8 @@ class GiBUUOutput:
def _read_flux_file(self):
fpath = join(self._data_path, FLUXDESCR_FILENAME)
+ if not exists(fpath):
+ return False
self.flux_data = np.loadtxt(fpath,
dtype=FLUX_INFORMATION_DTYPE,
usecols=(
@@ -297,6 +298,7 @@ class GiBUUOutput:
self._energy_min = np.min(self.flux_data["energy"])
self._energy_max = np.max(self.flux_data["energy"])
self._generated_events = int(np.sum(self.flux_data["events"]))
+ return True
def _event_xsec(self, root_tupledata):
weights = np.array(root_tupledata.weight)
@@ -550,15 +552,19 @@ class GiBUUOutput:
def fluxfunc(x, a, b):
return a * x**b
- 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["flux"][mask],
- p0=[1, -1])
- self._flux_index = popt[1]
+ energy_mask = self.flux_data["flux"] > 0
+ lower_limit = np.min(self.flux_data["energy"][energy_mask]) * 1.2
+ upper_limit = np.max(self.flux_data["energy"][energy_mask]) * 0.8
+ mask = (self.flux_data["energy"]
+ > lower_limit) & (self.flux_data["energy"] < upper_limit)
+ try:
+ popt, pcov = curve_fit(fluxfunc,
+ self.flux_data["energy"][mask],
+ self.flux_data["flux"][mask],
+ p0=[1, -1])
+ self._flux_index = popt[1]
+ except:
+ self._flux_index = np.nan
@property
def flux_index(self):
diff --git a/km3buu/tests/test_output.py b/km3buu/tests/test_output.py
index 553bc56..a8a0de9 100644
--- a/km3buu/tests/test_output.py
+++ b/km3buu/tests/test_output.py
@@ -77,7 +77,7 @@ class TestGiBUUOutput(unittest.TestCase):
assert self.output.A == 16
def test_flux_index(self):
- assert np.isclose(self.output.flux_index, -0.004812255, rtol=1e-3)
+ assert np.isclose(self.output.flux_index, -0.00576924, rtol=1e-3)
def test_w2weights(self):
w2 = self.output.w2weights(123.0, 2.6e28, 4 * np.pi)
--
GitLab