diff --git a/CMakeLists.txt b/CMakeLists.txt
index c436d8906e31f56c3d010e255736cc16c87a7050..5a38c4ecb422e5ea823f864deaf0df7acceb5ac9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -14,7 +14,7 @@ set (CMAKE_INSTALL_MESSAGE LAZY)
set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/CMakeModules)
include (CorsikaUtilities) # a few cmake function
-# --std=c++14
+# --std=c++17
set (CMAKE_CXX_STANDARD 17)
enable_testing ()
set (CTEST_OUTPUT_ON_FAILURE 1)
diff --git a/Framework/Geometry/BaseTrajectory.h b/Framework/Geometry/BaseTrajectory.h
new file mode 100644
index 0000000000000000000000000000000000000000..ec5e4062cfe2be63d1f7be83bdb4a3e02cc1d596
--- /dev/null
+++ b/Framework/Geometry/BaseTrajectory.h
@@ -0,0 +1,23 @@
+#ifndef _include_BASETRAJECTORY_H
+#define _include_BASETRAJECTORY_H
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <string>
+
+namespace corsika::geometry {
+
+ /*!
+ * Base class for trajectories.
+ */
+ class BaseTrajectory {
+
+ public:
+ //!< t for \f$ t = 0 \f$, the starting Point shall be returned.
+ virtual Point GetPosition(corsika::units::TimeType t) const = 0;
+ };
+
+} // namespace corsika::geometry
+
+#endif
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index 36fab64cc742c8eb292ccb2420b0e7f3f6bebbd3..a4fc43cb3364881d7f3e42939b0985b69df4afab 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -13,7 +13,9 @@ set (
Helix.h
BaseVector.h
QuantityVector.h
+ BaseTrajectory.h
LineTrajectory.h
+ Trajectory.h
)
set (
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index ca97fa8bf8c4100cb9031607485ce291c76a636e..29bb49438f23387413094c8319486db8eea67e70 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -1,31 +1,38 @@
#ifndef _include_HELIX_H_
#define _include_HELIX_H_
+#include <corsika/geometry/BaseTrajectory.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
namespace corsika::geometry {
+ /*!
+ * A Helix is defined by the cyclotron frequency \f$ \omega_c \f$, the initial
+ * Point r0 and
+ * the velocity vectors \f$ \vec{v}_{\parallel} \f$ and \f$ \vec{v}_{\perp} \f$
+ * denoting the projections of the initial velocity \f$ \vec{v}_0 \f$ parallel
+ * and perpendicular to the axis \f$ \vec{B} \f$, respectively, i.e.
+ * \f{align*}{
+ \vec{v}_{\parallel} &= \frac{\vec{v}_0 \cdot \vec{B}}{\vec{B}^2} \vec{B} \\
+ \vec{v}_{\perp} &= \vec{v}_0 - \vec{v}_{\parallel}
+ \f}
+ */
+
+ class Helix : public BaseTrajectory {
+ using VelocityVec = Vector<corsika::units::SpeedType::dimension_type>;
- using corsika::units::frequency_d;
- using corsika::units::FrequencyType;
- using corsika::units::quantity;
- using corsika::units::SpeedType;
- using corsika::units::TimeType;
-
- class Helix // TODO: inherit from to-be-implemented "Trajectory"
- {
- using SpeedVec = Vector<SpeedType::dimension_type>;
Point const r0;
- FrequencyType const omegaC;
- SpeedVec const vPar;
- SpeedVec vPerp, uPerp;
- LengthType const radius;
+ corsika::units::FrequencyType const omegaC;
+ VelocityVec const vPar;
+ VelocityVec const vPerp, uPerp;
+
+ corsika::units::LengthType const radius;
public:
- Helix(Point const& pR0, quantity<frequency_d> pOmegaC, SpeedVec const& pvPar,
- SpeedVec const& pvPerp)
+ Helix(Point const& pR0, corsika::units::FrequencyType pOmegaC,
+ VelocityVec const& pvPar, VelocityVec const& pvPerp)
: r0(pR0)
, omegaC(pOmegaC)
, vPar(pvPar)
@@ -33,7 +40,7 @@ namespace corsika::geometry {
, uPerp(vPerp.cross(vPar.normalized()))
, radius(pvPar.norm() / abs(pOmegaC)) {}
- auto GetPosition(TimeType t) const {
+ Point GetPosition(corsika::units::TimeType t) const {
return r0 + vPar * t +
(vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
}
diff --git a/Framework/Geometry/LineTrajectory.h b/Framework/Geometry/LineTrajectory.h
index fc550d4bdd87c63af1bb307721ba4b675e245c57..556da4642c52689c492fbd619863a3bfc0729d99 100644
--- a/Framework/Geometry/LineTrajectory.h
+++ b/Framework/Geometry/LineTrajectory.h
@@ -1,25 +1,25 @@
#ifndef _include_LINETRAJECTORY_H
#define _include_LINETRAJECTORY_H
+#include <corsika/geometry/BaseTrajectory.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
- class LineTrajectory // TODO: inherit from Trajectory
- {
- using SpeedVec = Vector<corsika::units::SpeedType::dimension_type>;
+ class LineTrajectory : public BaseTrajectory {
+ using VelocityVec = Vector<corsika::units::SpeedType::dimension_type>;
Point const r0;
- SpeedVec const v0;
+ VelocityVec const v0;
public:
- LineTrajectory(Point const& pR0, SpeedVec const& pV0)
- : r0(r0)
+ LineTrajectory(Point const& pR0, VelocityVec const& pV0)
+ : r0(pR0)
, v0(pV0) {}
- auto GetPosition(corsika::units::TimeType t) const { return r0 + v0 * t; }
+ Point GetPosition(corsika::units::TimeType t) const override { return r0 + v0 * t; }
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index c5af419ee4df97ec731a7fb0c95123fc2c843de5..671a6856092a0811aaf5d407ef68b0f7fbd37e4f 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -101,6 +101,8 @@ namespace corsika {
auto& operator-() const { return QuantityVector<dim>(-eVector); }
auto normalized() const { return (*this) * (1 / norm()); }
+
+ auto operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
};
} // end namespace corsika
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 552dfde6f1341a7c0bbfb3903e27c8fbab5e1e05..28b4e7dc63d5ded1ecd6d047df442af4557e4f20 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -15,6 +15,7 @@ namespace corsika::geometry {
: center(pCenter)
, radius(pRadius) {}
+ //! returns true if the Point p is within the sphere
auto isInside(Point const& p) const {
return radius * radius > (center - p).squaredNorm();
}
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a640cacd027510a64d0c574c48e66dbe11fd1d9
--- /dev/null
+++ b/Framework/Geometry/Trajectory.h
@@ -0,0 +1,31 @@
+#ifndef _include_TRAJECTORY_H
+#define _include_TRAJECTORY_H
+
+#include <corsika/geometry/BaseTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::geometry {
+
+ class Trajectory {
+ corsika::units::TimeType const fTStart, fTEnd;
+ BaseTrajectory const& fTrajectory;
+
+ public:
+ Trajectory(corsika::units::TimeType pTStart, corsika::units::TimeType pTEnd,
+ BaseTrajectory const& pTrajectory)
+ : fTStart(pTStart)
+ , fTEnd(pTEnd)
+ , fTrajectory(pTrajectory) {}
+
+ Point GetPosition(corsika::units::TimeType t) const {
+ return fTrajectory.GetPosition(t + fTStart);
+ }
+
+ Point GetPosition(double u) const {
+ return GetPosition(fTEnd * u + fTStart * (1 - u));
+ }
+ };
+
+} // namespace corsika::geometry
+
+#endif
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index a2f5181dc6f837b8608563541c3b6cca80b29248..469d2dec68643215ac7a0d33bb55edbee3e211ad 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -2,19 +2,166 @@
// cpp file
#include <catch2/catch.hpp>
+#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/Helix.h>
+#include <corsika/geometry/LineTrajectory.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
+#include <cmath>
-using namespace phys::units;
-using namespace phys::units::literals;
+using namespace corsika::geometry;
+using namespace corsika::units;
-TEST_CASE("Geometry", "[Geometry]") {
- SECTION("Coordinate Systems") {}
+double constexpr absMargin = 1.0e-8;
+;
- SECTION("Point") {}
+TEST_CASE("transformations between CoordinateSystems") {
+ CoordinateSystem rootCS;
- SECTION("Vector") {}
+ REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS)
+ .isApprox(EigenTransform::Identity()));
- SECTION("Helix") {}
+ corsika::QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+ Point p1(rootCS, coordinates);
- SECTION("Line") {}
+ corsika::QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla,
+ 0. * tesla};
+ Vector<magnetic_flux_density_d> v1(rootCS, components);
+
+ REQUIRE((p1.GetCoordinates() - coordinates).norm().magnitude() ==
+ Approx(0).margin(absMargin));
+ REQUIRE((p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() ==
+ Approx(0).margin(absMargin));
+
+ SECTION("unconnected CoordinateSystems") {
+ CoordinateSystem rootCS2;
+ REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2));
+ }
+
+ SECTION("translations") {
+ corsika::QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
+
+ CoordinateSystem translatedCS = rootCS.translate(translationVector);
+
+ REQUIRE(translatedCS.GetReference() == &rootCS);
+
+ REQUIRE((p1.GetCoordinates(translatedCS) + translationVector).norm().magnitude() ==
+ Approx(0).margin(absMargin));
+
+ // Vectors are not subject to translations
+ REQUIRE(
+ (v1.GetComponents(rootCS) - v1.GetComponents(translatedCS)).norm().magnitude() ==
+ Approx(0).margin(absMargin));
+
+ Point p2(translatedCS, {0_m, 0_m, 0_m});
+ REQUIRE(((p2 - p1).GetComponents() - translationVector).norm().magnitude() ==
+ Approx(0).margin(absMargin));
+ }
+
+ SECTION("multiple translations") {
+ corsika::QuantityVector<length_d> const tv1{0_m, 5_m, 0_m};
+ CoordinateSystem cs2 = rootCS.translate(tv1);
+
+ corsika::QuantityVector<length_d> const tv2{3_m, 0_m, 0_m};
+ CoordinateSystem cs3 = rootCS.translate(tv2);
+
+ corsika::QuantityVector<length_d> const tv3{0_m, 0_m, 2_m};
+ CoordinateSystem cs4 = cs3.translate(tv3);
+
+ REQUIRE(cs4.GetReference()->GetReference() == &rootCS);
+
+ REQUIRE(CoordinateSystem::GetTransformation(cs3, cs2).isApprox(
+ rootCS.translate({3_m, -5_m, 0_m}).GetTransform()));
+ REQUIRE(CoordinateSystem::GetTransformation(cs2, cs3).isApprox(
+ rootCS.translate({-3_m, +5_m, 0_m}).GetTransform()));
+ }
+
+ SECTION("rotations") {
+ corsika::QuantityVector<length_d> const axis{0_m, 0_m, 1_km};
+ double const angle = 90. / 180. * M_PI;
+
+ CoordinateSystem rotatedCS = rootCS.rotate(axis, angle);
+ REQUIRE(rotatedCS.GetReference() == &rootCS);
+
+ REQUIRE(v1.GetComponents(rotatedCS)[1].magnitude() ==
+ Approx((-1. * tesla).magnitude()));
+
+ // vector norm invariant under rotation
+ REQUIRE(v1.GetComponents(rotatedCS).norm().magnitude() ==
+ Approx(v1.GetComponents(rootCS).norm().magnitude()));
+ }
+
+ SECTION("multiple rotations") {
+ corsika::QuantityVector<length_d> const zAxis{0_m, 0_m, 1_km};
+ corsika::QuantityVector<length_d> const yAxis{0_m, 7_nm, 0_m};
+ corsika::QuantityVector<length_d> const xAxis{2_m, 0_nm, 0_m};
+
+ double const angle = 90. / 180. * M_PI;
+
+ CoordinateSystem rotated1 = rootCS.rotate(zAxis, angle);
+ CoordinateSystem rotated2 = rotated1.rotate(yAxis, angle);
+ CoordinateSystem rotated3 = rotated2.rotate(zAxis, -angle);
+
+ CoordinateSystem combined = rootCS.rotate(xAxis, -angle);
+
+ auto comp1 = v1.GetComponents(rootCS);
+ auto comp3 = v1.GetComponents(combined);
+ REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0).margin(absMargin));
+ }
+}
+
+TEST_CASE("Sphere") {
+ CoordinateSystem rootCS;
+ Point center(rootCS, {0_m, 3_m, 4_m});
+ Sphere sphere(center, 5_m);
+
+ SECTION("isInside") {
+ REQUIRE_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m})));
+ REQUIRE(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m})));
+ }
+}
+
+TEST_CASE("Trajectories") {
+ CoordinateSystem rootCS;
+ Point r0(rootCS, {0_m, 0_m, 0_m});
+
+ SECTION("Line") {
+ Vector<SpeedType::dimension_type> v0(rootCS,
+ {1_m / second, 0_m / second, 0_m / second});
+
+ LineTrajectory const lineTrajectory(r0, v0);
+ CHECK((lineTrajectory.GetPosition(2_s).GetCoordinates() -
+ corsika::QuantityVector<length_d>(2_m, 0_m, 0_m))
+ .norm()
+ .magnitude() == Approx(0).margin(absMargin));
+
+ BaseTrajectory const* base = &lineTrajectory;
+ CHECK(lineTrajectory.GetPosition(2_s).GetCoordinates() ==
+ base->GetPosition(2_s).GetCoordinates());
+ }
+
+ SECTION("Helix") {
+ Vector<SpeedType::dimension_type> const vPar(
+ rootCS, {0_m / second, 0_m / second, 4_m / second}),
+ vPerp(rootCS, {1_m / second, 0_m / second, 0_m / second});
+ auto const omegaC = 2 * M_PI / 1_s;
+
+ Helix const helix(r0, omegaC, vPar, vPerp);
+
+ CHECK((helix.GetPosition(1_s).GetCoordinates() -
+ corsika::QuantityVector<length_d>(0_m, 0_m, 4_m))
+ .norm()
+ .magnitude() == Approx(0).margin(absMargin));
+
+ CHECK((helix.GetPosition(0.25_s).GetCoordinates() -
+ corsika::QuantityVector<length_d>(-1_m / (2 * M_PI), -1_m / (2 * M_PI), 1_m))
+ .norm()
+ .magnitude() == Approx(0).margin(absMargin));
+
+ BaseTrajectory const* base = &helix;
+ CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
+ base->GetPosition(1234_s).GetCoordinates());
+ }
}
diff --git a/Framework/Particles/ParticleClassNames.xml b/Framework/Particles/ParticleClassNames.xml
index 98772d67e38f568a91cf1dea6cc4dc91198e4bb4..2c909fb1ae904e9f22d3af0da1682ef842aee54b 100644
--- a/Framework/Particles/ParticleClassNames.xml
+++ b/Framework/Particles/ParticleClassNames.xml
@@ -3,20 +3,20 @@
<!-- For selected particles it is possible to specify the C++ class
names for a specific unique PDG code -->
- <particle id="0" classname="unknown"> </particle> <!-- VOID in Pythia8 -->
+ <particle pdgID="0" classname="unknown"/> <!-- VOID in Pythia8 -->
- <particle id="11" classname="Electron"> </particle>
- <particle id="-11" classname="Positron"> </particle>
- <particle id="22" classname="Gamma"> </particle>
+ <particle pdgID="11" classname="Electron"/>
+ <particle pdgID="-11" classname="Positron"/>
+ <particle pdgID="22" classname="Gamma"/>
- <particle id="130" classname="K0Long"> </particle>
- <particle id="310" classname="K0Short"> </particle>
+ <particle pdgID="130" classname="K0Long"/>
+ <particle pdgID="310" classname="K0Short"/>
- <particle id="2112" classname="Neutron"> </particle>
- <particle id="-2112" classname="AntiNeutron"> </particle>
+ <particle pdgID="2112" classname="Neutron"/>
+ <particle pdgID="-2112" classname="AntiNeutron"/>
- <particle id="2212" classname="Proton"> </particle>
- <particle id="-2212" classname="AntiProton"> </particle>
+ <particle pdgID="2212" classname="Proton"/>
+ <particle pdgID="-2212" classname="AntiProton"/>
</chapter>
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index 7df92b661d891cf252cbef52205e1818e71a834d..23ef2369886f8b102f973c89dd921c0aab731e28 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -4,15 +4,16 @@
Interface to particle properties
*/
-#ifndef _include_Particle_h_
-#define _include_Particle_h_
+#ifndef _include_ParticleProperties_h_
+#define _include_ParticleProperties_h_
#include <array>
#include <cstdint>
#include <iostream>
+#include <type_traits>
+#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/particles/GeneratedParticleProperties.inc>
/**
* @namespace particle
@@ -23,31 +24,50 @@
*/
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
- *
- * return mass of particle
+ /*!
+ * returns 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) {
- return GetElectricChargeNumber(p) * (corsika::units::constants::e);
+ corsika::units::ElectricChargeType constexpr GetElectricCharge(Code const p) {
+ 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 {
std::ostream& operator<<(std::ostream& stream, Code const p) {
- stream << GetName(p);
- return stream;
+ return stream << GetName(p);
}
} // namespace io
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index 46958ea055785689515d4cb6aa6a127dbe125a7d..4c271e8fbbcf4b8030b7f6df86f6ae8831681f83 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/Framework/Particles/pdxml_reader.py
@@ -51,7 +51,7 @@ def class_names(filename):
for particle in root.iter("particle"):
name = particle.attrib["classname"]
- pdg_id = int(particle.attrib["id"])
+ pdg_id = int(particle.attrib["pdgID"])
map[pdg_id] = name
return map
@@ -97,9 +97,10 @@ def c_identifier(name):
#
# This function produces names of type "DeltaPlusPlus"
#
-def c_identifier_cammel(name):
+def c_identifier_camel(name):
orig = name
name = name[0].upper() + name[1:].lower() # all lower case
+
for c in "() /": # replace funny characters
name = name.replace(c, "_")
@@ -111,9 +112,9 @@ def c_identifier_cammel(name):
# move "Bar" to end of name
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')
-
+
# cleanup "_"s
while True:
tmp = name.replace("__", "_")
@@ -127,18 +128,18 @@ def c_identifier_cammel(name):
istart = 0
while True:
i = name.find('_', istart)
- if (i<1 or i>len(name)-1):
+ if i < 1 or i > len(name)-1:
break
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
break
name = name[:i] + name[i+1:]
# and last, for example: make NuE out of Nue
- if (name[i-1].islower() and name[i].islower()) :
- if (i<len(name)-1) :
+ if name[i-1].islower() and name[i].islower():
+ if i < len(name)-1:
name = name[:i] + name[i].upper() + name[i+1:]
- else :
+ else:
name = name[:i] + name[i].upper()
# check if name is valid C++ identifier
@@ -146,27 +147,25 @@ def c_identifier_cammel(name):
if pattern.match(name):
return name
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
#
def build_pythia_db(filename, classnames):
- particle_db = OrderedDict()
+ particle_db = OrderedDict()
for (pdg, name, mass, electric_charge, antiName) in parse(filename):
c_id = "unknown"
- if (pdg in classnames):
+ if pdg in classnames:
c_id = classnames[pdg]
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] = {
"name" : name,
"antiName" : antiName,
@@ -186,17 +185,20 @@ def build_pythia_db(filename, classnames):
#
def gen_internal_enum(pythia_db):
- string = "enum class Code : uint8_t {\n"
-
- 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
+ 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...
+
+
for k in filter(lambda k: "ngc_code" in pythia_db[k], pythia_db):
last_ngc_id = pythia_db[k]['ngc_code']
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
@@ -214,13 +216,13 @@ def gen_properties(pythia_db):
string += "\n"
# 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():
- 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"
# 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():
string += " {pdg:d}, // {name:s}\n".format(pdg = p['pdg'], name = p['name'])
string += "};\n"
@@ -232,7 +234,7 @@ def gen_properties(pythia_db):
string += "};\n"
# 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():
string += " {charge:d},\n".format(charge = p['electric_charge'])
string += "};\n"
@@ -269,23 +271,23 @@ def gen_classes(pythia_db):
string += "/** @class " + cname + "\n\n"
string += " * Particle properties are taken from the PYTHIA8 ParticleData.xml file:<br>\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 += " * - name=" + str(cname) + "\n"
string += " * - anti=" + str(antiP) + "\n"
string += "*/\n\n"
- string += "class " + cname + "{\n"
+ string += "class " + cname + " {\n"
string += " public:\n"
- string += " static Code GetCode() { return Type; }\n"
- string += " static quantity<energy_d> GetMass() { return masses[TypeIndex]; }\n"
- string += " static quantity<electric_charge_d> GetCharge() { return corsika::units::constants::e*electric_charge[TypeIndex]/3; }\n"
- string += " static int GetChargeNumber() { return electric_charge[TypeIndex]/3; }\n"
- string += " static std::string GetName() { return names[TypeIndex]; }\n"
- string += " static Code GetAntiParticle() { return AntiType; }\n"
- string += " static const Code Type = Code::" + cname + ";\n"
- string += " static const Code AntiType = Code::" + antiP + ";\n"
+ string += " static constexpr Code GetCode() { return Type; }\n"
+ string += " static constexpr corsika::units::MassType GetMass() { return corsika::particles::GetMass(Type); }\n"
+ string += " static constexpr corsika::units::ElectricChargeType GetCharge() { return corsika::particles::GetElectricCharge(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 constexpr Code GetAntiParticle() { return AntiType; }\n"
+ string += " static constexpr Code Type = Code::" + cname + ";\n"
+ string += " static constexpr Code AntiType = Code::" + antiP + ";\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"
return string
@@ -296,22 +298,8 @@ def gen_classes(pythia_db):
#
def inc_start():
- string = ""
- string += "#ifndef _include_GeneratedParticleDataTable_h_\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"
+ string = ('// generated by pdxml_reader.py\n'
+ '// MANUAL EDITS ON OWN RISK\n')
return string
@@ -321,9 +309,6 @@ def inc_start():
def inc_end():
string = ""
- string += "\n}\n\n"
- string += "\n}\n\n"
- string += "#endif\n"
return string
@@ -334,14 +319,14 @@ def inc_end():
if __name__ == "__main__":
- if (len(sys.argv)!=3) :
- print ("pdxml_reader.py Pythia8.xml ClassNames.xml")
- sys.exit(0)
+ if len(sys.argv) != 3:
+ print("usage: {:s} <Pythia8.xml> <ClassNames.xml>".format(sys.argv[0]))
+ sys.exit(1)
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("\n pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n")
counter = itertools.count(0)
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index 476264171bb8be56f715f0eee3ee9fe74f77116b..a04f36515c720c6f55938238c5a42f685b527f99 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -14,7 +14,7 @@ TEST_CASE("Particles", "[Particles]") {
SECTION("Types") { REQUIRE(Electron::GetCode() == Code::Electron); }
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::GetCharge() / constants::e == Approx(-1));
REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h
index a666f15ef3bfc61d7abfb589a9763597db4f7d62..5c1b3a91c8bce590416b1889e1021a3e00dcfca9 100644
--- a/Framework/Units/PhysicalConstants.h
+++ b/Framework/Units/PhysicalConstants.h
@@ -39,16 +39,17 @@ namespace corsika::units::constants {
// Avogadro constant
constexpr quantity<dimensions<0, 0, 0, 0, 0, -1> > N_sub_A{Rep(6.02214199e+23L) / mole};
// electronvolt
- constexpr quantity<energy_d> eV{Rep(1.60217733e-19L) * joule};
+ constexpr quantity<energy_d> eV{Rep(1.6021766208e-19L) * joule};
// 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
constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second};
// speed of light in a vacuum
constexpr quantity<speed_d> c{Rep(299792458L) * meter / second};
+ constexpr auto cSquared = c * c;
// unified atomic mass unit
constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 91f8161637a8a84b34e3218e54ff4a9165f2ecf9..c1cf07bb9c05e4b74ff26e7ad4da48d1a22e5ca1 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -25,7 +25,6 @@ namespace phys {
} // namespace phys
namespace corsika::units {
-
using namespace phys::units;
using namespace phys::units::literals;
// namespace literals = phys::units::literals;
@@ -37,6 +36,7 @@ namespace corsika::units {
using ElectricChargeType =
phys::units::quantity<phys::units::electric_charge_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