Resolve "ICHAN type"

Closes #6 (closed)

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