From bf87a9b887828ebd8c836fd7ea9f7f99a6a4ea3f Mon Sep 17 00:00:00 2001
From: Johannes Schumann <johannes.schumann@fau.de>
Date: Mon, 17 Jan 2022 11:28:07 +0100
Subject: [PATCH] Update w2list fields
---
km3buu/output.py | 36 +++++++++++++++++++++++++++++++-----
1 file changed, 31 insertions(+), 5 deletions(-)
diff --git a/km3buu/output.py b/km3buu/output.py
index bc978f8..1c0802e 100644
--- a/km3buu/output.py
+++ b/km3buu/output.py
@@ -23,6 +23,7 @@ import uproot
from scipy.interpolate import UnivariateSpline, interp1d
from scipy.spatial.transform import Rotation
import scipy.constants as constants
+from scipy.optimize import curve_fit
import mendeleev
from datetime import datetime
@@ -230,19 +231,18 @@ W2LIST_LOOKUP = {
"COLUMN_DEPTH": 3,
"P_EARTH": 4,
"WATER_INT_LEN": 5,
- "P_SCALE": 6,
"BX": 7,
"BY": 8,
"ICHAN": 9,
"CC": 10,
- "DISTAMAX": 11,
- "WATERXSEC": 12,
"XSEC": 13,
+ "DXSEC": 14,
"TARGETA": 15,
"TARGETZ": 16,
"VERINCAN": 17,
"LEPINCAN": 18,
- "GIBUU_WEIGHT": 19
+ "GIBUU_WEIGHT": 23,
+ "GIBUU_SCAT_TYPE": 24
}
W2LIST_LENGTH = max(W2LIST_LOOKUP.values()) + 1
@@ -300,9 +300,11 @@ class GiBUUOutput:
self._min_energy = np.nan
self._max_energy = np.nan
self._generated_events = -1
+ self._flux_index = np.nan
try:
self._read_flux_file()
+ self._determine_flux_index()
except OSError:
self._read_single_energy()
@@ -577,6 +579,22 @@ class GiBUUOutput:
def generated_events(self):
return self._generated_events
+ def _determine_flux_index(self):
+ def fluxfunc(x, a, b):
+ return a * x**b
+
+ mask = ~np.isclose(self.flux_data["events"], 0)
+ popt, pcov = curve_fit(fluxfunc,
+ self.flux_data["energy"][mask],
+ self.flux_data["events"][mask],
+ p0=[1, -1])
+
+ self._flux_index = popt[1]
+
+ @property
+ def flux_index(self):
+ return self._flux_index
+
def write_detector_file(gibuu_output,
ofile="gibuu.offline.root",
@@ -703,6 +721,7 @@ def write_detector_file(gibuu_output,
# Event Information (w2list)
evt.w2list.resize(W2LIST_LENGTH)
evt.w2list[W2LIST_LOOKUP["PS"]] = global_generation_weight
+ evt.w2list[W2LIST_LOOKUP["EG"]] = gibuu_output.flux_index
evt.w2list[W2LIST_LOOKUP["XSEC_MEAN"]] = mean_xsec_func(event.lepIn_E)
evt.w2list[W2LIST_LOOKUP["XSEC"]] = event.xsec
evt.w2list[W2LIST_LOOKUP["TARGETA"]] = gibuu_output.A
@@ -710,10 +729,17 @@ def write_detector_file(gibuu_output,
evt.w2list[W2LIST_LOOKUP["BX"]] = bjorkenx[mc_event_id]
evt.w2list[W2LIST_LOOKUP["BY"]] = bjorkeny[mc_event_id]
evt.w2list[W2LIST_LOOKUP["CC"]] = ichan
- evt.w2list[W2LIST_LOOKUP["ICHAN"]] = event.evType
+ evt.w2list[W2LIST_LOOKUP["ICHAN"]] = SCATTERING_TYPE_TO_GENIE[
+ event.evType]
evt.w2list[W2LIST_LOOKUP["VERINCAN"]] = 1
evt.w2list[W2LIST_LOOKUP["LEPINCAN"]] = 1
evt.w2list[W2LIST_LOOKUP["GIBUU_WEIGHT"]] = event.weight
+ evt.w2list[W2LIST_LOOKUP["GIBUU_SCAT_TYPE"]] = event.evType
+ #TODO
+ evt.w2list[W2LIST_LOOKUP["DXSEC"]] = np.nan
+ evt.w2list[W2LIST_LOOKUP["COLUMN_DEPTH"]] = np.nan
+ evt.w2list[W2LIST_LOOKUP["P_EARTH"]] = np.nan
+ evt.w2list[W2LIST_LOOKUP["WATER_INT_LEN"]] = np.nan
# Vertex Position
vtx_pos = np.array(geometry.random_pos())
--
GitLab