Resolve "ICHAN type"

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

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

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,

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())

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)

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")

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:

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()
 

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",

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")

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": [