diff --git a/Framework/Particles/NuclearData.xml b/Framework/Particles/NuclearData.xml
index 0fff50c33839093fc8bfacb9a7632f36110a2492..73eda0ef2350aa6849df60c61a34231dd9fff869 100644
--- a/Framework/Particles/NuclearData.xml
+++ b/Framework/Particles/NuclearData.xml
@@ -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" >
diff --git a/Framework/Particles/ParticleProperties.cc b/Framework/Particles/ParticleProperties.cc
index 4932492ddf33028b456670e9ca2da0a5b8d0efe2..b90a744bd38f5dc2755a2f10753af2e7c7386081 100644
--- a/Framework/Particles/ParticleProperties.cc
+++ b/Framework/Particles/ParticleProperties.cc
@@ -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
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index ad1a14e620899089e94093722c630e218b8db4c2..2f5365eb0525183b2810da96e3c7b85fe9766570 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -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
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index 99544cf5029c933ef9a43478dae73d28c2efbb56..946940f7dcf54551e2888e438c6eeeae0e1659e1 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/Framework/Particles/pdxml_reader.py
@@ -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)
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index 1ca8b2188b5004c22c3168b5ad90290aa10a8ffb..0603ea03575cb7668bed8e3e276b5c81f60e27c3 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -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);