From 7f37d1f9a2f0390be56ee465f58ae3c6120934d2 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 4 Dec 2018 22:20:59 +0100
Subject: [PATCH] added particle lifetimes
---
Framework/Particles/ParticleProperties.h | 8 ++++
Framework/Particles/pdxml_reader.py | 48 +++++++++++++++++-------
Framework/Particles/testParticles.cc | 8 ++++
3 files changed, 51 insertions(+), 13 deletions(-)
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index 330048b3..4bb3cda1 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -26,6 +26,7 @@
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
+
/**
* @namespace particle
*
@@ -36,6 +37,8 @@
namespace corsika::particles {
+ using corsika::units::si::second;
+
enum class Code : int16_t;
using PDGCodeType = int16_t;
@@ -47,6 +50,7 @@ namespace corsika::particles {
corsika::units::si::MassType constexpr GetMass(Code const);
PDGCodeType constexpr GetPDG(Code const);
constexpr std::string const& GetName(Code const);
+ corsika::units::si::TimeType constexpr GetLifetime(Code const);
#include <corsika/particles/GeneratedParticleProperties.inc>
@@ -76,6 +80,10 @@ namespace corsika::particles {
return names[static_cast<CodeIntType const>(p)];
}
+ corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
+ return lifetime[static_cast<CodeIntType const>(p)];
+ }
+
namespace io {
std::ostream& operator<<(std::ostream& stream, Code const p);
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index 96d50588..1483a6de 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/Framework/Particles/pdxml_reader.py
@@ -14,27 +14,37 @@ import pickle
def parse(filename):
tree = ET.parse(filename)
root = tree.getroot()
-
+
+ GeVfm = 0.19732696312541853
+ c_speed_of_light = 29.9792458e10 # mm / s
+
for particle in root.iter("particle"):
name = particle.attrib["name"]
antiName = "Unknown"
- # print (str(particle.attrib))
if ("antiName" in particle.attrib):
- antiName = particle.attrib["antiName"]
-# print ("found anti: " + name + " " + antiName + "\n" )
+ antiName = particle.attrib["antiName"]
pdg_id = int(particle.attrib["id"])
- mass = float(particle.attrib["m0"]) # GeV
- electric_charge = int(particle.attrib["chargeType"]) # in units of e/3
-
- decay_width = float(particle.attrib.get("mWidth", 0)) # GeV
- lifetime = float(particle.attrib.get("tau0", math.inf)) # mm / c
+ mass = float(particle.attrib["m0"]) # GeV
+ electric_charge = int(particle.attrib["chargeType"]) # in units of e/3
+ ctau = 0.
+ if pdg_id in (11, 12, 14, 16, 22, 2212): # these are the stable particles !
+ ctau = float('Inf')
+ elif 'tau0' in particle.attrib:
+ ctau = float(particle.attrib['tau0']) # mm / c
+ elif 'mWidth' in particle.attrib:
+ ctau = GeVfm / float(particle.attrib['mWidth']) * 1e-15 * 1000.0 # mm / s
+ elif pdg_id in (0, 423, 433, 4312, 4322, 5112, 5222): # those are certainly not stable....
+ ctau = 0.
+ else:
+ print ("missing lifetime: " + str(pdg_id) + " " + str(name))
+ sys.exit(0)
- yield (pdg_id, name, mass, electric_charge, antiName)
+ yield (pdg_id, name, mass, electric_charge, antiName, ctau/c_speed_of_light)
# TODO: read decay channels from child elements
if "antiName" in particle.attrib:
- yield (-pdg_id, antiName, mass, -electric_charge, name)
+ yield (-pdg_id, antiName, mass, -electric_charge, name, ctau/c_speed_of_light)
@@ -161,7 +171,7 @@ def build_pythia_db(filename, classnames):
counter = itertools.count(0)
- for (pdg, name, mass, electric_charge, antiName) in parse(filename):
+ for (pdg, name, mass, electric_charge, antiName, lifetime) in parse(filename):
c_id = "Unknown"
if pdg in classnames:
@@ -175,6 +185,7 @@ def build_pythia_db(filename, classnames):
"pdg" : pdg,
"mass" : mass, # in GeV
"electric_charge" : electric_charge, # in e/3
+ "lifetime" : lifetime,
"ngc_code" : next(counter)
}
@@ -249,6 +260,15 @@ def gen_properties(pythia_db):
# string += " {anti:d},\n".format(charge = p['anti_particle'])
# string += "};\n"
+ # lifetime
+ string += "static constexpr std::array<corsika::units::si::TimeType const, size> lifetime = {\n"
+ for p in pythia_db.values():
+ if p['lifetime'] == float("Inf") :
+ string += " std::numeric_limits<double>::infinity() * corsika::units::si::second, \n"
+ else :
+ string += " {tau:f} * corsika::units::si::second, \n".format(tau = p['lifetime'])
+ string += "};\n"
+
return string
@@ -335,10 +355,12 @@ if __name__ == "__main__":
print("usage: {:s} <Pythia8.xml> <ClassNames.xml>".format(sys.argv[0]), file=sys.stderr)
sys.exit(1)
+ print("\n pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n")
+
names = class_names(sys.argv[2])
pythia_db = build_pythia_db(sys.argv[1], names)
- print("\n pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n")
+ print (str(pythia_db))
with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f)
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index 443d167c..283bb48c 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -53,4 +53,12 @@ TEST_CASE("ParticleProperties", "[Particles]") {
REQUIRE(GetPDG(Code::NuE) == 12);
REQUIRE(GetPDG(Code::MuMinus) == 13);
}
+
+ SECTION("Lifetimes") {
+ REQUIRE(GetLifetime(Code::Electron) == std::numeric_limits<double>::infinity() * corsika::units::si::second);
+ REQUIRE(GetLifetime(Code::DPlus) < GetLifetime(Code::Gamma));
+ //REQUIRE(GetLifetime(Code::RhoPlus)/corsika::units::si::second == (Approx(4.414566727909413e-24).epsilon(1e-3)));
+ //REQUIRE(GetLifetime(Code::SigmaMinusBar)/corsika::units::si::second == (Approx(8.018880848563575e-11).epsilon(1e-5)));
+ //REQUIRE(GetLifetime(Code::MuPlus)/corsika::units::si::second == (Approx(2.1970332555864364e-06).epsilon(1e-5)));
+ }
}
--
GitLab