towards integration of environment in cascade

Framework/Cascade/Cascade.h

@@ -19,15 +19,11 @@
 
 #include <type_traits>
 
-using namespace corsika;
-using namespace corsika::units::si;
-
 namespace corsika::cascade {
 
   template <typename Tracking, typename ProcessList, typename Stack>
   class Cascade {
-
-    typedef typename Stack::ParticleType Particle;
+    using Particle = typename Stack::ParticleType;
 
     Cascade() = delete;
 
@@ -57,39 +53,35 @@ namespace corsika::cascade {
 
     void Step(Particle& particle) {
 
-      // get access to random number generator
-      static corsika::random::RNG& rmng =
-          corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-
       // determine geometric tracking
       corsika::setup::Trajectory step = fTracking.GetTrack(particle);
 
       // determine combined total interaction length (inverse)
-      const double total_inv_lambda =
+      InverseLengthType const total_inv_lambda =
           fProcesseList.GetTotalInverseInteractionLength(particle, step);
 
       // sample random exponential step length
-      const double sample_step = rmng() / (double)rmng.max();
-      const double next_interact = -log(sample_step) / total_inv_lambda;
+      std::exponential_distribution expDist(total_inv_lambda * 1_m);
+      LengthType const next_interact = 1_m * expDist(fRNG);
       std::cout << "total_inv_lambda=" << total_inv_lambda
                 << ", next_interact=" << next_interact << std::endl;
 
       // determine the maximum geometric step length
-      const double distance_max = fProcesseList.MaxStepLength(particle, step);
+      LengthType const distance_max = fProcesseList.MaxStepLength(particle, step);
       std::cout << "distance_max=" << distance_max << std::endl;
 
       // determine combined total inverse decay time
-      const double total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
+      InverseTimeType const total_inv_lifetime = fProcesseList.GetTotalInverseLifetime(particle);
 
       // sample random exponential decay time
-      const double sample_decay = rmng() / (double)rmng.max();
-      const double next_decay = -log(sample_decay) / total_inv_lifetime;
+      std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s);
+      TimeType const next_decay = 1_s * expDistDecay(fRNG);
       std::cout << "total_inv_lifetime=" << total_inv_lifetime
                 << ", next_decay=" << next_decay << std::endl;
 
       // convert next_step from grammage to length [m]
       // Environment::GetDistance(step, next_step);
-      const double distance_interact = next_interact;
+      const GrammageType distance_interact = ;
       // ....
 
       // convert next_decay from time to length [m]
@@ -144,6 +136,10 @@ namespace corsika::cascade {
     Tracking& fTracking;
     ProcessList& fProcesseList;
     Stack& fStack;
+    corsika::environment::Environment const& fEnvironment;
+    corsika::random::RNG& fRNG =
+          corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
+
   };
 
 } // namespace corsika::cascade

Framework/Geometry/Helix.h

@@ -56,6 +56,10 @@ namespace corsika::geometry {
              (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
     }
 
+    Point PositionFromArclength(corsika::units::si::LengthType l) const {
+      return GetPosition(TimeFromArclength(l));
+    }
+
     auto GetRadius() const { return radius; }
 
     corsika::units::si::LengthType ArcLength(corsika::units::si::TimeType t1,
@@ -64,8 +68,8 @@ namespace corsika::geometry {
     }
 
     corsika::units::si::TimeType TimeFromArclength(
-        corsika::units::si::LengthType t) const {
-      return t / (vPar + vPerp).norm();
+        corsika::units::si::LengthType l) const {
+      return l / (vPar + vPerp).norm();
     }
   };
 

Framework/Geometry/Line.h

@@ -31,6 +31,8 @@ namespace corsika::geometry {
         , v0(pV0) {}
 
     Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
+    
+    Point PositionFromArclength(corsika::units::si::LengthType l) const { return r0 + v0.normalized() * l; }
 
     LengthType ArcLength(corsika::units::si::TimeType t1,
                          corsika::units::si::TimeType t2) const {

Framework/Geometry/testGeometry.cc

@@ -152,11 +152,21 @@ TEST_CASE("Trajectories") {
 
   SECTION("Line") {
     Vector<SpeedType::dimension_type> v0(rootCS,
-                                         {1_m / second, 0_m / second, 0_m / second});
+                                         {3_m / second, 0_m / second, 0_m / second});
 
     Line const line(r0, v0);
     CHECK(
-        (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(2_m, 0_m, 0_m))
+        (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(6_m, 0_m, 0_m))
+            .norm()
+            .magnitude() == Approx(0).margin(absMargin));
+
+    CHECK((line.PositionFromArclength(4_m).GetCoordinates() -
+           QuantityVector<length_d>(4_m, 0_m, 0_m))
+              .norm()
+              .magnitude() == Approx(0).margin(absMargin));
+              
+    CHECK(
+        (line.GetPosition(7_s) - line.PositionFromArclength(line.ArcLength(0_s, 7_s)))
             .norm()
             .magnitude() == Approx(0).margin(absMargin));
 
@@ -164,7 +174,7 @@ TEST_CASE("Trajectories") {
     Trajectory<Line> base(line, t);
     CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
 
-    CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(1));
+    CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3));
   }
 
   SECTION("Helix") {
@@ -174,7 +184,8 @@ TEST_CASE("Trajectories") {
     Vector<SpeedType::dimension_type> const vPerp(
         rootCS, {3_m / second, 0_m / second, 0_m / second});
 
-    auto const omegaC = 2 * M_PI / 1_s;
+    auto const T = 1_s;
+    auto const omegaC = 2 * M_PI / T;
 
     Helix const helix(r0, omegaC, vPar, vPerp);
 
@@ -188,10 +199,15 @@ TEST_CASE("Trajectories") {
               .norm()
               .magnitude() == Approx(0).margin(absMargin));
 
+    CHECK(
+        (helix.GetPosition(7_s) - helix.PositionFromArclength(helix.ArcLength(0_s, 7_s)))
+            .norm()
+            .magnitude() == Approx(0).margin(absMargin));
+
     auto const t = 1234_s;
     Trajectory<Helix> const base(helix, t);
     CHECK(helix.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
 
-    CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(5));
+    CHECK(base.ArcLength(0_s, 1_s) / 1_m == Approx(5));
   }
 }

Framework/Units/PhysicalUnits.h

@@ -39,8 +39,7 @@ namespace corsika::units::si {
   constexpr phys::units::quantity<momentum_d> newton_second{meter * kilogram / second};
 
   /// defining cross section
-  using sigma_d = phys::units::dimensions<2, 0, 0>;
-  constexpr phys::units::quantity<sigma_d> barn{Rep(1.e-28L) * meter * meter};
+  constexpr phys::units::quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
 
   /// add the unit-types
   using LengthType = phys::units::quantity<phys::units::length_d, double>;
@@ -54,7 +53,9 @@ namespace corsika::units::si {
   using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
   using GrammageType = phys::units::quantity<phys::units::dimensions<-2, 1, 0>, double>;
   using MomentumType = phys::units::quantity<momentum_d, double>;
-  using CrossSectionType = phys::units::quantity<sigma_d, double>;
+  using CrossSectionType = phys::units::quantity<area_d, double>;
+  using InverseLengthType = phys::units::quantity<phys::units::dimensions<-1, 0, 0>, double>;
+  using InverseTimeType = phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
 
 } // end namespace corsika::units::si
 
@@ -76,11 +77,8 @@ namespace phys {
       QUANTITY_DEFINE_SCALING_LITERALS(barn, corsika::units::si::sigma_d,
                                        magnitude(corsika::units::si::constants::barn))
 
-      QUANTITY_DEFINE_SCALING_LITERALS(meter, length_d,
-                                       magnitude(corsika::units::si::constants::meter))
-
-      QUANTITY_DEFINE_SCALING_LITERALS(newton_second, corsika::units::si::momentum_d,
-                                       magnitude(corsika::units::si::newton_second))
+      QUANTITY_DEFINE_SCALING_LITERALS(Ns, corsika::units::si::momentum_d,
+                                       magnitude(1_m * 1_kg / 1_s))
 
     } // namespace literals
   }   // namespace units

Processes/TrackingLine/TrackingLine.h

@@ -41,14 +41,14 @@ namespace corsika::process {
 
     public:
       std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
-      TimeOfIntersection(corsika::geometry::Line const& traj,
+      TimeOfIntersection(corsika::geometry::Line const& line,
                          geometry::Sphere const& sphere) {
         using namespace corsika::units::si;
         auto const& cs = fEnvironment.GetCoordinateSystem();
         geometry::Point const origin(cs, 0_m, 0_m, 0_m);
 
-        auto const r0 = (traj.GetR0() - origin);
-        auto const v0 = traj.GetV0();
+        auto const r0 = (line.GetR0() - origin);
+        auto const v0 = line.GetV0();
         auto const c0 = (sphere.GetCenter() - origin);
 
         auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
@@ -74,12 +74,15 @@ namespace corsika::process {
           : fEnvironment(pEnv) {}
       void Init() {}
 
-      auto GetTrack(Particle& p) {
+      auto GetTrack(Particle const& p) {
         using namespace corsika::units::si;
+        
         geometry::Vector<SpeedType::dimension_type> const velocity =
             p.GetMomentum() / p.GetEnergy() * corsika::units::si::constants::cSquared;
 
         auto const currentPosition = p.GetPosition();
+        
+        // to do: include effect of magnetic field
         geometry::Line line(currentPosition, velocity);
 
         auto const* currentVolumeNode =