diff --git a/README.rst b/README.rst
index 68cc90b100fb212213975f1a6147f44d3246ccf0..8a891bdad4b73329b8e145add7a10d82f64010d8 100644
--- a/README.rst
+++ b/README.rst
@@ -64,6 +64,27 @@ This specific command runs all jobcards within the ``jobcards/examples`` folder
and stores the output inside the folder ``output``. The folder structure
is applied from the ``jobcards``\ folder.
+Theory
+------
+
+In order to retrieve correct results and provide correct KM3NeT weights (w2)
+the treatment of the GiBUU weights is an important step. A brief description
+of the GiBUU weights and how to calculate actual cross sections is given on the
+`GiBUU Homepage <https://gibuu.hepforge.org/trac/wiki/perWeight>`__ and
+a more detailed description of the calculation can be found in the `PhD Thesis
+of Tina Leitner <https://inspirehep.net/literature/849921>`__ in Chapter 8.3.
+As it is mentioned in the description of the output flux file in the
+`documentation <https://gibuu.hepforge.org/Documentation/code/init/neutrino/initNeutrino_f90.html#robo1685>`__ this is not taken somehow into account inside the weights.
+Following the description the GiBUU event weight can be converted to a binned
+cross section via
+
+.. math::
+ \frac{d\sigma}{E} = \frac{\sum_{i\in I_\text{bin}} w_i}{\Delta E}\cdot\frac{1}{E\Phi},
+
+where :math:`\Phi`__ is the simulated flux
+As the weights are given for each run individually the weight also has to be divided
+by the number of runs.
+
Tutorial
--------
The python framework is build around the GiBUU workflow, i.e. a jobcard is
diff --git a/km3buu/output.py b/km3buu/output.py
index adf34cb360f88281fff6bcb894bf052f619ff0d6..06140a036c9e8c7d477ce777cff10eeb9b4fb07a 100644
--- a/km3buu/output.py
+++ b/km3buu/output.py
@@ -49,6 +49,7 @@ FLUXDESCR_FILENAME = "neutrino_initialized_energyFlux.dat"
XSECTION_FILENAMES = {"all": "neutrino_absorption_cross_section_ALL.dat"}
SECONDS_PER_YEAR = 365.25 * 24 * 60 * 60
+SECONDS_WEIGHT_TIMESPAN = 1
PARTICLE_COLUMNS = ["E", "Px", "Py", "Pz", "barcode"]
EVENTINFO_COLUMNS = [
@@ -204,7 +205,7 @@ W2LIST_LOOKUP = {
"GIBUU_WEIGHT": 19
}
-W2LIST_LENGTH = len(W2LIST_LOOKUP)
+W2LIST_LENGTH = max(W2LIST_LOOKUP.values()) + 1
GIBUU_FIELDNAMES = [
'weight', 'barcode', 'Px', 'Py', 'Pz', 'E', 'evType', 'lepIn_E',
@@ -352,6 +353,15 @@ class GiBUUOutput:
bounds_error=False)
return lambda e: xsec_interp(e) * e
+ def global_generation_weight(self, solid_angle):
+ # I_E * I_theta * t_gen (* #NuTypes)
+ if self.flux_data is not None:
+ energy_phase_space = self.flux_interpolation.integral(
+ self._energy_min, self._energy_max)
+ else:
+ energy_phase_space = 1
+ return solid_angle * energy_phase_space * SECONDS_WEIGHT_TIMESPAN
+
def w2weights(self, volume, target_density, solid_angle):
"""
Calculate w2weights
@@ -378,7 +388,7 @@ class GiBUUOutput:
energy_factor = energy_phase_space * inv_gen_flux
else:
energy_factor = 1
- env_factor = volume * SECONDS_PER_YEAR
+ env_factor = volume * SECONDS_WEIGHT_TIMESPAN
retval = env_factor * solid_angle * energy_factor * xsec * 10**-42 * target_density
return retval
@@ -503,7 +513,11 @@ class GiBUUOutput:
tmp = fobj["RootTuple"].arrays()
retval = np.concatenate((retval, tmp))
if len(retval) == 0:
- retval = ak.Array(np.recarray((0,), dtype=list(zip(GIBUU_FIELDNAMES, len(GIBUU_FIELDNAMES)*[float]))))
+ retval = ak.Array(
+ np.recarray((0, ),
+ dtype=list(
+ zip(GIBUU_FIELDNAMES,
+ len(GIBUU_FIELDNAMES) * [float]))))
# Calculate additional information
retval["xsec"] = self._event_xsec(retval)
retval["Bx"] = GiBUUOutput.bjorken_x(retval)
@@ -605,14 +619,15 @@ def write_detector_file(gibuu_output,
tau_secondaries = propagate_lepton(event_data, np.sign(nu_type) * 15)
media = read_default_media_compositions()
- density = media["SeaWater"]["density"]
+ density = media["SeaWater"]["density"] # [g/cm^3]
element = mendeleev.element(gibuu_output.Z)
target = media["SeaWater"]["elements"][element.symbol]
- target_density = 1e3 * density * target[1]
- targets_per_volume = target_density * (1e3 * constants.Avogadro /
- target[0].atomic_weight)
+ target_density = 1e3 * density * target[1] # [kg/m^3]
+ targets_per_volume = target_density / target[
+ 0].atomic_weight / constants.atomic_mass
w2 = gibuu_output.w2weights(geometry.volume, targets_per_volume, 4 * np.pi)
+ global_generation_weight = gibuu_output.global_generation_weight(4 * np.pi)
head = ROOT.Head()
header_dct = EMPTY_KM3NET_HEADER_DICT.copy()
@@ -647,6 +662,7 @@ def write_detector_file(gibuu_output,
evt.w.push_back(-1.0) #w3 (= w2*flux)
# Event Information (w2list)
evt.w2list.resize(W2LIST_LENGTH)
+ evt.w2list[W2LIST_LOOKUP["PS"]] = global_generation_weight
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
diff --git a/km3buu/tests/test_output.py b/km3buu/tests/test_output.py
index 3ea01d93876bff4fdce04442240a278aac7c6d08..cb8de46664731163fc919035027aa4c3bb0cd3f4 100644
--- a/km3buu/tests/test_output.py
+++ b/km3buu/tests/test_output.py
@@ -76,9 +76,14 @@ class TestGiBUUOutput(unittest.TestCase):
def test_w2weights(self):
w2 = self.output.w2weights(123.0, 2.6e28, 4 * np.pi)
np.testing.assert_array_almost_equal(
- w2[:3], [7.64011311e+01, 3.61305080e-01, 1.13369700e+03],
+ w2[:3], [2.42100575e-06, 1.14490671e-08, 3.59246902e-05],
decimal=5)
+ def test_global_generation_weight(self):
+ self.assertAlmostEqual(self.output.global_generation_weight(4 * np.pi),
+ 2511.13458,
+ places=2)
+
@pytest.mark.skipif(not KM3NET_LIB_AVAILABLE,
reason="KM3NeT dataformat required")
@@ -131,4 +136,4 @@ class TestAANET(unittest.TestCase):
# CC/NC
np.testing.assert_equal(evt.w2list[10], 2)
# GiBUU weight
- np.testing.assert_almost_equal(evt.w2list[19], 0.8167222969153614)
+ np.testing.assert_almost_equal(evt.w2list[19], 0.004062418521597373)