IAP GITLAB

Skip to content
Snippets Groups Projects
Commit b4819cae authored by ralfulrich's avatar ralfulrich
Browse files

merged tranjectory stuff

parents 67e78cad b4b2ffb5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 341 additions and 122 deletions
CMakeLists.txt
+ 1
− 1
View file @ b4819cae
...@@ -14,7 +14,7 @@ set (CMAKE_INSTALL_MESSAGE LAZY) ...@@ -14,7 +14,7 @@ set (CMAKE_INSTALL_MESSAGE LAZY)
set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/CMakeModules) set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/CMakeModules)
include (CorsikaUtilities) # a few cmake function include (CorsikaUtilities) # a few cmake function
# --std=c++14 # --std=c++17
set (CMAKE_CXX_STANDARD 17) set (CMAKE_CXX_STANDARD 17)
enable_testing () enable_testing ()
set (CTEST_OUTPUT_ON_FAILURE 1) set (CTEST_OUTPUT_ON_FAILURE 1)
......
Framework/Geometry/BaseTrajectory.h 0 → 100644
+ 23
− 0
View file @ b4819cae
#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
Framework/Geometry/CMakeLists.txt
+ 2
− 0
View file @ b4819cae
...@@ -13,7 +13,9 @@ set ( ...@@ -13,7 +13,9 @@ set (
Helix.h Helix.h
BaseVector.h BaseVector.h
QuantityVector.h QuantityVector.h
BaseTrajectory.h
LineTrajectory.h LineTrajectory.h
Trajectory.h
) )
set ( set (
......
Framework/Geometry/Helix.h
+ 23
− 16
View file @ b4819cae
#ifndef _include_HELIX_H_ #ifndef _include_HELIX_H_
#define _include_HELIX_H_ #define _include_HELIX_H_
#include <corsika/geometry/BaseTrajectory.h>
#include <corsika/geometry/Point.h> #include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h> #include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h> #include <corsika/units/PhysicalUnits.h>
#include <cmath> #include <cmath>
namespace corsika::geometry { 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; Point const r0;
FrequencyType const omegaC; corsika::units::FrequencyType const omegaC;
SpeedVec const vPar; VelocityVec const vPar;
SpeedVec vPerp, uPerp; VelocityVec const vPerp, uPerp;
LengthType const radius;
corsika::units::LengthType const radius;
public: public:
Helix(Point const& pR0, quantity<frequency_d> pOmegaC, SpeedVec const& pvPar, Helix(Point const& pR0, corsika::units::FrequencyType pOmegaC,
SpeedVec const& pvPerp) VelocityVec const& pvPar, VelocityVec const& pvPerp)
: r0(pR0) : r0(pR0)
, omegaC(pOmegaC) , omegaC(pOmegaC)
, vPar(pvPar) , vPar(pvPar)
...@@ -33,7 +40,7 @@ namespace corsika::geometry { ...@@ -33,7 +40,7 @@ namespace corsika::geometry {
, uPerp(vPerp.cross(vPar.normalized())) , uPerp(vPerp.cross(vPar.normalized()))
, radius(pvPar.norm() / abs(pOmegaC)) {} , radius(pvPar.norm() / abs(pOmegaC)) {}
auto GetPosition(TimeType t) const { Point GetPosition(corsika::units::TimeType t) const {
return r0 + vPar * t + return r0 + vPar * t +
(vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC; (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
} }
......
Framework/Geometry/LineTrajectory.h
+ 7
− 7
View file @ b4819cae
#ifndef _include_LINETRAJECTORY_H #ifndef _include_LINETRAJECTORY_H
#define _include_LINETRAJECTORY_H #define _include_LINETRAJECTORY_H
#include <corsika/geometry/BaseTrajectory.h>
#include <corsika/geometry/Point.h> #include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h> #include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h> #include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry { namespace corsika::geometry {
class LineTrajectory // TODO: inherit from Trajectory class LineTrajectory : public BaseTrajectory {
{ using VelocityVec = Vector<corsika::units::SpeedType::dimension_type>;
using SpeedVec = Vector<corsika::units::SpeedType::dimension_type>;
Point const r0; Point const r0;
SpeedVec const v0; VelocityVec const v0;
public: public:
LineTrajectory(Point const& pR0, SpeedVec const& pV0) LineTrajectory(Point const& pR0, VelocityVec const& pV0)
: r0(r0) : r0(pR0)
, v0(pV0) {} , 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 } // namespace corsika::geometry
......
Framework/Geometry/QuantityVector.h
+ 2
− 0
View file @ b4819cae
...@@ -101,6 +101,8 @@ namespace corsika { ...@@ -101,6 +101,8 @@ namespace corsika {
auto& operator-() const { return QuantityVector<dim>(-eVector); } auto& operator-() const { return QuantityVector<dim>(-eVector); }
auto normalized() const { return (*this) * (1 / norm()); } auto normalized() const { return (*this) * (1 / norm()); }
auto operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
}; };
} // end namespace corsika } // end namespace corsika
......
Framework/Geometry/Sphere.h
+ 1
− 0
View file @ b4819cae
...@@ -15,6 +15,7 @@ namespace corsika::geometry { ...@@ -15,6 +15,7 @@ namespace corsika::geometry {
: center(pCenter) : center(pCenter)
, radius(pRadius) {} , radius(pRadius) {}
//! returns true if the Point p is within the sphere
auto isInside(Point const& p) const { auto isInside(Point const& p) const {
return radius * radius > (center - p).squaredNorm(); return radius * radius > (center - p).squaredNorm();
} }
......
Framework/Geometry/Trajectory.h 0 → 100644
+ 31
− 0
View file @ b4819cae
#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
Framework/Geometry/testGeometry.cc
+ 155
− 8
View file @ b4819cae
...@@ -2,19 +2,166 @@ ...@@ -2,19 +2,166 @@
// cpp file // cpp file
#include <catch2/catch.hpp> #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 <corsika/units/PhysicalUnits.h>
#include <cmath>
using namespace phys::units; using namespace corsika::geometry;
using namespace phys::units::literals; using namespace corsika::units;
TEST_CASE("Geometry", "[Geometry]") { double constexpr absMargin = 1.0e-8;
SECTION("Coordinate Systems") {} ;
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());
}
} }
Framework/Particles/ParticleClassNames.xml
+ 10
− 10
View file @ b4819cae
...@@ -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 id="0" classname="unknown"> </particle> <!-- VOID in Pythia8 --> <particle pdgID="0" classname="unknown"/> <!-- VOID in Pythia8 -->
<particle id="11" classname="Electron"> </particle> <particle pdgID="11" classname="Electron"/>
<particle id="-11" classname="Positron"> </particle> <particle pdgID="-11" classname="Positron"/>
<particle id="22" classname="Gamma"> </particle> <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 @ b4819cae
...@@ -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 @ b4819cae
...@@ -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 @ b4819cae
...@@ -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 @ b4819cae
...@@ -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 @ b4819cae
...@@ -25,7 +25,6 @@ namespace phys { ...@@ -25,7 +25,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;
...@@ -37,6 +36,7 @@ namespace corsika::units { ...@@ -37,6 +36,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