diff --git a/km3buu/config.py b/km3buu/config.py
index b9c7a0c6eff1fb4d6f5b4aa25b3cdb1e4895c162..9059d6b94d377e4079c3dd688ca87b5138fb391f 100644
--- a/km3buu/config.py
+++ b/km3buu/config.py
@@ -45,6 +45,7 @@ GSEAGEN_MEDIA_COMPOSITION_FILE = "MediaComposition.xml"
class Config(object):
+
def __init__(self, config_path=CONFIG_PATH):
self.config = ConfigParser()
self._config_path = config_path
diff --git a/km3buu/geometry.py b/km3buu/geometry.py
index 1ae11351e13f342f8c686177d149f0eba03f1466..42378c3b4d6457f2f0530fa009a08186c1db1131 100644
--- a/km3buu/geometry.py
+++ b/km3buu/geometry.py
@@ -19,6 +19,7 @@ class DetectorVolume(ABC):
"""
Detector geometry class
"""
+
def __init__(self):
self._volume = -1.0
self._coord_origin = (0., 0., 0.)
@@ -82,6 +83,7 @@ class CanVolume(DetectorVolume):
zmax: float [m] (default: 476.5)
Cylinder top z position
"""
+
def __init__(self, radius=403.4, zmin=0.0, zmax=476.5):
super().__init__()
self._radius = radius
@@ -118,6 +120,7 @@ class SphericalVolume(DetectorVolume):
Coordinate center of the sphere
(x, y, z)
"""
+
def __init__(self, radius, coord_origin=(0, 0, 0)):
super().__init__()
self._radius = radius
diff --git a/km3buu/jobcard.py b/km3buu/jobcard.py
index a03f217c936a97dfac60369874d56046135eb491..15c99f6b7b1f26d9eec1910273785ca2517022c5 100644
--- a/km3buu/jobcard.py
+++ b/km3buu/jobcard.py
@@ -58,6 +58,7 @@ class Jobcard(f90nml.Namelist):
input_path: str
The input path pointing to the GiBUU lookup data which should be used
"""
+
def __init__(self,
*args,
filename=DEFAULT_JOBCARD_FILENAME,
diff --git a/km3buu/output.py b/km3buu/output.py
index cbd51c2c9bdd61ca6ffd529495999bf74a40eda6..61e5f75391d57e5c87bcbb4086b1e6dc372c5544 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
@@ -79,7 +80,7 @@ FLUX_INFORMATION_DTYPE = np.dtype([("energy", np.float64),
("flux", np.float64),
("events", np.float64)])
-EVENT_TYPE = {
+SCATTERING_TYPE = {
1: "QE",
2: "P33(1232)",
3: "P11(1440)",
@@ -119,6 +120,46 @@ EVENT_TYPE = {
37: "2pi background",
}
+SCATTERING_TYPE_TO_GENIE = {
+ 1: 1, # QE -> kScQuasiElastic
+ 2: 4, # P33(1232) -> kScResonant
+ 3: 4, # P11(1440) -> kScResonant
+ 4: 4, # S11(1535) -> kScResonant
+ 5: 4, # S11(1650) -> kScResonant
+ 6: 4, # S11(2090) -> kScResonant
+ 7: 4, # D13(1520) -> kScResonant
+ 8: 4, # D13(1700) -> kScResonant
+ 9: 4, # D13(2080) -> kScResonant
+ 10: 4, # D15(1675) -> kScResonant
+ 11: 4, # G17(2190) -> kScResonant
+ 12: 4, # P11(1710) -> kScResonant
+ 13: 4, # P11(2100) -> kScResonant
+ 14: 4, # P13(1720) -> kScResonant
+ 15: 4, # P13(1900) -> kScResonant
+ 16: 4, # F15(1680) -> kScResonant
+ 17: 4, # F15(2000) -> kScResonant
+ 18: 4, # F17(1990) -> kScResonant
+ 19: 4, # S31(1620) -> kScResonant
+ 20: 4, # S31(1900) -> kScResonant
+ 21: 4, # D33(1700) -> kScResonant
+ 22: 4, # D33(1940) -> kScResonant
+ 23: 4, # D35(1930) -> kScResonant
+ 24: 4, # D35(2350) -> kScResonant
+ 25: 4, # P31(1750) -> kScResonant
+ 26: 4, # P31(1910) -> kScResonant
+ 27: 4, # P33(1600) -> kScResonant
+ 28: 4, # P33(1920) -> kScResonant
+ 29: 4, # F35(1750) -> kScResonant
+ 30: 4, # F35(1905) -> kScResonant
+ 31: 4, # F37(1950) -> kScResonant
+ 32: 0, # pi neutron-background -> kScNull
+ 33: 0, # pi proton-background -> kScNull
+ 34: 3, # DIS -> kScDeepInelastic
+ 35: 0, # 2p2h QE -> kScNull
+ 36: 0, # 2p2h Delta -> kScNull
+ 37: 0, # 2pi background -> kScNull
+}
+
ROOTTUPLE_KEY = "RootTuple"
EMPTY_KM3NET_HEADER_DICT = {
@@ -190,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
@@ -234,6 +274,7 @@ def read_nu_abs_xsection(filepath):
class GiBUUOutput:
+
def __init__(self, data_dir):
"""
Class for parsing GiBUU output files
@@ -260,9 +301,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()
@@ -537,6 +580,26 @@ class GiBUUOutput:
def generated_events(self):
return self._generated_events
+ def _determine_flux_index(self):
+
+ 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]
+
+ @property
+ def flux_index(self):
+ return self._flux_index
+
def write_detector_file(gibuu_output,
ofile="gibuu.offline.root",
@@ -663,6 +726,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
@@ -670,10 +734,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())
diff --git a/km3buu/tests/test_config.py b/km3buu/tests/test_config.py
index 18e2a13ec87d2122279a0a15a349da9e2a446e9c..65436a054e56e37945cb31052ebd130f6d406112 100644
--- a/km3buu/tests/test_config.py
+++ b/km3buu/tests/test_config.py
@@ -32,6 +32,7 @@ path=/tmp/gseagen
class TestConfig(unittest.TestCase):
+
def setUp(self):
self.cfg_tmpfile = NamedTemporaryFile(delete=False)
self.mock_image_file = NamedTemporaryFile(delete=False)
diff --git a/km3buu/tests/test_ctrl.py b/km3buu/tests/test_ctrl.py
index 34d4faa8e300471b005b04c82637d8099384ca06..23090173209c3eebc8e86f3f51d9445af8a65862 100644
--- a/km3buu/tests/test_ctrl.py
+++ b/km3buu/tests/test_ctrl.py
@@ -25,6 +25,7 @@ TESTDATA_DIR = data_path("gibuu")
class TestCTRLbyJobcardFile(unittest.TestCase):
+
def setUp(self):
self.filename = join(TESTDATA_DIR, "km3net_testdata.job")
self.output_dir = TemporaryDirectory()
@@ -48,6 +49,7 @@ class TestCTRLbyJobcardFile(unittest.TestCase):
class TestCTRLbyJobcardObject(unittest.TestCase):
+
def setUp(self):
log = get_logger("ctrl.py")
log.setLevel("INFO")
diff --git a/km3buu/tests/test_environment.py b/km3buu/tests/test_environment.py
index 247690b1e66adda9d5c3c93be01c24a19f5738d8..b5d76bf3b4819b7aecf1d7fe1452ee86ee2d1c58 100644
--- a/km3buu/tests/test_environment.py
+++ b/km3buu/tests/test_environment.py
@@ -19,6 +19,7 @@ from km3buu import DOCKER_URL, IMAGE_NAME
class TestBuild(unittest.TestCase):
+
def test_wrong_dir_path(self):
wrong_path = "foobar"
try:
diff --git a/km3buu/tests/test_geometry.py b/km3buu/tests/test_geometry.py
index 292cac2dd3b709a8279c02976db9148471de7b7a..10a22bda09609c05816a60379fdbd402080f8b2a 100644
--- a/km3buu/tests/test_geometry.py
+++ b/km3buu/tests/test_geometry.py
@@ -17,12 +17,14 @@ import numpy as np
class TestGeneralGeometry(unittest.TestCase):
+
def test_abstract_init(self):
with self.assertRaises(TypeError) as ctx:
d = DetectorVolume()
class TestSphere(unittest.TestCase):
+
def setUp(self):
self.detector_geometry = SphericalVolume(20, (2, 2, 2))
@@ -44,6 +46,7 @@ class TestSphere(unittest.TestCase):
class TestCan(unittest.TestCase):
+
def setUp(self):
self.detector_geometry = CanVolume()
diff --git a/km3buu/tests/test_jobcard.py b/km3buu/tests/test_jobcard.py
index 6fc45b08949ab74b1a53076233d43e7375dec8fa..69dc53951bfd9945d4fe165bf57434752013838d 100644
--- a/km3buu/tests/test_jobcard.py
+++ b/km3buu/tests/test_jobcard.py
@@ -20,6 +20,7 @@ from tempfile import TemporaryFile, TemporaryDirectory
class TestJobcard(unittest.TestCase):
+
def setUp(self):
self.test_jobcard = Jobcard()
# Insert some test elements
@@ -49,6 +50,7 @@ class TestJobcard(unittest.TestCase):
class TestNeutrinoEnergyRangeJobcard(unittest.TestCase):
+
def setUp(self):
self.test_fluxfile = TemporaryFile()
self.test_Z = np.random.randint(1, 100)
@@ -99,6 +101,7 @@ class TestNeutrinoEnergyRangeJobcard(unittest.TestCase):
class TestNeutrinoSingleEnergyJobcard(unittest.TestCase):
+
def setUp(self):
self.test_fluxfile = TemporaryFile()
self.test_Z = np.random.randint(1, 100)
@@ -133,6 +136,7 @@ class TestNeutrinoSingleEnergyJobcard(unittest.TestCase):
class TestJobcardSeed(unittest.TestCase):
+
def setUp(self):
jc = generate_neutrino_jobcard(100,
"CC",
diff --git a/km3buu/tests/test_output.py b/km3buu/tests/test_output.py
index cb8de46664731163fc919035027aa4c3bb0cd3f4..02c9aba5e0206dd9c8658e843c976c9afd9e9609 100644
--- a/km3buu/tests/test_output.py
+++ b/km3buu/tests/test_output.py
@@ -37,6 +37,7 @@ except ModuleNotFoundError:
class TestXSection(unittest.TestCase):
+
def test_xsection_all(self):
filename = join(TESTDATA_DIR, XSECTION_FILENAMES["all"])
xsection = read_nu_abs_xsection(filename)
@@ -47,6 +48,7 @@ class TestXSection(unittest.TestCase):
class TestGiBUUOutput(unittest.TestCase):
+
def setup_class(self):
self.output = GiBUUOutput(TESTDATA_DIR)
@@ -73,6 +75,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(
@@ -88,6 +93,7 @@ class TestGiBUUOutput(unittest.TestCase):
@pytest.mark.skipif(not KM3NET_LIB_AVAILABLE,
reason="KM3NeT dataformat required")
class TestAANET(unittest.TestCase):
+
def setUp(self):
output = GiBUUOutput(TESTDATA_DIR)
datafile = NamedTemporaryFile(suffix=".root")
diff --git a/km3buu/tests/test_propagation.py b/km3buu/tests/test_propagation.py
index 27d17044e99561100afe8ee46058e0732a1d1420..b1cbea92d9438be21a0338bf6da515b32b386950 100644
--- a/km3buu/tests/test_propagation.py
+++ b/km3buu/tests/test_propagation.py
@@ -28,6 +28,7 @@ pp.RandomGenerator.get().set_seed(1234)
@pytest.mark.skip(reason="CI boost lib problem")
class TestTauPropagation(unittest.TestCase):
+
def setUp(self):
data = ak.Array({
"lepOut_E": [