IAP GITLAB

Skip to content
Snippets Groups Projects
Commit c2cef4eb authored by Ralf Ulrich's avatar Ralf Ulrich
Browse files

Merge branch '35-handle-nuclei' into 'master' 

Resolve "Handle nuclei"

Closes #35

See merge request AirShowerPhysics/corsika!33
parents edbf26d5 989f9ef2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/Particles/CMakeLists.txt
+ 3
− 1
View file @ c2cef4eb
add_custom_command ( add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc OUTPUT ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
COMMAND ${PROJECT_SOURCE_DIR}/Framework/Particles/pdxml_reader.py COMMAND ${PROJECT_SOURCE_DIR}/Framework/Particles/pdxml_reader.py
${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleData.xml ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleData.xml
${PROJECT_SOURCE_DIR}/Framework/Particles/NuclearData.xml
${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleClassNames.xml ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleClassNames.xml
DEPENDS pdxml_reader.py DEPENDS pdxml_reader.py
ParticleData.xml ParticleData.xml
NuclearData.xml
ParticleClassNames.xml ParticleClassNames.xml
WORKING_DIRECTORY WORKING_DIRECTORY
${PROJECT_BINARY_DIR}/Framework/Particles/ ${PROJECT_BINARY_DIR}/Framework/Particles/
......
Framework/Particles/NuclearData.xml 0 → 100644
+ 27
− 0
View file @ c2cef4eb
<chapter name="Nuclear Data">
<particle id="1000010010" name="hydrogen" A="1" Z="1" >
</particle>
<particle id="1000010020" name="deuterium" A="2" Z="1" >
</particle>
<particle id="1000010030" name="tritium" A="3" Z="1" >
</particle>
<particle id="1000020040" name="helium" A="4" Z="2" >
</particle>
<particle id="1000060120" name="carbon" A="12" Z="6" >
</particle>
<particle id="1000070140" name="nitrogen" A="14" Z="7" >
</particle>
<particle id="1000080160" name="oxygen" A="16" Z="8" >
</particle>
<particle id="1000180400" name="argon" A="40" Z="18" >
</particle>
</chapter>
Framework/Particles/ParticleProperties.h
+ 23
− 6
View file @ c2cef4eb
...@@ -40,7 +40,7 @@ namespace corsika::particles { ...@@ -40,7 +40,7 @@ namespace corsika::particles {
enum class Code : int16_t; enum class Code : int16_t;
using PDGCodeType = int16_t; using PDGCodeType = int32_t;
using CodeIntType = std::underlying_type<Code>::type; using CodeIntType = std::underlying_type<Code>::type;
// forward declarations to be used in GeneratedParticleProperties // forward declarations to be used in GeneratedParticleProperties
...@@ -51,24 +51,28 @@ namespace corsika::particles { ...@@ -51,24 +51,28 @@ namespace corsika::particles {
constexpr std::string const& GetName(Code const); constexpr std::string const& GetName(Code const);
corsika::units::si::TimeType constexpr GetLifetime(Code const); corsika::units::si::TimeType constexpr GetLifetime(Code const);
bool constexpr IsNucleus(Code const);
int constexpr GetNucleusA(Code const);
int constexpr GetNucleusZ(Code const);
#include <corsika/particles/GeneratedParticleProperties.inc> #include <corsika/particles/GeneratedParticleProperties.inc>
/*! /*!
* returns mass of particle * returns mass of particle
*/ */
corsika::units::hep::MassType constexpr GetMass(Code const p) { corsika::units::hep::MassType constexpr GetMass(Code const p) {
return masses[static_cast<CodeIntType const>(p)]; return detail::masses[static_cast<CodeIntType const>(p)];
} }
PDGCodeType constexpr GetPDG(Code const p) { PDGCodeType constexpr GetPDG(Code const p) {
return pdg_codes[static_cast<CodeIntType const>(p)]; return detail::pdg_codes[static_cast<CodeIntType const>(p)];
} }
/*! /*!
* returns electric charge of particle / (e/3). * returns electric charge of particle / (e/3).
*/ */
int16_t constexpr GetElectricChargeNumber(Code const p) { int16_t constexpr GetElectricChargeNumber(Code const p) {
return electric_charges[static_cast<CodeIntType const>(p)]; return detail::electric_charges[static_cast<CodeIntType const>(p)];
} }
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) { corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) {
...@@ -76,13 +80,26 @@ namespace corsika::particles { ...@@ -76,13 +80,26 @@ namespace corsika::particles {
} }
constexpr std::string const& GetName(Code const p) { constexpr std::string const& GetName(Code const p) {
return names[static_cast<CodeIntType const>(p)]; return detail::names[static_cast<CodeIntType const>(p)];
} }
corsika::units::si::TimeType constexpr GetLifetime(Code const p) { corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
return lifetime[static_cast<CodeIntType const>(p)]; return detail::lifetime[static_cast<CodeIntType const>(p)] * corsika::units::si::second;
} }
bool constexpr IsNucleus(Code const p) {
return detail::isNucleus[static_cast<CodeIntType const>(p)];
}
int constexpr GetNucleusA(Code const p) {
return detail::nucleusA[static_cast<CodeIntType const>(p)];
}
int constexpr GetNucleusZ(Code const p) {
return detail::nucleusZ[static_cast<CodeIntType const>(p)];
}
namespace io { namespace io {
std::ostream& operator<<(std::ostream& stream, Code const p); std::ostream& operator<<(std::ostream& stream, Code const p);
......
Framework/Particles/pdxml_reader.py
+ 176
− 65
View file @ c2cef4eb
...@@ -5,21 +5,23 @@ import xml.etree.ElementTree as ET ...@@ -5,21 +5,23 @@ import xml.etree.ElementTree as ET
from collections import OrderedDict from collections import OrderedDict
import pickle import pickle
GeVfm = 0.19732696312541853
c_speed_of_light = 29.9792458e10 # mm / s
# for nuclear masses
mneutron = 0.9395654133 # GeV
mproton = 0.9382720813 # GeV
############################################################## ##############################################################
# #
# reading xml input data, return line by line particle data # reading xml input data, return line by line particle data
# #
def parsePythia(filename):
def parse(filename):
tree = ET.parse(filename) tree = ET.parse(filename)
root = tree.getroot() root = tree.getroot()
GeVfm = 0.19732696312541853
c_speed_of_light = 29.9792458e10 # mm / s
for particle in root.iter("particle"): for particle in root.iter("particle"):
name = particle.attrib["name"] name = particle.attrib["name"]
antiName = "Unknown" antiName = "Unknown"
if ("antiName" in particle.attrib): if ("antiName" in particle.attrib):
antiName = particle.attrib["antiName"] antiName = particle.attrib["antiName"]
...@@ -47,12 +49,39 @@ def parse(filename): ...@@ -47,12 +49,39 @@ def parse(filename):
yield (-pdg_id, antiName, mass, -electric_charge, name, ctau/c_speed_of_light) yield (-pdg_id, antiName, mass, -electric_charge, name, ctau/c_speed_of_light)
############################################################## ##############################################################
# #
# returns dict with particle codes and class names # reading xml input data, return line by line particle data
# #
def parseNuclei(filename):
tree = ET.parse(filename)
root = tree.getroot()
for particle in root.iter("particle"):
name = particle.attrib["name"]
antiName = "Unknown"
if ("antiName" in particle.attrib):
antiName = particle.attrib["antiName"]
pdg_id = int(particle.attrib["id"])
A = int(particle.attrib["A"])
Z = int(particle.attrib["Z"])
# mass in GeV
if ("mass" in particle.attrib):
mass = particle.attrib["mass"]
else:
mass = (A-Z)*mneutron + Z*mproton
electric_charge = Z*3 # in units of e/3
ctau = float('Inf')
yield (pdg_id, name, mass, electric_charge, antiName, ctau/c_speed_of_light, A, Z)
##############################################################
#
# returns dict with particle codes and class names
#
def class_names(filename): def class_names(filename):
tree = ET.parse(filename) tree = ET.parse(filename)
root = tree.getroot() root = tree.getroot()
...@@ -164,14 +193,12 @@ def c_identifier_camel(name): ...@@ -164,14 +193,12 @@ def c_identifier_camel(name):
########################################################## ##########################################################
# #
# returns dict containing all data from pythia-xml input # returns dict containing all data from pythia-xml input
# #
def read_pythia_db(filename, particle_db, classnames):
def build_pythia_db(filename, classnames): counter = itertools.count(len(particle_db))
particle_db = OrderedDict()
counter = itertools.count(0) for (pdg, name, mass, electric_charge, antiName, lifetime) in parsePythia(filename):
for (pdg, name, mass, electric_charge, antiName, lifetime) in parse(filename):
c_id = "Unknown" c_id = "Unknown"
if pdg in classnames: if pdg in classnames:
...@@ -186,7 +213,41 @@ def build_pythia_db(filename, classnames): ...@@ -186,7 +213,41 @@ def build_pythia_db(filename, classnames):
"mass" : mass, # in GeV "mass" : mass, # in GeV
"electric_charge" : electric_charge, # in e/3 "electric_charge" : electric_charge, # in e/3
"lifetime" : lifetime, "lifetime" : lifetime,
"ngc_code" : next(counter) "ngc_code" : next(counter),
"isNucleus" : False,
}
return particle_db
##########################################################
#
# returns dict containing all data from pythia-xml input
#
def read_nuclei_db(filename, particle_db, classnames):
counter = itertools.count(len(particle_db))
for (pdg, name, mass, electric_charge, antiName, lifetime, A, Z) in parseNuclei(filename):
c_id = "Unknown"
if pdg in classnames:
c_id = classnames[pdg]
else:
c_id = c_identifier_camel(name)
particle_db[c_id] = {
"name" : name,
"antiName" : antiName,
"pdg" : pdg,
"mass" : mass, # in GeV
"electric_charge" : electric_charge, # in e/3
"lifetime" : lifetime,
"ngc_code" : next(counter),
"A" : A,
"Z" : Z,
"isNucleus" : True,
} }
return particle_db return particle_db
...@@ -198,14 +259,13 @@ def build_pythia_db(filename, classnames): ...@@ -198,14 +259,13 @@ def build_pythia_db(filename, classnames):
# #
# return string with enum of all internal particle codes # return string with enum of all internal particle codes
# #
def gen_internal_enum(particle_db):
def gen_internal_enum(pythia_db):
string = ("enum class Code : int16_t {\n" string = ("enum class Code : int16_t {\n"
" FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops... " FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops...
for k in filter(lambda k: "ngc_code" in pythia_db[k], pythia_db): for k in filter(lambda k: "ngc_code" in particle_db[k], particle_db):
last_ngc_id = pythia_db[k]['ngc_code'] last_ngc_id = particle_db[k]['ngc_code']
string += " {key:s} = {code:d},\n".format(key = k, code = last_ngc_id) string += " {key:s} = {code:d},\n".format(key = k, code = last_ngc_id)
string += (" LastParticle = {:d},\n" # identifier for eventual loops... string += (" LastParticle = {:d},\n" # identifier for eventual loops...
...@@ -223,52 +283,83 @@ def gen_internal_enum(pythia_db): ...@@ -223,52 +283,83 @@ def gen_internal_enum(pythia_db):
# #
# return string with all data arrays # return string with all data arrays
# #
def gen_properties(particle_db):
def gen_properties(pythia_db):
# number of particles, size of tables # number of particles, size of tables
string = "static constexpr std::size_t size = {size:d};\n".format(size = len(pythia_db)) string = "static constexpr std::size_t size = {size:d};\n".format(size = len(particle_db))
string += "\n" string += "\n"
# particle masses table # particle masses table
string += "static constexpr std::array<corsika::units::hep::MassType const, size> masses = {\n" string += "static constexpr std::array<corsika::units::hep::MassType const, size> masses = {\n"
for p in pythia_db.values(): for p in particle_db.values():
string += " {mass:f} * (1e9 * corsika::units::si::constants::eV), // {name:s}\n".format(mass = p['mass'], name = p['name']) string += " {mass:e} * (1e9 * corsika::units::si::constants::eV), // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n" string += "};\n\n"
# PDG code table # PDG code table
string += "static constexpr std::array<PDGCodeType const, size> pdg_codes = {\n" string += "static constexpr std::array<PDGCodeType const, size> pdg_codes = {\n"
for p in pythia_db.values(): for p in particle_db.values():
string += " {pdg:d}, // {name:s}\n".format(pdg = p['pdg'], name = p['name']) string += " {pdg:d}, // {name:s}\n".format(pdg = p['pdg'], name = p['name'])
string += "};\n" string += "};\n"
# name string table # name string table
string += "static const std::array<std::string const, size> names = {\n" string += "static const std::array<std::string const, size> names = {\n"
for p in pythia_db.values(): for p in particle_db.values():
string += " \"{name:s}\",\n".format(name = p['name']) string += " \"{name:s}\",\n".format(name = p['name'])
string += "};\n" string += "};\n"
# electric charges table # electric charges table
string += "static constexpr std::array<int16_t, size> electric_charges = {\n" string += "static constexpr std::array<int16_t, size> electric_charges = {\n"
for p in pythia_db.values(): for p in particle_db.values():
string += " {charge:d},\n".format(charge = p['electric_charge']) string += " {charge:d},\n".format(charge = p['electric_charge'])
string += "};\n" string += "};\n"
# anti-particle table # anti-particle table
# string += "static constexpr std::array<size, size> anti_particle = {\n" # string += "static constexpr std::array<size, size> anti_particle = {\n"
# for p in pythia_db.values(): # for p in particle_db.values():
# string += " {anti:d},\n".format(charge = p['anti_particle']) # string += " {anti:d},\n".format(charge = p['anti_particle'])
# string += "};\n" # string += "};\n"
# lifetime # lifetime
string += "static constexpr std::array<corsika::units::si::TimeType const, size> lifetime = {\n" #string += "static constexpr std::array<corsika::units::si::TimeType const, size> lifetime = {\n"
for p in pythia_db.values(): string += "static constexpr std::array<double const, size> lifetime = {\n"
for p in particle_db.values():
if p['lifetime'] == float("Inf") : if p['lifetime'] == float("Inf") :
string += " std::numeric_limits<double>::infinity() * corsika::units::si::second, \n" string += " std::numeric_limits<double>::infinity(), \n" # * corsika::units::si::second, \n"
else : else :
string += " {tau:f} * corsika::units::si::second, \n".format(tau = p['lifetime']) string += " {tau:e}, \n".format(tau = p['lifetime'])
#string += " {tau:e} * corsika::units::si::second, \n".format(tau = p['lifetime'])
string += "};\n"
### nuclear data ###
# is nucleus flag
string += "static constexpr std::array<bool, size> isNucleus = {\n"
for p in particle_db.values():
value = 'false'
if p['isNucleus']:
value = 'true'
string += " {val},\n".format(val = value)
string += "};\n"
# nucleus mass number A
string += "static constexpr std::array<int, size> nucleusA = {\n"
for p in particle_db.values():
A = 0
if p['isNucleus']:
A = p['A']
string += " {val},\n".format(val = A)
string += "};\n" string += "};\n"
# nucleus charge number Z
string += "static constexpr std::array<int, size> nucleusZ = {\n"
for p in particle_db.values():
Z = 0
if p['isNucleus']:
Z = p['Z']
string += " {val},\n".format(val = Z)
string += "};\n"
return string return string
...@@ -278,27 +369,29 @@ def gen_properties(pythia_db): ...@@ -278,27 +369,29 @@ def gen_properties(pythia_db):
# #
# return string with a list of classes for all particles # return string with a list of classes for all particles
# #
def gen_classes(particle_db):
def gen_classes(pythia_db):
string = "// list of C++ classes to access particle properties" string = "// list of C++ classes to access particle properties"
for cname in pythia_db: for cname in particle_db:
antiP = 'Unknown' antiP = 'Unknown'
for cname_anti in pythia_db: for cname_anti in particle_db:
if (pythia_db[cname_anti]['name'] == pythia_db[cname]['antiName']): if (particle_db[cname_anti]['name'] == particle_db[cname]['antiName']):
antiP = cname_anti antiP = cname_anti
break break
string += "\n"; string += "\n";
string += "/** @class " + cname + "\n\n" string += "/** @class " + cname + "\n\n"
string += " * Particle properties are taken from the PYTHIA8 ParticleData.xml file:<br>\n" string += " * Particle properties are taken from the PYTHIA8 ParticleData.xml file:<br>\n"
string += " * - pdg=" + str(pythia_db[cname]['pdg']) +"\n" string += " * - pdg=" + str(particle_db[cname]['pdg']) +"\n"
string += " * - mass=" + str(pythia_db[cname]['mass']) + " GeV \n" string += " * - mass=" + str(particle_db[cname]['mass']) + " GeV \n"
string += " * - charge= " + str(pythia_db[cname]['electric_charge']/3) + " \n" string += " * - charge= " + str(particle_db[cname]['electric_charge']/3) + " \n"
string += " * - name=" + str(cname) + "\n" string += " * - name=" + str(cname) + "\n"
string += " * - anti=" + str(antiP) + "\n" string += " * - anti=" + str(antiP) + "\n"
if (particle_db[cname]['isNucleus']):
string += " * - nuclear A=" + str(particle_db[cname]['A']) + "\n"
string += " * - nuclear Z=" + str(particle_db[cname]['Z']) + "\n"
string += "*/\n\n" string += "*/\n\n"
string += "class " + cname + " {\n" string += "class " + cname + " {\n"
string += " public:\n" string += " public:\n"
...@@ -308,6 +401,9 @@ def gen_classes(pythia_db): ...@@ -308,6 +401,9 @@ def gen_classes(pythia_db):
string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetElectricChargeNumber(Type); }\n" string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetElectricChargeNumber(Type); }\n"
string += " static std::string const& GetName() { return corsika::particles::GetName(Type); }\n" string += " static std::string const& GetName() { return corsika::particles::GetName(Type); }\n"
string += " static constexpr Code GetAntiParticle() { return AntiType; }\n" string += " static constexpr Code GetAntiParticle() { return AntiType; }\n"
string += " static constexpr bool IsNucleus() { return corsika::particles::IsNucleus(Type); }\n"
string += " static constexpr int GetNucleusA() { return corsika::particles::GetNucleusA(Type); }\n"
string += " static constexpr int GetNucleusZ() { return corsika::particles::GetNucleusZ(Type); }\n"
string += " static constexpr Code Type = Code::" + cname + ";\n" string += " static constexpr Code Type = Code::" + cname + ";\n"
string += " static constexpr Code AntiType = Code::" + antiP + ";\n" string += " static constexpr Code AntiType = Code::" + antiP + ";\n"
string += " private:\n" string += " private:\n"
...@@ -320,17 +416,30 @@ def gen_classes(pythia_db): ...@@ -320,17 +416,30 @@ def gen_classes(pythia_db):
############################################################### ###############################################################
# #
# #
def inc_start(): def inc_start():
string = ('// generated by pdxml_reader.py\n' string = ('// generated by pdxml_reader.py\n'
'// MANUAL EDITS ON OWN RISK\n') '// MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN. \n')
return string return string
############################################################### ###############################################################
# #
# #
def detail_start():
string = ('namespace detail {\n\n')
return string
###############################################################
#
#
def detail_end():
string = "\n}//end namespace detail\n"
return string
###############################################################
#
#
def inc_end(): def inc_end():
string = "" string = ""
return string return string
...@@ -338,36 +447,38 @@ def inc_end(): ...@@ -338,36 +447,38 @@ def inc_end():
################################################################### ###################################################################
# #
# Serialize pythia_db into file # Serialize particle_db into file
# #
def serialize_particle_db(particle_db, file):
pickle.dump(particle_db, file)
def serialize_pythia_db(pythia_db, file):
pickle.dump(pythia_db, file)
################################################################### ###################################################################
# #
# Main function # Main function
# #
if __name__ == "__main__": if __name__ == "__main__":
if len(sys.argv) != 3: if len(sys.argv) != 4:
print("usage: {:s} <Pythia8.xml> <ClassNames.xml>".format(sys.argv[0]), file=sys.stderr) print("usage: {:s} <Pythia8.xml> <Nuclei.xml> <ClassNames.xml>".format(sys.argv[0]), file=sys.stderr)
sys.exit(1) sys.exit(1)
print("\n pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n") print("\n pdxml_reader.py: Automatically produce particle-properties from input files\n")
names = class_names(sys.argv[2]) names = class_names(sys.argv[3])
pythia_db = build_pythia_db(sys.argv[1], names) particle_db = OrderedDict()
read_pythia_db(sys.argv[1], particle_db, names)
print (str(pythia_db)) read_nuclei_db(sys.argv[2], particle_db, names)
with open("GeneratedParticleProperties.inc", "w") as f: with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f) print(inc_start(), file=f)
print(gen_internal_enum(pythia_db), file=f) print(gen_internal_enum(particle_db), file=f)
print(gen_properties(pythia_db), file=f) print(detail_start(), file=f)
print(gen_classes(pythia_db), file=f) print(gen_properties(particle_db), file=f)
print(detail_end(), file=f)
print(gen_classes(particle_db), file=f)
print(inc_end(), file=f) print(inc_end(), file=f)
with open("pythia_db.pkl", "wb") as f: with open("particle_db.pkl", "wb") as f:
serialize_pythia_db(pythia_db, f) serialize_particle_db(particle_db, f)
Framework/Particles/testParticles.cc
+ 23
− 6
View file @ c2cef4eb
...@@ -58,11 +58,28 @@ TEST_CASE("ParticleProperties", "[Particles]") { ...@@ -58,11 +58,28 @@ TEST_CASE("ParticleProperties", "[Particles]") {
REQUIRE(GetLifetime(Code::Electron) == REQUIRE(GetLifetime(Code::Electron) ==
std::numeric_limits<double>::infinity() * corsika::units::si::second); std::numeric_limits<double>::infinity() * corsika::units::si::second);
REQUIRE(GetLifetime(Code::DPlus) < GetLifetime(Code::Gamma)); REQUIRE(GetLifetime(Code::DPlus) < GetLifetime(Code::Gamma));
// REQUIRE(GetLifetime(Code::RhoPlus)/corsika::units::si::second == REQUIRE(GetLifetime(Code::RhoPlus)/corsika::units::si::second ==
// (Approx(4.414566727909413e-24).epsilon(1e-3))); (Approx(4.414566727909413e-24).epsilon(1e-3)));
// REQUIRE(GetLifetime(Code::SigmaMinusBar)/corsika::units::si::second == REQUIRE(GetLifetime(Code::SigmaMinusBar)/corsika::units::si::second ==
// (Approx(8.018880848563575e-11).epsilon(1e-5))); (Approx(8.018880848563575e-11).epsilon(1e-5)));
// REQUIRE(GetLifetime(Code::MuPlus)/corsika::units::si::second == REQUIRE(GetLifetime(Code::MuPlus)/corsika::units::si::second ==
// (Approx(2.1970332555864364e-06).epsilon(1e-5))); (Approx(2.1970332555864364e-06).epsilon(1e-5)));
} }
SECTION("Nuclei") {
REQUIRE(IsNucleus(Code::Gamma) == false);
REQUIRE(IsNucleus(Code::Argon) == true);
REQUIRE(IsNucleus(Code::Proton) == false);
REQUIRE(IsNucleus(Code::Hydrogen) == true);
REQUIRE(Argon::IsNucleus() == true);
REQUIRE(EtaC::IsNucleus() == false);
REQUIRE(GetNucleusA(Code::Hydrogen) == 1);
REQUIRE(GetNucleusA(Code::Tritium) == 3);
REQUIRE(Hydrogen::GetNucleusZ() == 1);
REQUIRE(Tritium::GetNucleusA() == 3);
}
} }
Processes/Sibyll/CMakeLists.txt
+ 2
− 2
View file @ c2cef4eb
add_custom_command ( add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc OUTPUT ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
COMMAND ${PROJECT_SOURCE_DIR}/Processes/Sibyll/code_generator.py COMMAND ${PROJECT_SOURCE_DIR}/Processes/Sibyll/code_generator.py
${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
${PROJECT_SOURCE_DIR}/Processes/Sibyll/sibyll_codes.dat ${PROJECT_SOURCE_DIR}/Processes/Sibyll/sibyll_codes.dat
DEPENDS code_generator.py DEPENDS code_generator.py
sibyll_codes.dat sibyll_codes.dat
${PROJECT_BINARY_DIR}/Framework/Particles/pythia_db.pkl ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
WORKING_DIRECTORY WORKING_DIRECTORY
${PROJECT_BINARY_DIR}/Processes/Sibyll/ ${PROJECT_BINARY_DIR}/Processes/Sibyll/
COMMENT "Generate conversion tables for particle codes SIBYLL <-> CORSIKA" COMMENT "Generate conversion tables for particle codes SIBYLL <-> CORSIKA"
......
Processes/Sibyll/code_generator.py
+ 35
− 35
View file @ c2cef4eb
...@@ -13,16 +13,16 @@ import pickle, sys, itertools ...@@ -13,16 +13,16 @@ import pickle, sys, itertools
# loads the pickled pythia_db (which is an OrderedDict) # loads the pickled particle_db (which is an OrderedDict)
def load_pythiadb(filename): def load_particledb(filename):
with open(filename, "rb") as f: with open(filename, "rb") as f:
pythia_db = pickle.load(f) particle_db = pickle.load(f)
return pythia_db return particle_db
# #
def read_sibyll_codes(filename, pythia_db): def read_sibyll_codes(filename, particle_db):
with open(filename) as f: with open(filename) as f:
for line in f: for line in f:
line = line.strip() line = line.strip()
...@@ -30,19 +30,19 @@ def read_sibyll_codes(filename, pythia_db): ...@@ -30,19 +30,19 @@ def read_sibyll_codes(filename, pythia_db):
continue continue
identifier, sib_code, canInteractFlag, xsType = line.split() identifier, sib_code, canInteractFlag, xsType = line.split()
try: try:
pythia_db[identifier]["sibyll_code"] = int(sib_code) particle_db[identifier]["sibyll_code"] = int(sib_code)
pythia_db[identifier]["sibyll_canInteract"] = int(canInteractFlag) particle_db[identifier]["sibyll_canInteract"] = int(canInteractFlag)
pythia_db[identifier]["sibyll_xsType"] = int(xsType) particle_db[identifier]["sibyll_xsType"] = int(xsType)
except KeyError as e: except KeyError as e:
raise Exception("Identifier '{:s}' not found in pythia_db".format(identifier)) raise Exception("Identifier '{:s}' not found in particle_db".format(identifier))
# generates the enum to access sibyll particles by readable names # generates the enum to access sibyll particles by readable names
def generate_sibyll_enum(pythia_db): def generate_sibyll_enum(particle_db):
output = "enum class SibyllCode : int8_t {\n" output = "enum class SibyllCode : int8_t {\n"
for identifier, pData in pythia_db.items(): for identifier, pData in particle_db.items():
if pData.get('sibyll_code') != None: if pData.get('sibyll_code') != None:
output += " {:s} = {:d},\n".format(identifier, pData['sibyll_code']) output += " {:s} = {:d},\n".format(identifier, pData['sibyll_code'])
output += "};\n" output += "};\n"
...@@ -51,9 +51,9 @@ def generate_sibyll_enum(pythia_db): ...@@ -51,9 +51,9 @@ def generate_sibyll_enum(pythia_db):
# generates the look-up table to convert corsika codes to sibyll codes # generates the look-up table to convert corsika codes to sibyll codes
def generate_corsika2sibyll(pythia_db): def generate_corsika2sibyll(particle_db):
string = "std::array<SibyllCodeIntType, {:d}> constexpr corsika2sibyll = {{\n".format(len(pythia_db)) string = "std::array<SibyllCodeIntType, {:d}> constexpr corsika2sibyll = {{\n".format(len(particle_db))
for identifier, pData in pythia_db.items(): for identifier, pData in particle_db.items():
sibCode = pData.get("sibyll_code", 0) sibCode = pData.get("sibyll_code", 0)
string += " {:d}, // {:s}\n".format(sibCode, identifier if sibCode else identifier + " (not implemented in SIBYLL)") string += " {:d}, // {:s}\n".format(sibCode, identifier if sibCode else identifier + " (not implemented in SIBYLL)")
string += "};\n" string += "};\n"
...@@ -62,9 +62,9 @@ def generate_corsika2sibyll(pythia_db): ...@@ -62,9 +62,9 @@ def generate_corsika2sibyll(pythia_db):
# generates the look-up table to convert corsika codes to sibyll codes # generates the look-up table to convert corsika codes to sibyll codes
def generate_corsika2sibyll_xsType(pythia_db): def generate_corsika2sibyll_xsType(particle_db):
string = "std::array<int, {:d}> constexpr corsika2sibyllXStype = {{\n".format(len(pythia_db)) string = "std::array<int, {:d}> constexpr corsika2sibyllXStype = {{\n".format(len(particle_db))
for identifier, pData in pythia_db.items(): for identifier, pData in particle_db.items():
sibCodeXS = pData.get("sibyll_xsType", -1) sibCodeXS = pData.get("sibyll_xsType", -1)
string += " {:d}, // {:s}\n".format(sibCodeXS, identifier if sibCodeXS else identifier + " (not implemented in SIBYLL)") string += " {:d}, // {:s}\n".format(sibCodeXS, identifier if sibCodeXS else identifier + " (not implemented in SIBYLL)")
string += "};\n" string += "};\n"
...@@ -72,18 +72,18 @@ def generate_corsika2sibyll_xsType(pythia_db): ...@@ -72,18 +72,18 @@ def generate_corsika2sibyll_xsType(pythia_db):
# generates the look-up table to convert sibyll codes to corsika codes # generates the look-up table to convert sibyll codes to corsika codes
def generate_sibyll2corsika(pythia_db) : def generate_sibyll2corsika(particle_db) :
string = "" string = ""
minID = 0 minID = 0
for identifier, pData in pythia_db.items() : for identifier, pData in particle_db.items() :
if 'sibyll_code' in pData: if 'sibyll_code' in pData:
minID = min(minID, pData['sibyll_code']) minID = min(minID, pData['sibyll_code'])
string += "SibyllCodeIntType constexpr minSibyll = {:d};\n\n".format(minID) string += "SibyllCodeIntType constexpr minSibyll = {:d};\n\n".format(minID)
pDict = {} pDict = {}
for identifier, pData in pythia_db.items() : for identifier, pData in particle_db.items() :
if 'sibyll_code' in pData: if 'sibyll_code' in pData:
sib_code = pData['sibyll_code'] - minID sib_code = pData['sibyll_code'] - minID
pDict[sib_code] = identifier pDict[sib_code] = identifier
...@@ -104,13 +104,13 @@ def generate_sibyll2corsika(pythia_db) : ...@@ -104,13 +104,13 @@ def generate_sibyll2corsika(pythia_db) :
# generates the bitset for the flag whether Sibyll knows the particle # generates the bitset for the flag whether Sibyll knows the particle
def generate_known_particle(pythia_db): def generate_known_particle(particle_db):
num_particles = len(pythia_db) num_particles = len(particle_db)
num_bytes = num_particles // 32 + 1 num_bytes = num_particles // 32 + 1
string = "Bitset2::bitset2<{:d}> constexpr isKnown{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes) string = "Bitset2::bitset2<{:d}> constexpr isKnown{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes)
numeric = 0 numeric = 0
for identifier, pData in reversed(pythia_db.items()): for identifier, pData in reversed(particle_db.items()):
handledBySibyll = int("sibyll_code" in pData) & 0x1 handledBySibyll = int("sibyll_code" in pData) & 0x1
numeric = (numeric << 1) | handledBySibyll numeric = (numeric << 1) | handledBySibyll
...@@ -125,14 +125,14 @@ def generate_known_particle(pythia_db): ...@@ -125,14 +125,14 @@ def generate_known_particle(pythia_db):
# generates the bitset for the flag whether Sibyll can use particle as projectile # generates the bitset for the flag whether Sibyll can use particle as projectile
def generate_interacting_particle(pythia_db): def generate_interacting_particle(particle_db):
num_particles = len(pythia_db) num_particles = len(particle_db)
num_bytes = num_particles // 32 + 1 num_bytes = num_particles // 32 + 1
string = "Bitset2::bitset2<{:d}> constexpr canInteract{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes) string = "Bitset2::bitset2<{:d}> constexpr canInteract{{ std::array<uint32_t, {:d}>{{{{\n".format(num_particles, num_bytes)
#string = "std::array<bool, {:d}> constexpr corsika2sibyll = {{\n".format(num_particles) #string = "std::array<bool, {:d}> constexpr corsika2sibyll = {{\n".format(num_particles)
numeric = 0 numeric = 0
for identifier, pData in reversed(pythia_db.items()): for identifier, pData in reversed(particle_db.items()):
can = 0 can = 0
if 'sibyll_canInteract' in pData: if 'sibyll_canInteract' in pData:
if pData['sibyll_canInteract'] > 0: if pData['sibyll_canInteract'] > 0:
...@@ -151,19 +151,19 @@ def generate_interacting_particle(pythia_db): ...@@ -151,19 +151,19 @@ def generate_interacting_particle(pythia_db):
if __name__ == "__main__": if __name__ == "__main__":
if len(sys.argv) != 3: if len(sys.argv) != 3:
print("usage: {:s} <pythia_db.pkl> <sibyll_codes.dat>".format(sys.argv[0]), file=sys.stderr) print("usage: {:s} <particle_db.pkl> <sibyll_codes.dat>".format(sys.argv[0]), file=sys.stderr)
sys.exit(1) sys.exit(1)
print("code_generator.py for SIBYLL") print("code_generator.py for SIBYLL")
pythia_db = load_pythiadb(sys.argv[1]) particle_db = load_particledb(sys.argv[1])
read_sibyll_codes(sys.argv[2], pythia_db) read_sibyll_codes(sys.argv[2], particle_db)
with open("Generated.inc", "w") as f: with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=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_sibyll_enum(particle_db), file=f)
print(generate_corsika2sibyll(pythia_db), file=f) print(generate_corsika2sibyll(particle_db), file=f)
print(generate_known_particle(pythia_db), file=f) print(generate_known_particle(particle_db), file=f)
print(generate_sibyll2corsika(pythia_db), file=f) print(generate_sibyll2corsika(particle_db), file=f)
print(generate_interacting_particle(pythia_db), file=f) print(generate_interacting_particle(particle_db), file=f)
print(generate_corsika2sibyll_xsType(pythia_db), file=f) print(generate_corsika2sibyll_xsType(particle_db), file=f)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum