diff --git a/km3buu/physics.py b/km3buu/physics.py
new file mode 100644
index 0000000000000000000000000000000000000000..b1e7e06aba4b24687482c52a13e04b009c1ab6f6
--- /dev/null
+++ b/km3buu/physics.py
@@ -0,0 +1,247 @@
+# Filename: physics.py
+"""
+Additional physics functionality
+
+"""
+
+__author__ = "Johannes Schumann"
+__copyright__ = "Copyright 2021, Johannes Schumann and the KM3NeT collaboration."
+__credits__ = []
+__license__ = "MIT"
+__maintainer__ = "Johannes Schumann"
+__email__ = "jschumann@km3net.de"
+__status__ = "Development"
+
+ELEC_PARAMS = {
+ "ELECa": 1.33356e5,
+ "ELECb": 1.66113e2,
+ "ELECc": 16.4949,
+ "ELECd": 1.5385e5,
+ "ELECe": 6.04871e5,
+}
+
+PION_PARAMS = {
+ "PIa": 0.538346,
+ "PIb": 1.32041,
+ "PIc": 0.737415,
+ "PId": -0.813861,
+ "PIe": -2.22444,
+ "PIf": -2.15795,
+ "PIg": -6.47242,
+ "PIh": -2.7567,
+ "PIx": 8.83498,
+}
+
+KAON_PARAMS = {
+ "Ka": 12.7537,
+ "Kb": -1.052,
+ "Kc": 3.49559,
+ "Kd": 16.7914,
+ "Ke": -0.355066,
+ "Kf": 2.77116,
+}
+
+KSHORT_PARAMS = {
+ "K0sa": 0.351242,
+ "K0sb": 0.613076,
+ "K0sc": 6.24682,
+ "K0sd": 2.8858,
+}
+
+KLONG_PARAMS = {
+ "K0la": 2.18152,
+ "K0lb": -0.632798,
+ "K0lc": 0.999514,
+ "K0ld": 1.36052,
+ "K0le": 4.22484,
+ "K0lf": -0.307859,
+}
+
+PROTON_PARAMS = {
+ "Pa": 12.1281,
+ "Pb": -17.1528,
+ "Pc": 0.573158,
+ "Pd": 34.1436,
+ "Pe": -0.28944,
+ "Pf": -13.2619,
+ "Pg": 24.1357,
+}
+
+NEUTRON_PARAMS = {
+ "Na": 1.24605,
+ "Nb": 0.63819,
+ "Nc": -0.802822,
+ "Nd": -0.935327,
+ "Ne": -6.1126,
+ "Nf": -1.96894,
+ "Ng": 0.716954,
+ "Nh": 2.68246,
+ "Ni": -9.39464,
+ "Nj": 15.0737,
+}
+
+HE_PARAMS = {
+ "Ma": 72.425,
+ "Mb": -49.417,
+ "Mc": 5.858,
+ "Md": 207.252,
+ "Me": 132.784,
+ "Mf": -10.277,
+ "Mg": -19.441,
+ "Mh": 58.598,
+ "Mi": 53.161,
+ "Mkref": 2.698,
+}
+
+
+def visible_energy_fraction(pdgid, energy):
+ """
+ Returns the visible energy fraction in the one particle approximation (OPA)
+
+ Parameters
+ ----------
+ pdgid: int
+ PDGID of the given particle
+ energy: float
+ Kinetic energy of the given particle
+ """
+
+ # Cover trivial cases, i.e. 'non-visible' particles and ultra-high energy
+ if pdgid in [11, -11, 22, 111, 221]:
+ return 1.0
+ elif energy > 1e7:
+ return high_energy_weight(energy)
+
+ tmp = energy if energy < 40. else 40.
+ weight = 0.0
+
+ if abs(pdgid) == 211:
+ weight = pion_weight(tmp)
+ elif pdgid == 130:
+ weight = klong_weight(tmp)
+ elif pdgid == 310:
+ weight = kshort_weight(tmp)
+ elif abs(pdgid) == 321:
+ weight = kaon_weight(tmp)
+ elif pdgid == 2112:
+ weight = neutron_weight(tmp)
+ elif pdgid == 2122:
+ weight = proton_weight(tmp)
+ elif pdgid in [12, -13, 13, 113, 213]:
+ weight = 0.0
+ else:
+ weight = proton_weight(tmp)
+
+ if energy < 40.:
+ return weight
+ else:
+ he_weight = high_energy_weight(energy)
+ he40GeV_weight = high_energy_weight(40.)
+ return he_weight - (he40GeV_weight -
+ weight) * (7. - np.log10(energy)) / 5.398
+
+
+def number_photons_per_electron(energy):
+ """Expected number of photons for electrons as function of energy"""
+ exp_coeff = np.exp(-energy / ELEC_PARAMS["ELECc"])
+ n = (ELEC_PARAMS["ELECa"] * energy + ELEC_PARAMS["ELECb"]) * exp_coeff + (
+ ELEC_PARAMS["ELECd"] * energy + ELEC_PARAMS["ELECe"]) * (1 - exp_coeff)
+ return n
+
+
+def pion_weight(energy):
+ norm = nphotons_per_electron(energy)
+ if energy < 6e-2:
+ return 1e4 * PION_PARAMS["PIa"] / norm
+ elif energy < 1.5e-1:
+ return 1e4 * PION_PARAMS["PIx"] * energy / norm
+ else:
+ logE = np.log(energy)
+ return (1e5 * PION_PARAMS["PIb"] * energy +
+ (energy**(1. - 1. / PION_PARAMS["PIc"])) *
+ (PION_PARAMS["PId"] * 1e4 + 1e4 * PION_PARAMS["PIe"] * logE +
+ 1e4 * PION_PARAMS["PIf"] * logE**2 + 1e3 * PION_PARAMS["PIg"]
+ * logE**3 + 1e2 * PION_PARAMS["PIh"] * logE**4)) / norm
+
+
+def kaon_weight(energy):
+ norm = number_photons_per_electron(energy)
+ if energy > 0.26:
+ exp_coeff = np.exp(-energy / KAON_PARAMS["Kc"])
+ return (1e4 * KAON_PARAMS["Ka"] * (energy + KAON_PARAMS["Kb"]) *
+ (1. - exp_coeff) + 1e4 *
+ (KAON_PARAMS["Kd"] * energy + KAON_PARAMS["Ke"]) *
+ exp_coeff) / norm
+ else:
+ return KAON_PARAMS["Kf"] * 1e4 / norm
+
+
+def kshort_weight(energy):
+ norm = number_photons_per_electron(energy)
+ return (KSHORT_PARAMS["K0sa"] * 1e5 + KSHORT_PARAMS["K0sb"] * 1e5 * energy
+ + energy * KSHORT_PARAMS["K0sc"] * 1e4 *
+ np.log(KSHORT_PARAMS["K0sd"] + 1. / energy)) / norm
+
+
+def klong_weight(energy):
+ norm = number_photons_per_electron(energy)
+ if energy < 1.5:
+ return (1e4 * KLONG_PARAMS["K0la"] +
+ energy * 1e5 * KLONG_PARAMS["K0lb"] * np.log(energy) +
+ 1e5 * KLONG_PARAMS["K0lc"] * energy**1.5) / norm
+ else:
+ return (energy * KLONG_PARAMS["K0ld"] * 1e5 +
+ energy**(1. - 1. / KLONG_PARAMS["K0le"]) *
+ KLONG_PARAMS["K0lf"] * 1e5) / norm
+
+
+def proton_weight(energy):
+ exp_coeff = np.exp(-energy / PROTON_PARAMS["Pc"])
+ norm = number_photons_per_electron(energy)
+ weight = 1e4 * (PROTON_PARAMS["Pa"] * energy + PROTON_PARAMS["Pb"]) * (
+ 1 - exp_coeff) + 1e4 * (PROTON_PARAMS["Pd"] * energy +
+ PROTON_PARAMS["Pe"] +
+ PROTON_PARAMS["Pf"] * energy**2 +
+ PROTON_PARAMS["Pg"] * energy**3) * exp_coeff
+ return weight / norm
+
+
+def neutron_weight(energy):
+ norm = number_photons_per_electron(energy)
+ if energy > 0.5:
+ logE = np.log(energy)
+ return (NEUTRON_PARAMS["Na"] * 1e5 * energy +
+ 1e3 * energy**(1. - 1. / NEUTRON_PARAMS["Nb"]) *
+ (100 * NEUTRON_PARAMS["Nc"] + 100 * NEUTRON_PARAMS["Nd"] * logE
+ + 10 * NEUTRON_PARAMS["Ne"] * logE**2 +
+ 11 * NEUTRON_PARAMS["Nf"] * logE**3)) / norm
+ else:
+ return (1e3 * NEUTRON_PARAMS["Ng"] + 1e4 * NEUTRON_PARAMS["Nh"] *
+ energy + 1e4 * NEUTRON_PARAMS["Ni"] * energy**2 +
+ 1e4 * NEUTRON_PARAMS["Nj"] * energy**3) / norm
+
+
+def high_energy_weight(energy):
+ """
+ High energy weight (valid above 40 GeV)
+
+ Parameters
+ ----------
+ energy: float
+ Kinetic energy of the given particle
+ """
+ if energy < 0.2:
+ return 0.292
+
+ logE = np.log10(energy)
+ mlc = (HE_PARAMS["Ma"] - HE_PARAMS["Mf"]) / HE_PARAMS["Mkref"]
+ denom = HE_PARAMS["Mi"] + HE_PARAMS["Mh"] * logE + HE_PARAMS[
+ "Mg"] * logE**2 + HE_PARAMS["Mf"] * logE**3 + mlc * logE**4
+ if denom <= 0:
+ return 0
+ num = HE_PARAMS[
+ "Me"] + HE_PARAMS["Md"] * logE + HE_PARAMS["Mc"] * logE**2 + HE_PARAMS[
+ "Mb"] * logE**3 + HE_PARAMS["Ma"] * logE**4 + mlc * logE**5
+ lognp = num / denom
+ Ee = 10.**(lognp - HE_PARAMS["Mkref"])
+ return Ee / energy