IAP GITLAB

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

Merge branch '125-add-pdg-to-corsika-id-conversion' into 'master' 

Resolve "add PDG to CORSIKA id conversion"

Closes #125

See merge request AirShowerPhysics/corsika!74
parents 1ce78161 582c6a05
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/Particles/NuclearData.xml
+ 6
− 6
View file @ 7a2ba843
......@@ -18,19 +18,19 @@
<particle id="1000020030" name="helium3" A="3" Z="2" >
</particle>
<particle id="1000020030" name="lithium" A="6" Z="3" >
<particle id="1000030070" name="lithium7" A="7" Z="3" >
</particle>
<particle id="1000020030" name="beryllium" A="9" Z="4" >
<particle id="1000040090" name="beryllium9" A="9" Z="4" >
</particle>
<particle id="1000020030" name="boron" A="10" Z="5" >
<particle id="1000110050" name="boron11" A="11" Z="5" >
</particle>
<particle id="1000060120" name="carbon" A="12" Z="6" >
</particle>
<particle id="1000060120" name="carbon13" A="13" Z="6" >
<particle id="1000060130" name="carbon13" A="13" Z="6" >
</particle>
<particle id="1000070140" name="nitrogen" A="14" Z="7" >
......@@ -39,10 +39,10 @@
<particle id="1000080160" name="oxygen" A="16" Z="8" >
</particle>
<particle id="1000080160" name="fluor" A="18" Z="9" >
<particle id="1000080180" name="fluor" A="18" Z="9" >
</particle>
<particle id="1000100220" name="neon21" A="21" Z="10" >
<particle id="1000100210" name="neon21" A="21" Z="10" >
</particle>
<particle id="1000100220" name="neon" A="22" Z="10" >
......
Framework/Particles/ParticleProperties.cc
+ 11
− 0
View file @ 7a2ba843
......@@ -27,5 +27,16 @@ namespace corsika::particles {
std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p) {
return stream << corsika::particles::GetName(p);
}
Code ConvertFromPDG(PDGCode p) {
static_assert(detail::conversionArray.size() % 2 == 1);
auto constexpr maxPDG{(detail::conversionArray.size() - 1) >> 1};
auto k = static_cast<PDGCodeType>(p);
if (abs(k) <= maxPDG) {
return detail::conversionArray[k + maxPDG];
} else {
return detail::conversionMap.at(p);
}
}
} // namespace corsika::particles
Framework/Particles/ParticleProperties.h
+ 23
− 15
View file @ 7a2ba843
......@@ -40,14 +40,15 @@ namespace corsika::particles {
* The Code enum is the actual place to define CORSIKA 8 particle codes.
*/
enum class Code : int16_t;
enum class PDGCode : int32_t;
using CodeIntType = std::underlying_type<Code>::type;
using PDGCodeType = int32_t;
using PDGCodeType = std::underlying_type<PDGCode>::type;
// forward declarations to be used in GeneratedParticleProperties
int16_t constexpr GetElectricChargeNumber(Code const);
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const);
corsika::units::si::HEPMassType constexpr GetMass(Code const);
PDGCodeType constexpr GetPDG(Code const);
PDGCode constexpr GetPDG(Code const);
constexpr std::string const& GetName(Code const);
corsika::units::si::TimeType constexpr GetLifetime(Code const);
......@@ -58,46 +59,52 @@ namespace corsika::particles {
#include <corsika/particles/GeneratedParticleProperties.inc>
/*!
* returns mass of particle
* returns mass of particle in natural units
*/
corsika::units::si::HEPMassType constexpr GetMass(Code const p) {
return detail::masses[static_cast<CodeIntType const>(p)];
return detail::masses[static_cast<CodeIntType>(p)];
}
PDGCodeType constexpr GetPDG(Code const p) {
return detail::pdg_codes[static_cast<CodeIntType const>(p)];
/*!
* returns PDG id
*/
PDGCode constexpr GetPDG(Code const p) {
return detail::pdg_codes[static_cast<CodeIntType>(p)];
}
/*!
* returns electric charge of particle / (e/3).
* returns electric charge of particle / (e/3), e.g. return 3 for a proton.
*/
int16_t constexpr GetElectricChargeNumber(Code const p) {
return detail::electric_charges[static_cast<CodeIntType const>(p)];
return detail::electric_charges[static_cast<CodeIntType>(p)];
}
/*!
* returns electric charge of particle, e.g. return 1.602e-19_C for a proton.
*/
corsika::units::si::ElectricChargeType constexpr GetElectricCharge(Code const p) {
return GetElectricChargeNumber(p) * (corsika::units::constants::e / 3.);
return GetElectricChargeNumber(p) * (corsika::units::constants::e * (1. / 3.));
}
constexpr std::string const& GetName(Code const p) {
return detail::names[static_cast<CodeIntType const>(p)];
return detail::names[static_cast<CodeIntType>(p)];
}
corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
return detail::lifetime[static_cast<CodeIntType const>(p)] *
return detail::lifetime[static_cast<CodeIntType>(p)] *
corsika::units::si::second;
}
bool constexpr IsNucleus(Code const p) {
return detail::isNucleus[static_cast<CodeIntType const>(p)];
return detail::isNucleus[static_cast<CodeIntType>(p)];
}
int constexpr GetNucleusA(Code const p) {
return detail::nucleusA[static_cast<CodeIntType const>(p)];
return detail::nucleusA[static_cast<CodeIntType>(p)];
}
int constexpr GetNucleusZ(Code const p) {
return detail::nucleusZ[static_cast<CodeIntType const>(p)];
return detail::nucleusZ[static_cast<CodeIntType>(p)];
}
/**
......@@ -105,7 +112,8 @@ namespace corsika::particles {
**/
std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p);
Code ConvertFromPDG(PDGCode);
} // namespace corsika::particles
#endif
Framework/Particles/pdxml_reader.py
+ 65
− 42
View file @ 7a2ba843
......@@ -4,6 +4,7 @@ import sys, math, itertools, re, csv, pprint
import xml.etree.ElementTree as ET
from collections import OrderedDict
import pickle
import io
GeVfm = 0.19732696312541853
c_speed_of_light = 29.9792458e10 # mm / s
......@@ -39,7 +40,7 @@ def parsePythia(filename):
ctau = 0.
else:
print ("missing lifetime: " + str(pdg_id) + " " + str(name))
sys.exit(0)
sys.exit(1)
yield (pdg_id, name, mass, electric_charge, antiName, ctau/c_speed_of_light)
......@@ -93,42 +94,7 @@ def class_names(filename):
pdg_id = int(particle.attrib["pdgID"])
map[pdg_id] = name
return map
##############################################################
#
# Automatically produce a string qualifying as C++ class name
#
# This function produces names of type "DELTA_PLUS_PLUS"
#
def c_identifier(name):
orig = name
name = name.upper()
for c in "() /":
name = name.replace(c, "_")
name = name.replace("BAR", "_BAR")
name = name.replace("0", "_0")
name = name.replace("*", "_STAR")
name = name.replace("'", "_PRIME")
name = name.replace("+", "_PLUS")
name = name.replace("-", "_MINUS")
while True:
tmp = name.replace("__", "_")
if tmp == name:
break
else:
name = tmp
pattern = re.compile(r'^[A-Z_][A-Z_0-9]*$')
if pattern.match(name):
return name.strip("_")
else:
raise Exception("could not generate C identifier for '{:s}'".format(orig))
return map
##############################################################
#
......@@ -163,7 +129,7 @@ def c_identifier_camel(name):
name = tmp
name.strip("_")
# remove all "_", if this does not by accident concatenate two number
# remove all "_", if this does not by accident concatenate two numbers
istart = 0
while True:
i = name.find('_', istart)
......@@ -253,6 +219,47 @@ def read_nuclei_db(filename, particle_db, classnames):
return particle_db
###############################################################
#
# build conversion table PDG -> ngc
#
def gen_conversion_PDG_ngc(particle_db):
# todo: find a optimum value, think about cache miss with array vs lookup time with map
P_MAX = 500 # the maximum PDG code that is still filled into the table
conversionDict = dict()
conversionTable = [None] * (2*P_MAX + 1)
for cId, p in particle_db.items():
pdg = p['pdg']
if abs(pdg) < P_MAX:
if conversionTable[pdg + P_MAX]:
raise Exception("table entry already occupied")
else:
conversionTable[pdg + P_MAX] = cId
else:
if pdg in conversionDict.keys():
raise Exception(f"map entry {pdg} already occupied")
else:
conversionDict[pdg] = cId
output = io.StringIO()
def oprint(*args, **kwargs):
print(*args, **kwargs, file=output)
oprint(f"static std::array<Code, {len(conversionTable)}> constexpr conversionArray {{")
for ngc in conversionTable:
oprint(" Code::{0},".format(ngc if ngc else "Unknown"))
oprint("};")
oprint()
oprint("static std::map<PDGCode, Code> const conversionMap {")
for ngc in conversionDict.values():
oprint(f" {{PDGCode::{ngc}, Code::{ngc}}},")
oprint("};")
oprint()
return output.getvalue()
###############################################################
......@@ -260,7 +267,7 @@ def read_nuclei_db(filename, particle_db, classnames):
# return string with enum of all internal particle codes
#
def gen_internal_enum(particle_db):
string = ("enum class Code : int16_t {\n"
string = ("enum class Code : CodeIntType {\n"
" FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops...
......@@ -277,6 +284,20 @@ def gen_internal_enum(particle_db):
return string
###############################################################
#
# return string with enum of all PDG particle codes
#
def gen_pdg_enum(particle_db):
string = "enum class PDGCode : PDGCodeType {\n"
for cId in particle_db:
pdgCode = particle_db[cId]['pdg']
string += " {key:s} = {code:d},\n".format(key = cId, code = pdgCode)
string += " };\n"
return string
###############################################################
......@@ -296,9 +317,9 @@ def gen_properties(particle_db):
string += "};\n\n"
# PDG code table
string += "static constexpr std::array<PDGCodeType const, size> pdg_codes = {\n"
for p in particle_db.values():
string += " {pdg:d}, // {name:s}\n".format(pdg = p['pdg'], name = p['name'])
string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n"
for p in particle_db.keys():
string += f" PDGCode::{p},\n"
string += "};\n"
# name string table
......@@ -473,8 +494,10 @@ if __name__ == "__main__":
with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f)
print(gen_internal_enum(particle_db), file=f)
print(gen_pdg_enum(particle_db), file=f)
print(detail_start(), file=f)
print(gen_properties(particle_db), file=f)
print(gen_conversion_PDG_ngc(particle_db), file=f)
print(detail_end(), file=f)
print(gen_classes(particle_db), file=f)
print(inc_end(), file=f)
......
Framework/Particles/testParticles.cc
+ 23
− 11
View file @ 7a2ba843
......@@ -50,11 +50,23 @@ TEST_CASE("ParticleProperties", "[Particles]") {
}
SECTION("PDG") {
REQUIRE(GetPDG(Code::PiPlus) == 211);
REQUIRE(GetPDG(Code::DPlus) == 411);
REQUIRE(GetPDG(Code::NuMu) == 14);
REQUIRE(GetPDG(Code::NuE) == 12);
REQUIRE(GetPDG(Code::MuMinus) == 13);
REQUIRE(GetPDG(Code::PiPlus) == PDGCode::PiPlus);
REQUIRE(GetPDG(Code::DPlus) == PDGCode::DPlus);
REQUIRE(GetPDG(Code::NuMu) == PDGCode::NuMu);
REQUIRE(GetPDG(Code::NuE) == PDGCode::NuE);
REQUIRE(GetPDG(Code::MuMinus) == PDGCode::MuMinus);
REQUIRE(static_cast<int>(GetPDG(Code::PiPlus)) == 211);
REQUIRE(static_cast<int>(GetPDG(Code::DPlus)) == 411);
REQUIRE(static_cast<int>(GetPDG(Code::NuMu)) == 14);
REQUIRE(static_cast<int>(GetPDG(Code::NuEBar)) == -12);
REQUIRE(static_cast<int>(GetPDG(Code::MuMinus)) == 13);
}
SECTION("Conversion PDG -> internal") {
REQUIRE(ConvertFromPDG(PDGCode::KStarMinus) == Code::KStarMinus);
REQUIRE(ConvertFromPDG(PDGCode::MuPlus) == Code::MuPlus);
REQUIRE(ConvertFromPDG(PDGCode::SigmaStarCMinusBar) == Code::SigmaStarCMinusBar);
}
SECTION("Lifetimes") {
......@@ -70,12 +82,12 @@ TEST_CASE("ParticleProperties", "[Particles]") {
}
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_FALSE(IsNucleus(Code::Gamma));
REQUIRE(IsNucleus(Code::Argon));
REQUIRE_FALSE(IsNucleus(Code::Proton));
REQUIRE(IsNucleus(Code::Hydrogen));
REQUIRE(Argon::IsNucleus());
REQUIRE_FALSE(EtaC::IsNucleus());
REQUIRE(GetNucleusA(Code::Hydrogen) == 1);
REQUIRE(GetNucleusA(Code::Tritium) == 3);
......
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