diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 00a152fffb194ca1599ce5ce7180feb099618ef7..cfd96473717afafb31fd4cc52fe266fe48fef024 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -1,3 +1,6 @@
+add_executable (helix_example helix_example.cc)
+target_link_libraries(helix_example CORSIKAgeometry CORSIKAunits)
+install (TARGETS helix_example DESTINATION share/examples)
add_executable (geometry_example geometry_example.cc)
target_link_libraries (geometry_example CORSIKAgeometry CORSIKAunits)
diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc
index 742de1202d9a02c5e5c7b2889d46cae72793fbab..6927bd8380e3968a4f19c531e61474324ea850a7 100644
--- a/Documentation/Examples/geometry_example.cc
+++ b/Documentation/Examples/geometry_example.cc
@@ -14,7 +14,7 @@ using namespace phys::units::literals; // support unit literals like 5_m;
int main()
{
- //~ // define the root coordinate system
+ // define the root coordinate system
CoordinateSystem root;
// another CS defined by a translation relative to the root CS
@@ -37,10 +37,10 @@ int main()
std::cout << "p2-p1 norm^2: " << norm << std::endl;
Sphere s(p1, 10_m); // define a sphere around a point with a radius
- //~ std::cout << "p1 inside s: " << s.isInside(p2) << std::endl;
+ std::cout << "p1 inside s: " << s.isInside(p2) << std::endl;
Sphere s2(p1, 3_um); // another sphere
- //~ std::cout << "p1 inside s2: " << s2.isInside(p2) << std::endl;
+ std::cout << "p1 inside s2: " << s2.isInside(p2) << std::endl;
// let's try parallel projections:
@@ -49,14 +49,15 @@ int main()
auto const v3 = v1.parallelProjectionOnto(v2);
- auto const cross = v1.cross(v2).normalized();
+ // cross product
+ auto const cross = v1.cross(v2).normalized(); // normalized() returns dimensionless, normalized vectors
// if a CS is not given as parameter for getComponents(), the components
// in the "home" CS are returned
- std::cout << v1.getComponents() << std::endl;
- std::cout << v2.getComponents() << std::endl;
- std::cout << v3.getComponents() << std::endl;
- std::cout << cross.getComponents() << std::endl;
+ std::cout << "v1: " << v1.getComponents() << std::endl;
+ std::cout << "v2: " <<v2.getComponents() << std::endl;
+ std::cout << "parallel projection of v1 onto v2: " << v3.getComponents() << std::endl;
+ std::cout << "normalized cross product of v1 x v2" << cross.getComponents() << std::endl;
return EXIT_SUCCESS;
}
diff --git a/Documentation/Examples/helix_example.cc b/Documentation/Examples/helix_example.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ebb5cb4e9d2e4da785d76d960496f9d2074d214e
--- /dev/null
+++ b/Documentation/Examples/helix_example.cc
@@ -0,0 +1,46 @@
+#include <Units/PhysicalUnits.h>
+#include <Geometry/Vector.h>
+#include <Geometry/CoordinateSystem.h>
+#include <Geometry/Helix.h>
+#include <Geometry/Point.h>
+#include <cstdlib>
+#include <iostream>
+#include <array>
+
+using namespace phys::units;
+using namespace phys::units::io; // support stream << unit
+using namespace phys::units::literals; // support unit literals like 5_m;
+
+int main()
+{
+ CoordinateSystem root;
+
+ Point const r0(root, {0_m, 0_m, 0_m});
+ auto const omegaC = 2 * M_PI * 1_Hz;
+ Vector<speed_d> vPar(root, {0_m / second, 0_m / second, 10_cm / second});
+ Vector<speed_d> vPerp(root, {1_m / second, 0_m / second, 0_m / second});
+
+ Helix h(r0, omegaC, vPar, vPerp);
+
+ auto constexpr t0 = 0_s;
+ auto constexpr t1 = 1_s;
+ auto constexpr dt = 1_us;
+ auto constexpr n = long((t1 - t0) / dt) + 1;
+
+ auto arr = std::make_unique<std::array<std::array<double, 4>, n>>();
+ auto& positions = *arr;
+
+ for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i)
+ {
+ auto const r = h.getPosition(t).getCoordinates();
+
+ positions[i][0] = t / 1_s;
+ positions[i][1] = r[0] / 1_m;
+ positions[i][2] = r[1] / 1_m;
+ positions[i][3] = r[2] / 1_m;
+ }
+
+ std::cout << positions[n-2][0] << " " << positions[n-2][1] << " " << positions[n-2][2] << " " << positions[n-2][3] << std::endl;
+
+ return EXIT_SUCCESS;
+}
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 9dacfb9dd6920d7ffe71c2e0d982dce23339d941..a1ba107519a5e485f2644fc3501e3d9314f12c75 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -7,25 +7,33 @@
class Helix // TODO: inherit from to-be-implemented "Trajectory"
{
- using SpeedVec = Vector<phys::units::speed_d>;
+ using SpeedVec = Vector<Speed>;
Point const r0;
- phys::units::quantity<phys::units::frequency_d> const omegaC;
+ Frequency const omegaC;
SpeedVec const vPar;
SpeedVec vPerp, uPerp;
+ Length const radius;
+
public:
Helix(Point const pR0, phys::units::quantity<phys::units::frequency_d> pOmegaC,
SpeedVec const pvPar, SpeedVec const pvPerp) :
- r0(pR0), omegaC(pOmegaC), vPar(pvPar), vPerp(pvPerp), uPerp(vPar.normalized().cross(vPerp))
+ r0(pR0), omegaC(pOmegaC), vPar(pvPar), vPerp(pvPerp), uPerp(vPar.normalized().cross(vPerp)),
+ radius(pvPar.norm() / pOmegaC)
{
}
- Point getPosition(phys::units::quantity<phys::units::time_interval_d> t) const
+ auto getPosition(Time t) const
{
return r0 + vPar * t + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
}
+
+ auto getRadius() const
+ {
+ return radius;
+ }
};
#endif
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index 8176715122102302637e332f24c519325c3c29f0..26e6c2a4104dfe60cca55b48edfd9792b30c9810 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -86,7 +86,7 @@ public:
return QuantityVector<dim>(eVector / p);
}
- auto& operator/=(double const p) const
+ auto& operator/=(double const p)
{
eVector /= p;
return *this;
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 9c21910cf188ba58f43bb17b6288811af4a94ee1..cc87e1088161c0ea6800bad116c154e23deff015 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -6,7 +6,6 @@
class Sphere
{
- using Length = phys::units::quantity<phys::units::length_d, double>;
Point center;
Length const radius;
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index b85d35327cf2796f1d690d3a33c5cc5e94c9a19c..648e68f00afec5c462ef9d70da0435e834f7c4d7 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -4,6 +4,7 @@
#include <Units/PhysicalUnits.h>
using namespace phys::units;
+using namespace phys::units::literals;
TEST_CASE( "PhysicalUnits", "[Units]" )
{
diff --git a/Framework/Logging/Makefile b/Framework/Logging/Makefile
deleted file mode 100644
index c6c84cac0c0fd96166af6eb86fa9e390a2ad3811..0000000000000000000000000000000000000000
--- a/Framework/Logging/Makefile
+++ /dev/null
@@ -1,10 +0,0 @@
-CXXFLAGS+=-I. --std=c++14 -O3
-
-all: test test_off
-
-#test: test.o
-# $(CXX) $(CXXFLAGS) $(LDFLAGS) $^ -o $@
-
-clean:
- rm -rf *.o test test_off *.dat *.log
-
diff --git a/Framework/ParticleProperties/ParticleProperties.hpp b/Framework/ParticleProperties/ParticleProperties.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f627966915e96af83dd282bc75cf16e65ca1b280
--- /dev/null
+++ b/Framework/ParticleProperties/ParticleProperties.hpp
@@ -0,0 +1,45 @@
+#ifndef PARTICLEPROPERTIES_HPP
+#define PARTICLEPROPERTIES_HPP
+#include <array>
+#include <cstdint>
+#include <iostream>
+
+namespace ParticleProperties {
+
+typedef int16_t PDGCode;
+
+#include "generated_particle_properties.inc"
+
+auto constexpr getMass(InternalParticleCode const p)
+{
+ return masses[static_cast<uint8_t const>(p)];
+}
+
+auto constexpr getPDG(InternalParticleCode const p)
+{
+ return pdg_codes[static_cast<uint8_t const>(p)];
+}
+
+auto constexpr getElectricChargeQN(InternalParticleCode const p)
+{
+ return electric_charge[static_cast<uint8_t const>(p)];
+}
+
+auto constexpr getElectricCharge(InternalParticleCode const p)
+{
+ return getElectricChargeQN(p) * (phys::units::e / 3.);
+}
+
+auto const getName(InternalParticleCode const p)
+{
+ return names[static_cast<uint8_t const>(p)];
+}
+
+std::ostream& operator<< (std::ostream& stream, InternalParticleCode const p)
+{
+ stream << getName(p);
+ return stream;
+}
+
+}
+#endif
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 83ed903f9f2ffd7da290349f69daa059aa8d18e8..d9bb61ac4d78826cc4a5c2bfdbc2c90ee661e89e 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -11,8 +11,8 @@
Define _XeV literals, alowing 10_GeV in the code.
*/
-using namespace phys::units::io;
-using namespace phys::units::literals;
+/*using namespace phys::units::io;
+using namespace phys::units::literals;*/
namespace phys {
namespace units {
@@ -22,5 +22,11 @@ namespace phys {
}
}
+using Length = phys::units::quantity<phys::units::length_d, double>;
+using Time = phys::units::quantity<phys::units::time_interval_d, double>;
+using Speed = phys::units::quantity<phys::units::speed_d, double>;
+using Frequency = phys::units::quantity<phys::units::frequency_d, double>;
+using ElectricCharge = phys::units::quantity<phys::units::electric_charge_d, double>;
+
#endif
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index 3206c50951328502350614c244fd5e0e7e35c3b6..305d06fbd0093b32a341fe093be8ee47aae3ab52 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -4,6 +4,7 @@
#include <Units/PhysicalUnits.h>
using namespace phys::units;
+using namespace phys::units::literals;
TEST_CASE( "PhysicalUnits", "[Units]" ) {
REQUIRE( 1_m/1_m == 1 );
diff --git a/ThirdParty/phys/units/quantity.hpp b/ThirdParty/phys/units/quantity.hpp
index 29ab27f65e487a088a71e0ac068b539dbf608f80..5f764409d6939ea43c393e254f2f62208c192f87 100644
--- a/ThirdParty/phys/units/quantity.hpp
+++ b/ThirdParty/phys/units/quantity.hpp
@@ -340,7 +340,7 @@ private:
enum { has_dimension = ! Dims::is_all_zero };
- // static_assert( has_dimension, "quantity dimensions must not all be zero" );
+ // static_assert( has_dimension, "quantity dimensions must not all be zero" ); // MR: removed
private:
// friends:
diff --git a/Tools/pdxml_reader.py b/Tools/pdxml_reader.py
index b17c7b4903f8f6acc55ea1db66db188f76720d6e..fe3d06e974b00219be7afdff4dfc1cb927cc3f20 100755
--- a/Tools/pdxml_reader.py
+++ b/Tools/pdxml_reader.py
@@ -42,13 +42,13 @@ def parse(filename):
decay_width = float(particle.attrib.get("mWidth", 0)) # GeV
lifetime = float(particle.attrib.get("tau0", math.inf)) # mm / c
- yield (pdg_id, name, mass)
+ yield (pdg_id, name, mass, electric_charge)
# TODO: read decay channels from child elements
if "antiName" in particle.attrib:
name = particle.attrib['antiName']
- yield (-pdg_id, name, mass)
+ yield (-pdg_id, name, mass, -electric_charge)
def c_identifier(name):
@@ -82,7 +82,7 @@ def c_identifier(name):
def build_pythia_db(filename):
particle_db = OrderedDict()
- for (pdg, name, mass) in parse(filename):
+ for (pdg, name, mass, electric_charge) in parse(filename):
c_id = c_identifier(name)
#~ print(name, c_id, sep='\t', file=sys.stderr)
@@ -91,6 +91,7 @@ def build_pythia_db(filename):
"name" : name,
"pdg" : pdg,
"mass" : mass, # in GeV
+ "electric_charge" : electric_charge # in e/3
}
return particle_db
@@ -163,6 +164,8 @@ def gen_internal_enum(pythia_db):
def gen_properties(pythia_db):
+
+ # masses
string = "static constexpr std::size_t size = {size:d};\n".format(size = len(pythia_db))
string += "static constexpr std::array<double const, size> masses{{\n"
@@ -171,24 +174,34 @@ def gen_properties(pythia_db):
string += " {mass:f}, // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += ("}};\n"
+
+ # PDG codes
"static constexpr std::array<PDGCode 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"
+
+ # name strings
"static const std::array<std::string const, size> names{{\n")
for p in pythia_db.values():
string += " \"{name:s}\",\n".format(name = p['name'])
- string += "}};\n"
+ string += ("}};\n"
+
+ # electric charges
+ "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'])
return string
if __name__ == "__main__":
- pythia_db = build_pythia_db("ParticleData.xml")
+ pythia_db = build_pythia_db("Tools/ParticleData.xml")
for c_id, sib_info in read_sibyll("sibyll_codes.dat"):
#~ print(c_id, sib_info)
@@ -216,7 +229,7 @@ if __name__ == "__main__":
#~ if table != pdg:
#~ raise Exception(p, sib_db, pdg, table)
- with open("generated_particle_properties.inc", "w") as f:
+ with open("Framework/ParticleProperties/generated_particle_properties.inc", "w") as f:
print(gen_internal_enum(pythia_db), file=f)
print(gen_properties(pythia_db), file=f)