IAP GITLAB

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

reworked particle properties

parent e23d8a3f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/Particles/ParticleClassNames.xml
+ 10
− 10
View file @ b4b2ffb5
...@@ -3,20 +3,20 @@ ...@@ -3,20 +3,20 @@
<!-- For selected particles it is possible to specify the C++ class <!-- For selected particles it is possible to specify the C++ class
names for a specific unique PDG code --> names for a specific unique PDG code -->
<particle pdgID="0" classname="unknown" /> <!-- VOID in Pythia8 --> <particle pdgID="0" classname="unknown"/> <!-- VOID in Pythia8 -->
<particle pdgID="11" classname="Electron" /> <particle pdgID="11" classname="Electron"/>
<particle pdgID="-11" classname="Positron" /> <particle pdgID="-11" classname="Positron"/>
<particle pdgID="22" classname="Gamma" /> <particle pdgID="22" classname="Gamma"/>
<particle id="130" classname="K0Long"> </particle> <particle pdgID="130" classname="K0Long"/>
<particle id="310" classname="K0Short"> </particle> <particle pdgID="310" classname="K0Short"/>
<particle id="2112" classname="Neutron"> </particle> <particle pdgID="2112" classname="Neutron"/>
<particle id="-2112" classname="AntiNeutron"> </particle> <particle pdgID="-2112" classname="AntiNeutron"/>
<particle id="2212" classname="Proton"> </particle> <particle pdgID="2212" classname="Proton"/>
<particle id="-2212" classname="AntiProton"> </particle> <particle pdgID="-2212" classname="AntiProton"/>
</chapter> </chapter>
Framework/Particles/ParticleProperties.h
+ 36
− 16
View file @ b4b2ffb5
...@@ -4,15 +4,16 @@ ...@@ -4,15 +4,16 @@
Interface to particle properties Interface to particle properties
*/ */
#ifndef _include_Particle_h_ #ifndef _include_ParticleProperties_h_
#define _include_Particle_h_ #define _include_ParticleProperties_h_
#include <array> #include <array>
#include <cstdint> #include <cstdint>
#include <iostream> #include <iostream>
#include <type_traits>
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h> #include <corsika/units/PhysicalUnits.h>
#include <corsika/particles/GeneratedParticleProperties.inc>
/** /**
* @namespace particle * @namespace particle
...@@ -23,31 +24,50 @@ ...@@ -23,31 +24,50 @@
*/ */
namespace corsika::particles { namespace corsika::particles {
enum class Code : int16_t;
using PDGCodeType = int16_t;
using CodeIntType = std::underlying_type<Code>::type;
// forward declarations to be used in GeneratedParticleProperties
int16_t constexpr GetElectricChargeNumber(Code const);
corsika::units::ElectricChargeType constexpr GetElectricCharge(Code const);
corsika::units::MassType constexpr GetMass(Code const);
PDGCodeType constexpr GetPDG(Code const);
std::string const GetName(Code const);
#include <corsika/particles/GeneratedParticleProperties.inc>
/** /*!
* @function GetMass * returns mass of particle
*
* return mass of particle
*/ */
auto constexpr GetMass(Code const p) { return masses[static_cast<uint8_t const>(p)]; } corsika::units::MassType constexpr GetMass(Code const p) {
return masses[static_cast<CodeIntType const>(p)];
}
auto constexpr GetPDG(Code const p) { return pdg_codes[static_cast<uint8_t const>(p)]; } PDGCodeType constexpr GetPDG(Code const p) {
return pdg_codes[static_cast<CodeIntType const>(p)];
}
auto constexpr GetElectricChargeNumber(Code const p) { /*!
return electric_charge[static_cast<uint8_t const>(p)] / 3; * returns electric charge of particle / (e/3).
*/
int16_t constexpr GetElectricChargeNumber(Code const p) {
return electric_charges[static_cast<CodeIntType const>(p)];
} }
auto constexpr GetElectricCharge(Code const p) { corsika::units::ElectricChargeType constexpr GetElectricCharge(Code const p) {
return GetElectricChargeNumber(p) * (corsika::units::constants::e); return GetElectricChargeNumber(p) * (corsika::units::constants::e / 3.);
} }
auto const GetName(Code const p) { return names[static_cast<uint8_t const>(p)]; } std::string const GetName(Code const p) {
return names[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) {
stream << GetName(p); return stream << GetName(p);
return stream;
} }
} // namespace io } // namespace io
......
Framework/Particles/pdxml_reader.py
+ 45
− 60
View file @ b4b2ffb5
...@@ -51,7 +51,7 @@ def class_names(filename): ...@@ -51,7 +51,7 @@ def class_names(filename):
for particle in root.iter("particle"): for particle in root.iter("particle"):
name = particle.attrib["classname"] name = particle.attrib["classname"]
pdg_id = int(particle.attrib["id"]) pdg_id = int(particle.attrib["pdgID"])
map[pdg_id] = name map[pdg_id] = name
return map return map
...@@ -97,9 +97,10 @@ def c_identifier(name): ...@@ -97,9 +97,10 @@ def c_identifier(name):
# #
# This function produces names of type "DeltaPlusPlus" # This function produces names of type "DeltaPlusPlus"
# #
def c_identifier_cammel(name): def c_identifier_camel(name):
orig = name orig = name
name = name[0].upper() + name[1:].lower() # all lower case name = name[0].upper() + name[1:].lower() # all lower case
for c in "() /": # replace funny characters for c in "() /": # replace funny characters
name = name.replace(c, "_") name = name.replace(c, "_")
...@@ -111,9 +112,9 @@ def c_identifier_cammel(name): ...@@ -111,9 +112,9 @@ def c_identifier_cammel(name):
# move "Bar" to end of name # move "Bar" to end of name
ibar = name.find('Bar') ibar = name.find('Bar')
if (ibar>0 and ibar<len(name)-3) : if ibar > 0 and ibar < len(name)-3:
name = name[:ibar] + name[ibar+3:] + str('Bar') name = name[:ibar] + name[ibar+3:] + str('Bar')
# cleanup "_"s # cleanup "_"s
while True: while True:
tmp = name.replace("__", "_") tmp = name.replace("__", "_")
...@@ -127,18 +128,18 @@ def c_identifier_cammel(name): ...@@ -127,18 +128,18 @@ def c_identifier_cammel(name):
istart = 0 istart = 0
while True: while True:
i = name.find('_', istart) i = name.find('_', istart)
if (i<1 or i>len(name)-1): if i < 1 or i > len(name)-1:
break break
istart = i istart = i
if (name[i-1].isdigit() and name[i+1].isdigit()): if name[i-1].isdigit() and name[i+1].isdigit():
# there is a number on both sides # there is a number on both sides
break break
name = name[:i] + name[i+1:] name = name[:i] + name[i+1:]
# and last, for example: make NuE out of Nue # and last, for example: make NuE out of Nue
if (name[i-1].islower() and name[i].islower()) : if name[i-1].islower() and name[i].islower():
if (i<len(name)-1) : if i < len(name)-1:
name = name[:i] + name[i].upper() + name[i+1:] name = name[:i] + name[i].upper() + name[i+1:]
else : else:
name = name[:i] + name[i].upper() name = name[:i] + name[i].upper()
# check if name is valid C++ identifier # check if name is valid C++ identifier
...@@ -146,27 +147,25 @@ def c_identifier_cammel(name): ...@@ -146,27 +147,25 @@ def c_identifier_cammel(name):
if pattern.match(name): if pattern.match(name):
return name return name
else: else:
raise Exception("could not generate C identifier for '{:s}' name='{:s}'".format(orig, name)) raise Exception("could not generate C identifier for '{:s}': result '{:s}'".format(orig, name))
# ######################################################## ##########################################################
# #
# returns dict containing all data from pythia-xml input # returns dict containing all data from pythia-xml input
# #
def build_pythia_db(filename, classnames): def build_pythia_db(filename, classnames):
particle_db = OrderedDict() particle_db = OrderedDict()
for (pdg, name, mass, electric_charge, antiName) in parse(filename): for (pdg, name, mass, electric_charge, antiName) in parse(filename):
c_id = "unknown" c_id = "unknown"
if (pdg in classnames): if pdg in classnames:
c_id = classnames[pdg] c_id = classnames[pdg]
else: else:
c_id = c_identifier_cammel(name) # the cammel case names c_id = c_identifier_camel(name) # the camel case names
#~ print(name, c_id, sep='\t', file=sys.stderr)
#~ enums += "{:s} = {:d}, ".format(c_id, corsika_id)
particle_db[c_id] = { particle_db[c_id] = {
"name" : name, "name" : name,
"antiName" : antiName, "antiName" : antiName,
...@@ -186,17 +185,20 @@ def build_pythia_db(filename, classnames): ...@@ -186,17 +185,20 @@ def build_pythia_db(filename, classnames):
# #
def gen_internal_enum(pythia_db): def gen_internal_enum(pythia_db):
string = "enum class Code : uint8_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...
string += " FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n" # identifier for eventual loops...
last_ngc_id = 0
for k in filter(lambda k: "ngc_code" in pythia_db[k], pythia_db): for k in filter(lambda k: "ngc_code" in pythia_db[k], pythia_db):
last_ngc_id = pythia_db[k]['ngc_code'] last_ngc_id = pythia_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 = " + str(last_ngc_id+1) + ",\n" # identifier for eventual loops... string += (" LastParticle = {:d},\n" # identifier for eventual loops...
"}};").format(last_ngc_id + 1)
if last_ngc_id > 0x7fff: # does not fit into int16_t
raise Exception("Integer overflow in internal particle code definition prevented!")
string += "};"
return string return string
...@@ -214,13 +216,13 @@ def gen_properties(pythia_db): ...@@ -214,13 +216,13 @@ def gen_properties(pythia_db):
string += "\n" string += "\n"
# particle masses table # particle masses table
string += "static constexpr std::array<quantity<energy_d> const, size> masses = {\n" string += "static constexpr std::array<corsika::units::MassType const, size> masses = {\n"
for p in pythia_db.values(): for p in pythia_db.values():
string += " {mass:f}_GeV, // {name:s}\n".format(mass = p['mass'], name = p['name']) string += " {mass:f} * (1e9 * corsika::units::constants::eV / corsika::units::constants::cSquared), // {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<PDGCode 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 pythia_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"
...@@ -232,7 +234,7 @@ def gen_properties(pythia_db): ...@@ -232,7 +234,7 @@ def gen_properties(pythia_db):
string += "};\n" string += "};\n"
# electric charges table # electric charges table
string += "static constexpr std::array<int16_t, size> electric_charge = {\n" string += "static constexpr std::array<int16_t, size> electric_charges = {\n"
for p in pythia_db.values(): for p in pythia_db.values():
string += " {charge:d},\n".format(charge = p['electric_charge']) string += " {charge:d},\n".format(charge = p['electric_charge'])
string += "};\n" string += "};\n"
...@@ -269,23 +271,23 @@ def gen_classes(pythia_db): ...@@ -269,23 +271,23 @@ def gen_classes(pythia_db):
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(pythia_db[cname]['pdg']) +"\n"
string += " * - mass=" + str(pythia_db[cname]['mass']) + " GeV \n" string += " * - mass=" + str(pythia_db[cname]['mass']) + " GeV/c² \n"
string += " * - charge= " + str(pythia_db[cname]['electric_charge']/3) + " \n" string += " * - charge= " + str(pythia_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"
string += "*/\n\n" string += "*/\n\n"
string += "class " + cname + "{\n" string += "class " + cname + " {\n"
string += " public:\n" string += " public:\n"
string += " static Code GetCode() { return Type; }\n" string += " static constexpr Code GetCode() { return Type; }\n"
string += " static quantity<energy_d> GetMass() { return masses[TypeIndex]; }\n" string += " static constexpr corsika::units::MassType GetMass() { return corsika::particles::GetMass(Type); }\n"
string += " static quantity<electric_charge_d> GetCharge() { return corsika::units::constants::e*electric_charge[TypeIndex]/3; }\n" string += " static constexpr corsika::units::ElectricChargeType GetCharge() { return corsika::particles::GetElectricCharge(Type); }\n"
string += " static int GetChargeNumber() { return electric_charge[TypeIndex]/3; }\n" string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetElectricChargeNumber(Type); }\n"
string += " static std::string GetName() { return names[TypeIndex]; }\n" string += " static std::string const GetName() { return corsika::particles::GetName(Type); }\n"
string += " static Code GetAntiParticle() { return AntiType; }\n" string += " static constexpr Code GetAntiParticle() { return AntiType; }\n"
string += " static const Code Type = Code::" + cname + ";\n" string += " static constexpr Code Type = Code::" + cname + ";\n"
string += " static const Code AntiType = Code::" + antiP + ";\n" string += " static constexpr Code AntiType = Code::" + antiP + ";\n"
string += " private:\n" string += " private:\n"
string += " static const uint8_t TypeIndex = static_cast<uint8_t const>(Type);\n" string += " static constexpr CodeIntType TypeIndex = static_cast<CodeIntType const>(Type);\n"
string += "};\n" string += "};\n"
return string return string
...@@ -296,22 +298,8 @@ def gen_classes(pythia_db): ...@@ -296,22 +298,8 @@ def gen_classes(pythia_db):
# #
def inc_start(): def inc_start():
string = "" string = ('// generated by pdxml_reader.py\n'
string += "#ifndef _include_GeneratedParticleDataTable_h_\n" '// MANUAL EDITS ON OWN RISK\n')
string += "#define _include_GeneratedParticleDataTable_h_\n\n"
string += "#include <corsika/units/PhysicalUnits.h>\n"
string += "#include <corsika/units/PhysicalConstants.h>\n"
string += "#include <array>\n"
string += "#include <cstdint>\n"
# string += "#include <iostream>\n\n"
string += "namespace corsika { \n\n"
# string += "using namespace literals; \n"
string += "namespace particles { \n\n"
string += "using corsika::units::energy_d;\n"
string += "using corsika::units::electric_charge_d;\n"
string += "using corsika::units::quantity;\n"
string += "using corsika::units::operator\"\"_GeV;\n"
string += "typedef int16_t PDGCode;\n\n"
return string return string
...@@ -321,9 +309,6 @@ def inc_start(): ...@@ -321,9 +309,6 @@ def inc_start():
def inc_end(): def inc_end():
string = "" string = ""
string += "\n}\n\n"
string += "\n}\n\n"
string += "#endif\n"
return string return string
...@@ -334,14 +319,14 @@ def inc_end(): ...@@ -334,14 +319,14 @@ def inc_end():
if __name__ == "__main__": if __name__ == "__main__":
if (len(sys.argv)!=3) : if len(sys.argv) != 3:
print ("pdxml_reader.py Pythia8.xml ClassNames.xml") print("usage: {:s} <Pythia8.xml> <ClassNames.xml>".format(sys.argv[0]))
sys.exit(0) sys.exit(1)
names = class_names(sys.argv[2]) names = class_names(sys.argv[2])
pythia_db = build_pythia_db(sys.argv[1], names) pythia_db = build_pythia_db(sys.argv[1], names)
print ("\n pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n") print("\n pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n")
counter = itertools.count(0) counter = itertools.count(0)
......
Framework/Particles/testParticles.cc
+ 1
− 1
View file @ b4b2ffb5
...@@ -14,7 +14,7 @@ TEST_CASE("Particles", "[Particles]") { ...@@ -14,7 +14,7 @@ TEST_CASE("Particles", "[Particles]") {
SECTION("Types") { REQUIRE(Electron::GetCode() == Code::Electron); } SECTION("Types") { REQUIRE(Electron::GetCode() == Code::Electron); }
SECTION("Data") { SECTION("Data") {
REQUIRE(Electron::GetMass() / 0.511_MeV == Approx(1)); REQUIRE(Electron::GetMass() / (511_keV / constants::cSquared) == Approx(1));
REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1)); REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
REQUIRE(Electron::GetCharge() / constants::e == Approx(-1)); REQUIRE(Electron::GetCharge() / constants::e == Approx(-1));
REQUIRE(Positron::GetCharge() / constants::e == Approx(+1)); REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
......
Framework/Units/PhysicalConstants.h
+ 3
− 2
View file @ b4b2ffb5
...@@ -39,16 +39,17 @@ namespace corsika::units::constants { ...@@ -39,16 +39,17 @@ namespace corsika::units::constants {
// Avogadro constant // Avogadro constant
constexpr quantity<dimensions<0, 0, 0, 0, 0, -1> > N_sub_A{Rep(6.02214199e+23L) / mole}; constexpr quantity<dimensions<0, 0, 0, 0, 0, -1> > N_sub_A{Rep(6.02214199e+23L) / mole};
// electronvolt // electronvolt
constexpr quantity<energy_d> eV{Rep(1.60217733e-19L) * joule}; constexpr quantity<energy_d> eV{Rep(1.6021766208e-19L) * joule};
// elementary charge // elementary charge
constexpr quantity<electric_charge_d> e{Rep(1.602176462e-19L) * coulomb}; constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb};
// Planck constant // Planck constant
constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second}; constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second};
// speed of light in a vacuum // speed of light in a vacuum
constexpr quantity<speed_d> c{Rep(299792458L) * meter / second}; constexpr quantity<speed_d> c{Rep(299792458L) * meter / second};
constexpr auto cSquared = c * c;
// unified atomic mass unit // unified atomic mass unit
constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram}; constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
......
Framework/Units/PhysicalUnits.h
+ 1
− 1
View file @ b4b2ffb5
...@@ -22,7 +22,6 @@ namespace phys { ...@@ -22,7 +22,6 @@ namespace phys {
} // namespace phys } // namespace phys
namespace corsika::units { namespace corsika::units {
using namespace phys::units; using namespace phys::units;
using namespace phys::units::literals; using namespace phys::units::literals;
// namespace literals = phys::units::literals; // namespace literals = phys::units::literals;
...@@ -34,6 +33,7 @@ namespace corsika::units { ...@@ -34,6 +33,7 @@ namespace corsika::units {
using ElectricChargeType = using ElectricChargeType =
phys::units::quantity<phys::units::electric_charge_d, double>; phys::units::quantity<phys::units::electric_charge_d, double>;
using EnergyType = phys::units::quantity<phys::units::energy_d, double>; using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
using MassType = phys::units::quantity<phys::units::mass_d, double>;
} // end namespace corsika::units } // end namespace corsika::units
......
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