huge commit before IKP long shutdown

Documentation/Examples/CMakeLists.txt

+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)

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;
 }

Documentation/Examples/helix_example.cc

+#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;
+}

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

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;

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;
     

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]" )
 {  

Framework/Logging/Makefile

-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
-

Framework/ParticleProperties/ParticleProperties.hpp

+#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

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
 

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 );

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:

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)