diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index fe594f3afa5c20989b0dfc58a32a94f986b5fe7c..26d6814423d8948e08241babec5e1f646f0ba06a 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -1,3 +1,3 @@
add_subdirectory (NullModel)
-add_subdirectory (Sibyll)
+#add_subdirectory (Sibyll)
add_subdirectory (StackInspector)
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index 041fdb8c34e8947d9ecc10b19961b7a2a653519b..cd745af06c6fcff22152df354eb72a7684c98958 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -1,4 +1,17 @@
-
+add_custom_command (
+ OUTPUT ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
+ COMMAND ${PROJECT_SOURCE_DIR}/Processes/Sibyll/code_generator.py
+ ${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl
+ ${PROJECT_SOURCE_DIR}/Processes/Sibyll/sibyll_codes.dat
+ DEPENDS code_generator.py
+ sibyll_codes.dat
+ ${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl
+ WORKING_DIRECTORY
+ ${PROJECT_BINARY_DIR}/Processes/Sibyll/Particles/
+ COMMENT "Generate conversion tables for particle codes SIBYLL <-> CORSIKA"
+ VERBATIM
+ )
+
set (
MODEL_SOURCES
ParticleConversion.cc
@@ -7,6 +20,7 @@ set (
set (
MODEL_HEADERS
ParticleConversion.h
+ ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
set (
@@ -48,7 +62,10 @@ install (
# --------------------
# code unit testing
-add_executable (testSibyll testSibyll.cc)
+add_executable (testSibyll
+ testSibyll.cc
+ ${MODEL_HEADERS}
+ )
target_link_libraries (
testSibyll
diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h
index f39c1297964d33846bedcb6160203f11f882c95d..e4461a24acacf3f86516a39f8787750730ceca34 100644
--- a/Processes/Sibyll/ParticleConversion.h
+++ b/Processes/Sibyll/ParticleConversion.h
@@ -5,195 +5,25 @@
#include <map>
-namespace corsika::process {
+namespace corsika::process::sibyll {
+ enum class Code : int8_t;
+ using PIDIntType = std::underlying_type<Code>::type;
- namespace sibyll {
+#include <corsika/processes/sibyll/Generated.inc>
- enum class PID : int {
- E_MINUS = 3,
- E_PLUS = 2,
- NU_E = 15,
- NU_E_BAR = 16,
- MU_MINUS = 5,
- MU_PLUS = 4,
- NU_MU = 17,
- NU_MU_BAR = 18,
- TAU_MINUS = 91,
- TAU_PLUS = 90,
- NU_TAU = 92,
- NU_TAU_BAR = 93,
- GAMMA = 1,
- PI_0 = 6,
- RHO_0 = 27,
- K_L_0 = 11,
- PI_PLUS = 7,
- PI_MINUS = 8,
- RHO_PLUS = 25,
- RHO_MINUS = 26,
- ETA = 23,
- OMEGA = 32,
- K_S_0 = 12,
- K_STAR_0 = 30,
- K_STAR_BAR_0 = 31,
- K_PLUS = 9,
- K_MINUS = 10,
- K_STAR_PLUS = 28,
- K_STAR_MINUS = 29,
- D_PLUS = 59,
- D_MINUS = 60,
- D_STAR_PLUS = 78,
- D_STAR_MINUS = 79,
- D_0 = 71,
- D_BAR_0 = 72,
- D_STAR_0 = 80,
- D_STAR_BAR_0 = 81,
- D_S_PLUS = 74,
- D_S_MINUS = 75,
- D_STAR_S_PLUS = 76,
- D_STAR_S_MINUS = 77,
- ETA_C = 73,
- N_0 = 14,
- N_BAR_0 = -14,
- DELTA_0 = 42,
- DELTA_BAR_0 = -42,
- P_PLUS = 13,
- P_BAR_MINUS = -13,
- DELTA_PLUS = 41,
- DELTA_BAR_MINUS = -41,
- DELTA_PLUS_PLUS = 40,
- DELTA_BAR_MINUS_MINUS = -40,
- SIGMA_MINUS = 36,
- SIGMA_BAR_PLUS = -36,
- LAMBDA_0 = 39,
- LAMBDA_BAR_0 = -39,
- SIGMA_0 = 35,
- SIGMA_BAR_0 = -35,
- SIGMA_PLUS = 34,
- SIGMA_BAR_MINUS = -34,
- XI_MINUS = 38,
- XI_BAR_PLUS = -38,
- XI_0 = 37,
- XI_BAR_0 = -37,
- OMEGA_MINUS = 49,
- OMEGA_BAR_PLUS = -49,
- SIGMA_C_0 = 86,
- SIGMA_C_BAR_0 = -86,
- SIGMA_STAR_C_0 = 96,
- SIGMA_STAR_C_BAR_0 = -96,
- LAMBDA_C_PLUS = 89,
- LAMBDA_C_BAR_MINUS = -89,
- XI_C_0 = 88,
- XI_C_BAR_0 = -88,
- SIGMA_C_PLUS = 85,
- SIGMA_C_BAR_MINUS = -85,
- SIGMA_STAR_C_PLUS = 95,
- SIGMA_STAR_C_BAR_MINUS = -95,
- SIGMA_C_PLUS_PLUS = 84,
- SIGMA_C_BAR_MINUS_MINUS = -84,
- SIGMA_STAR_C_PLUS_PLUS = 94,
- SIGMA_STAR_C_BAR_MINUS_MINUS = -94,
- XI_C_PLUS = 87,
- XI_C_BAR_MINUS = -87,
- OMEGA_C_0 = 99,
- OMEGA_C_BAR_0 = -99,
- J_PSI = 83,
- VOID = 0,
- };
+ bool handledBySibyll(corsika::particles::Code pCode) {
+ return handleable[static_cast<CodeIntType>(pCode)];
+ }
- static const std::map<sibyll::PID, corsika::particles::Code> Sibyll2Corsika = {
- {PID::E_MINUS, corsika::particles::Code::Electron},
- {PID::E_PLUS, corsika::particles::Code::Positron},
- {PID::NU_E, corsika::particles::Code::NuE},
- {PID::NU_E_BAR, corsika::particles::Code::NuEBar},
- {PID::MU_MINUS, corsika::particles::Code::MuMinus},
- {PID::MU_PLUS, corsika::particles::Code::MuPlus},
- {PID::NU_MU, corsika::particles::Code::NuMu},
- {PID::NU_MU_BAR, corsika::particles::Code::NuMuBar},
- {PID::TAU_MINUS, corsika::particles::Code::TauMinus},
- /*
- TAU_PLUS = 90,
- NU_TAU = 92,
- NU_TAU_BAR = 93,
- GAMMA = 1,
- PI_0 = 6,
- RHO_0 = 27,
- K_L_0 = 11,
- PI_PLUS = 7,
- PI_MINUS = 8,
- RHO_PLUS = 25,
- RHO_MINUS = 26,
- ETA = 23,
- OMEGA = 32,
- K_S_0 = 12,
- K_STAR_0 = 30,
- K_STAR_BAR_0 = 31,
- K_PLUS = 9,
- K_MINUS = 10,
- K_STAR_PLUS = 28,
- K_STAR_MINUS = 29,
- D_PLUS = 59,
- D_MINUS = 60,
- D_STAR_PLUS = 78,
- D_STAR_MINUS = 79,
- D_0 = 71,
- D_BAR_0 = 72,
- D_STAR_0 = 80,
- D_STAR_BAR_0 = 81,
- D_S_PLUS = 74,
- D_S_MINUS = 75,
- D_STAR_S_PLUS = 76,
- D_STAR_S_MINUS = 77,
- ETA_C = 73,
- N_0 = 14,
- N_BAR_0 = -14,
- DELTA_0 = 42,
- DELTA_BAR_0 = -42,
- P_PLUS = 13,
- P_BAR_MINUS = -13,
- DELTA_PLUS = 41,
- DELTA_BAR_MINUS = -41,
- DELTA_PLUS_PLUS = 40,
- DELTA_BAR_MINUS_MINUS = -40,
- SIGMA_MINUS = 36,
- SIGMA_BAR_PLUS = -36,
- LAMBDA_0 = 39,
- LAMBDA_BAR_0 = -39,
- SIGMA_0 = 35,
- SIGMA_BAR_0 = -35,
- SIGMA_PLUS = 34,
- SIGMA_BAR_MINUS = -34,
- XI_MINUS = 38,
- XI_BAR_PLUS = -38,
- XI_0 = 37,
- XI_BAR_0 = -37,
- OMEGA_MINUS = 49,
- OMEGA_BAR_PLUS = -49,
- SIGMA_C_0 = 86,
- SIGMA_C_BAR_0 = -86,
- SIGMA_STAR_C_0 = 96,
- SIGMA_STAR_C_BAR_0 = -96,
- LAMBDA_C_PLUS = 89,
- LAMBDA_C_BAR_MINUS = -89,
- XI_C_0 = 88,
- XI_C_BAR_0 = -88,
- SIGMA_C_PLUS = 85,
- SIGMA_C_BAR_MINUS = -85,
- SIGMA_STAR_C_PLUS = 95,
- SIGMA_STAR_C_BAR_MINUS = -95,
- SIGMA_C_PLUS_PLUS = 84,
- SIGMA_C_BAR_MINUS_MINUS = -84,
- SIGMA_STAR_C_PLUS_PLUS = 94,
- SIGMA_STAR_C_BAR_MINUS_MINUS = -94,
- XI_C_PLUS = 87,
- XI_C_BAR_MINUS = -87,
- OMEGA_C_0 = 99,
- OMEGA_C_BAR_0 = -99,
- J_PSI = 83,
- VOID = 0,*/
- };
+ Code constexpr ConvertToSibyll(corsika::particles::Code pCode) {
+ //~ assert(handledBySibyll(pCode));
+ return static_cast<Code>(corsika2sibyll[static_cast<CodeIntType>(pCode)]);
+ }
- } // namespace sibyll
+ corsika::particles::Code constexpr ConvertFromSibyll(Code pCode) {
+ return sibyll2corsika[static_cast<PIDIntType>(pCode) - minSibyll];
+ }
-} // namespace corsika::process
+} // namespace corsika::process::sibyll
#endif
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
new file mode 100755
index 0000000000000000000000000000000000000000..a51deef031daabd25c3f39cf5841eb730fec6b37
--- /dev/null
+++ b/Processes/Sibyll/code_generator.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+
+# (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+#
+# See file AUTHORS for a list of contributors.
+#
+# This software is distributed under the terms of the GNU General Public
+# Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+# the license.
+
+
+import pickle, sys, itertools
+
+
+# loads the pickled pythia_db (which is an OrderedDict)
+def load_pythiadb(filename):
+ with open(filename, "rb") as f:
+ pythia_db = pickle.load(f)
+ return pythia_db
+
+
+#
+def read_sibyll_codes(filename, pythia_db):
+ with open(filename) as f:
+ for line in f:
+ line = line.strip()
+ if line[0] == '#':
+ continue
+ identifier, sib_code = line.split()
+ try:
+ pythia_db[identifier]["sibyll_code"] = int(sib_code)
+ except KeyError as e:
+ raise Exception("Identifier '{:s}' not found in pythia_db".format(identifier))
+
+
+# generates the enum to access sibyll particles by readable names
+def generate_sibyll_enum(pythia_db):
+ output = "enum class PID : int16_t {\n"
+ for identifier, d in pythia_db.items():
+ if d.get('sibyll_code') != None:
+ output += " {:s} = {:d},\n".format(identifier, d['sibyll_code'])
+ output += "};\n"
+ return output
+
+
+# generates the look-up table to convert corsika codes to sibyll codes
+def generate_corsika2sibyll(pythia_db):
+ string = "std::array<PIDIntType, {:d}> constexpr corsika2sibyll = {{".format(len(pythia_db))
+ for identifier, d in pythia_db.items():
+ sibCode = d.get("sibyll_code", 0)
+ string += " {:d}, // {:s}\n".format(sibCode, identifier if sibCode else identifier + " (not implemented in SIBYLL)")
+ string += "};\n"
+ return string
+
+
+# generates the look-up table to convert sibyll codes to corsika codes
+def generate_sibyll2corsika(pythia_db):
+ d = {}
+ for identifier, p in pythia_db.items():
+ if 'sibyll_code' in p:
+ sib_code = p['sibyll_code']
+ corsika_code = p['ngc_code']
+ d[sib_code] = (corsika_code, identifier)
+
+ string = "std::array<corsika::particles::CodeIntType, {:d}> sibyll2corsika = {{\n".format(len(d))
+
+ for k in range(min(d.keys()), max(d.keys())+1):
+ if k in d:
+ corsika_code = d[k][0]
+ identifier = d[k][1]
+ else:
+ corsika_code = 0
+ identifier = ""
+ string += " {:d}, // {:s}\n".format(corsika_code, identifier)
+
+ string += "};\n"
+ string += "PIDIntType constexpr minSibyll = {:d};\n".format(min(d.keys()))
+ return string
+
+
+# generates the bitset for the flag whether Sibyll knows the particle
+def generate_handles_particle(pythia_db):
+ num_particles = len(pythia_db)
+ num_bytes = num_particles // 32 + 1
+ string = "Bitset2::bitset2<{:d}> constexpr handleable{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes)
+
+ numeric = 0
+ for identifier, d in reversed(pythia_db.items()):
+ handledBySibyll = int("sibyll_code" in d) & 0x1
+ numeric = (numeric << 1) | handledBySibyll
+
+ while numeric != 0:
+ low = numeric & 0xFFFFFFFF
+ numeric = numeric >> 32
+ string += " 0x{:0x},\n".format(low)
+
+ string += "}}};\n"
+ return string
+
+
+if __name__ == "__main__":
+ if len(sys.argv) != 3:
+ print("usage: {:s} <pythia_db.pkl> <sibyll_codes.dat>".format(sys.argv[0]), file=sys.stderr)
+ sys.exit(1)
+
+ print("code_generator.py for SIBYLL")
+
+ pythia_db = load_pythiadb(sys.argv[1])
+
+ read_sibyll_codes(sys.argv[2], pythia_db)
+
+ with open("Generated.inc", "w") as f:
+ print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
+ print(generate_sibyll_enum(pythia_db), file=f)
+ print(generate_corsika2sibyll(pythia_db), file=f)
+ print(generate_handles_particle(pythia_db), file=f)
+ print(generate_sibyll2corsika(pythia_db), file=f)
diff --git a/Processes/Sibyll/sibyll_codes.dat b/Processes/Sibyll/sibyll_codes.dat
index edc282f17398c435f5a3b0c81d6a51daf0c027c6..880699309325533fecbb32a377982e26a2954585 100644
--- a/Processes/Sibyll/sibyll_codes.dat
+++ b/Processes/Sibyll/sibyll_codes.dat
@@ -1,88 +1,88 @@
-E_MINUS 3
-E_PLUS 2
-NU_E 15
-NU_E_BAR 16
-MU_MINUS 5
-MU_PLUS 4
-NU_MU 17
-NU_MU_BAR 18
-TAU_MINUS 91
-TAU_PLUS 90
-NU_TAU 92
-NU_TAU_BAR 93
-GAMMA 1
-PI_0 6
-RHO_0 27
-K_L_0 11
-PI_PLUS 7
-PI_MINUS 8
-RHO_PLUS 25
-RHO_MINUS 26
-ETA 23
-OMEGA 32
-K_S_0 12
-K_STAR_0 30
-K_STAR_BAR_0 31
-K_PLUS 9
-K_MINUS 10
-K_STAR_PLUS 28
-K_STAR_MINUS 29
-D_PLUS 59
-D_MINUS 60
-D_STAR_PLUS 78
-D_STAR_MINUS 79
-D_0 71
-D_BAR_0 72
-D_STAR_0 80
-D_STAR_BAR_0 81
-D_S_PLUS 74
-D_S_MINUS 75
-D_STAR_S_PLUS 76
-D_STAR_S_MINUS 77
-ETA_C 73
-N_0 14
-N_BAR_0 -14
-DELTA_0 42
-DELTA_BAR_0 -42
-P_PLUS 13
-P_BAR_MINUS -13
-DELTA_PLUS 41
-DELTA_BAR_MINUS -41
-DELTA_PLUS_PLUS 40
-DELTA_BAR_MINUS_MINUS -40
-SIGMA_MINUS 36
-SIGMA_BAR_PLUS -36
-LAMBDA_0 39
-LAMBDA_BAR_0 -39
-SIGMA_0 35
-SIGMA_BAR_0 -35
-SIGMA_PLUS 34
-SIGMA_BAR_MINUS -34
-XI_MINUS 38
-XI_BAR_PLUS -38
-XI_0 37
-XI_BAR_0 -37
-OMEGA_MINUS 49
-OMEGA_BAR_PLUS -49
-SIGMA_C_0 86
-SIGMA_C_BAR_0 -86
-SIGMA_STAR_C_0 96
-SIGMA_STAR_C_BAR_0 -96
-LAMBDA_C_PLUS 89
-LAMBDA_C_BAR_MINUS -89
-XI_C_0 88
-XI_C_BAR_0 -88
-SIGMA_C_PLUS 85
-SIGMA_C_BAR_MINUS -85
-SIGMA_STAR_C_PLUS 95
-SIGMA_STAR_C_BAR_MINUS -95
-SIGMA_C_PLUS_PLUS 84
-SIGMA_C_BAR_MINUS_MINUS -84
-SIGMA_STAR_C_PLUS_PLUS 94
-SIGMA_STAR_C_BAR_MINUS_MINUS -94
-XI_C_PLUS 87
-XI_C_BAR_MINUS -87
-OMEGA_C_0 99
-OMEGA_C_BAR_0 -99
-J_PSI 83
-VOID 0
+Electron 3
+Positron 2
+NuE 15
+NuEBar 16
+MuMinus 5
+MuPlus 4
+NuMu 17
+NuMuBar 18
+TauMinus 91
+TauPlus 90
+NuTau 92
+NuTauBar 93
+Gamma 1
+Pi0 6
+Rho0 27
+K0Long 11
+PiPlus 7
+PiMinus 8
+RhoPlus 25
+RhoMinus 26
+Eta 23
+Omega 32
+K0Short 12
+KStar0 30
+KStar0Bar 31
+KPlus 9
+KMinus 10
+KStarPlus 28
+KStarMinus 29
+DPlus 59
+DMinus 60
+DStarPlus 78
+DStarMinus 79
+D0 71
+D0Bar 72
+DStar0 80
+DStar0Bar 81
+DsPlus 74
+DsMinus 75
+DStarSPlus 76
+DStarSMinus 77
+EtaC 73
+Neutron 14
+AntiNeutron -14
+Delta0 42
+Delta0Bar -42
+Proton 13
+AntiProton -13
+DeltaPlus 41
+DeltaMinusBar -41
+DeltaPlusPlus 40
+DeltaMinusMinusBar -40
+SigmaMinus 36
+SigmaPlusBar -36
+Lambda0 39
+Lambda0Bar -39
+Sigma0 35
+Sigma0Bar -35
+SigmaPlus 34
+SigmaMinusBar -34
+XiMinus 38
+XiPlusBar -38
+Xi0 37
+Xi0Bar -37
+OmegaMinus 49
+OmegaPlusBar -49
+SigmaC0 86
+SigmaC0Bar -86
+SigmaStarC0 96
+SigmaStarC0Bar -96
+LambdaCPlus 89
+LambdaCMinusBar -89
+XiC0 88
+XiC0Bar -88
+SigmaCPlus 85
+SigmaCMinusBar -85
+SigmaStarCPlus 95
+SigmaStarCMinusBar -95
+SigmaCPlusPlus 84
+SigmaCMinusMinusBar -84
+SigmaStarCPlusPlus 94
+SigmaStarCMinusMinusBar -94
+XiCPlus 87
+XiCMinusBar -87
+OmegaC0 99
+OmegaC0Bar -99
+Jpsi 83
+#Unknown 0
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 4c5bed4df1b42c351587569ed5a0b192f19f1f2d..581e5a10b0e66e2c111a3958f51b23548e9d5197 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -11,17 +11,20 @@ using namespace corsika;
TEST_CASE("Sibyll", "[processes]") {
- SECTION("ParticleConversion") {
+ SECTION("Sibyll -> Corsika") {
REQUIRE(corsika::particles::Electron::GetCode() ==
- process::sibyll::Sibyll2Corsika.at(process::sibyll::PID::E_MINUS));
+ process::sibyll::ConvertFromSibyll(process::sibyll::Code::Electron));
}
- SECTION("Data") {
- REQUIRE(corsika::particles::GetName(process::sibyll::Sibyll2Corsika.at(
- process::sibyll::PID::E_PLUS)) == "e+");
+ SECTION("Corsika -> Sibyll") {
+ REQUIRE(process::sibyll::ConvertToSibyll(corsika::particles::Electron::GetCode()) ==
+ process::sibyll::Code::Electron);
}
- SECTION("bla") {}
+ SECTION("handledBySibyll") {
+ REQUIRE(process::sibyll::handledBySibyll(corsika::particles::Electron::GetCode()));
- SECTION("blubb") {}
+ REQUIRE_FALSE(
+ process::sibyll::handledBySibyll(corsika::particles::XiPrimeC0::GetCode()));
+ }
}