IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 829e0de6 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

conversion for Sibyll

parent e4c039cf
No related branches found
No related tags found
No related merge requests found
add_subdirectory (NullModel)
add_subdirectory (Sibyll)
#add_subdirectory (Sibyll)
add_subdirectory (StackInspector)
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
......
......@@ -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
#!/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)
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
......@@ -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()));
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment