Merge branch '97-cascade-example-debugging-2' into 'master'

Resolve "cascade example debugging"

See merge request !39

.gitlab-ci.yml

@@ -16,7 +16,7 @@ build:
     - cd build
     - cmake ..
     - cmake --build .
-    - ctest -V
+    - ctest
 
 # code_quality:
 #   image: docker:stable

Documentation/Examples/cascade_example.cc

@@ -45,24 +45,24 @@ using namespace corsika::environment;
 using namespace std;
 using namespace corsika::units::hep;
 
-static EnergyType fEnergy = 0. * 1_GeV;
+class ProcessCut : public corsika::process::ContinuousProcess<ProcessCut> {
 
-// FOR NOW: global static variables for ParticleCut process
-// this is just wrong...
-static EnergyType fEmEnergy;
-static int fEmCount;
+  EnergyType fECut;
 
-static EnergyType fInvEnergy;
-static int fInvCount;
+  mutable EnergyType fEnergy = 0_GeV;
+  mutable EnergyType fEmEnergy = 0_GeV;
+  mutable int fEmCount = 0;
+  mutable EnergyType fInvEnergy = 0_GeV;
+  mutable int fInvCount = 0;
 
-class ProcessEMCut : public corsika::process::ContinuousProcess<ProcessEMCut> {
 public:
-  ProcessEMCut() {}
+  ProcessCut(const EnergyType v)
+      : fECut(v) {}
   template <typename Particle>
   bool isBelowEnergyCut(Particle& p) const {
     // FOR NOW: center-of-mass energy hard coded
     const EnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
-    if (p.GetEnergy() < 50_GeV || Ecm < 10_GeV)
+    if (p.GetEnergy() < fECut || Ecm < 10_GeV)
       return true;
     else
       return false;
@@ -134,25 +134,27 @@ public:
 
   template <typename Particle, typename Stack>
   EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
-    cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
     const Code pid = p.GetPID();
+    EnergyType energy = p.GetEnergy();
+    cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
+         << ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
     EProcessReturn ret = EProcessReturn::eOk;
     if (isEmParticle(pid)) {
       cout << "removing em. particle..." << endl;
-      fEmEnergy += p.GetEnergy();
+      fEmEnergy += energy;
       fEmCount += 1;
-      p.Delete();
+      // p.Delete();
       ret = EProcessReturn::eParticleAbsorbed;
     } else if (isInvisible(pid)) {
       cout << "removing inv. particle..." << endl;
-      fInvEnergy += p.GetEnergy();
+      fInvEnergy += energy;
       fInvCount += 1;
-      p.Delete();
+      // p.Delete();
       ret = EProcessReturn::eParticleAbsorbed;
     } else if (isBelowEnergyCut(p)) {
       cout << "removing low en. particle..." << endl;
-      fEnergy += p.GetEnergy();
-      p.Delete();
+      fEnergy += energy;
+      // p.Delete();
       ret = EProcessReturn::eParticleAbsorbed;
     }
     return ret;
@@ -178,20 +180,21 @@ public:
          << " ******************************" << endl;
   }
 
-  EnergyType GetInvEnergy() { return fInvEnergy; }
-
-  EnergyType GetCutEnergy() { return fEnergy; }
-
-  EnergyType GetEmEnergy() { return fEmEnergy; }
-
-private:
+  EnergyType GetInvEnergy() const { return fInvEnergy; }
+  EnergyType GetCutEnergy() const { return fEnergy; }
+  EnergyType GetEmEnergy() const { return fEmEnergy; }
 };
 
+//
+// The example main program for a particle cascade
+//
 int main() {
 
+  // initialize random number sequence(s)
   corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
 
-  corsika::environment::Environment env; // dummy environment
+  // setup environment, geometry
+  corsika::environment::Environment env;
   auto& universe = *(env.GetUniverse());
 
   auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
@@ -201,47 +204,47 @@ int main() {
   using MyHomogeneousModel =
       corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
   theMedium->SetModelProperties<MyHomogeneousModel>(
-      1_g / (1_m * 1_m * 1_m),
+      1_kg / (1_m * 1_m * 1_m),
       corsika::environment::NuclearComposition(
-          std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
+          std::vector<corsika::particles::Code>{corsika::particles::Code::Oxygen},
           std::vector<float>{1.}));
 
   universe.AddChild(std::move(theMedium));
 
-  CoordinateSystem& rootCS =
-      RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+  const CoordinateSystem& rootCS = env.GetCoordinateSystem();
 
+  // setup processes, decays and interactions
   tracking_line::TrackingLine<setup::Stack> tracking(env);
   stack_inspector::StackInspector<setup::Stack> p0(true);
-
   corsika::process::sibyll::Interaction sibyll;
   corsika::process::sibyll::Decay decay;
-  ProcessEMCut cut;
+  ProcessCut cut(8_GeV);
+
+  // assemble all processes into an ordered process list
   const auto sequence = /*p0 <<*/ sibyll << decay << cut;
 
+  // setup particle stack, and add primary particle
   setup::Stack stack;
+  stack.Clear();
+  const hep::EnergyType E0 = 1_TeV;
+  {
+    auto particle = stack.NewParticle();
+    particle.SetPID(Code::Proton);
+    hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
+    auto plab = stack::super_stupid::MomentumVector(rootCS, 0_GeV, 0_GeV, -P0);
+    particle.SetEnergy(E0);
+    particle.SetMomentum(plab);
+    particle.SetTime(0_ns);
+    Point p(rootCS, 0_m, 0_m, 10_km);
+    particle.SetPosition(p);
+  }
 
+  // define air shower object, run simulation
   corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
-
-  stack.Clear();
-  auto particle = stack.NewParticle();
-  EnergyType E0 = 100_GeV;
-  hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
-  auto plab = stack::super_stupid::MomentumVector(rootCS, 0. * 1_GeV, 0. * 1_GeV, P0);
-  particle.SetEnergy(E0);
-  particle.SetMomentum(plab);
-  particle.SetPID(Code::Proton);
-  particle.SetTime(0_ns);
-  Point p(rootCS, 0_m, 0_m, 0_m);
-  particle.SetPosition(p);
   EAS.Init();
   EAS.Run();
-  cout << "Result: E0="
-       << E0 / 1_GeV
-       //<< "GeV, particles below energy threshold =" << p1.GetCount()
-       << endl;
-  cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
-       << std::endl;
+
+  cout << "Result: E0=" << E0 / 1_GeV << endl;
   cut.ShowResults();
   cout << "total energy (GeV): "
        << (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy()) / 1_GeV << endl;

Framework/Cascade/Cascade.h

@@ -71,9 +71,8 @@ namespace corsika::cascade {
           fProcessSequence.GetTotalInverseInteractionLength(particle, step);
 
       // sample random exponential step length in grammage
-      auto constexpr grammageConversion = 1_g / (1_m * 1_m);
-      std::exponential_distribution expDist(1 / (grammageConversion * total_inv_lambda));
-      GrammageType const next_interact = grammageConversion * expDist(fRNG);
+      std::exponential_distribution expDist(total_inv_lambda * (1_g / (1_m * 1_m)));
+      GrammageType const next_interact = (1_g / (1_m * 1_m)) * expDist(fRNG);
 
       std::cout << "total_inv_lambda=" << total_inv_lambda
                 << ", next_interact=" << next_interact << std::endl;
@@ -100,7 +99,6 @@ namespace corsika::cascade {
                 << ", next_decay=" << next_decay << std::endl;
 
       // convert next_decay from time to length [m]
-      // Environment::GetDistance(step, next_decay);
       LengthType const distance_decay = next_decay * particle.GetMomentum().norm() /
                                         particle.GetEnergy() *
                                         corsika::units::constants::c;
@@ -109,10 +107,11 @@ namespace corsika::cascade {
       auto const min_distance =
           std::min({distance_interact, distance_decay, distance_max});
 
+      std::cout << " move particle by : " << min_distance << std::endl;
+
       // here the particle is actually moved along the trajectory to new position:
       // std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
       particle.SetPosition(step.PositionFromArclength(min_distance));
-
       // .... also update time, momentum, direction, ...
 
       // apply all continuous processes on particle + track
@@ -120,39 +119,40 @@ namespace corsika::cascade {
           fProcessSequence.DoContinuous(particle, step, fStack);
 
       if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
-        // fStack.Delete(particle); // TODO: check if this is really needed
-      } else {
-
-        std::cout << "sth. happening before geometric limit ?"
-                  << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
-
-        if (min_distance < distance_max) { // interaction to happen within geometric limit
-          // check weather decay or interaction limits this step
-
-          if (min_distance == distance_interact) {
-            std::cout << "collide" << std::endl;
-
-            InverseGrammageType const actual_inv_length =
-                fProcessSequence.GetTotalInverseInteractionLength(particle, step);
-
-            corsika::random::UniformRealDistribution<InverseGrammageType> uniDist(
-                actual_inv_length);
-            const auto sample_process = uniDist(fRNG);
-            InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
-            fProcessSequence.SelectInteraction(particle, fStack, sample_process,
-                                               inv_lambda_count);
-          } else {
-            std::cout << "decay" << std::endl;
-            InverseTimeType const actual_decay_time =
-                fProcessSequence.GetTotalInverseLifetime(particle);
-
-            corsika::random::UniformRealDistribution<InverseTimeType> uniDist(
-                actual_decay_time);
-            const auto sample_process = uniDist(fRNG);
-            InverseTimeType inv_decay_count = 0 / second;
-            fProcessSequence.SelectDecay(particle, fStack, sample_process,
-                                         inv_decay_count);
-          }
+        std::cout << "Cascade: delete absorbed particle " << particle.GetPID() << " "
+                  << particle.GetEnergy() / 1_GeV << "GeV" << std::endl;
+        particle.Delete();
+        return;
+      }
+
+      std::cout << "sth. happening before geometric limit ?"
+                << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
+
+      if (min_distance < distance_max) { // interaction to happen within geometric limit
+        // check weather decay or interaction limits this step
+
+        if (min_distance == distance_interact) {
+          std::cout << "collide" << std::endl;
+
+          InverseGrammageType const actual_inv_length =
+              fProcessSequence.GetTotalInverseInteractionLength(particle, step);
+
+          corsika::random::UniformRealDistribution<InverseGrammageType> uniDist(
+              actual_inv_length);
+          const auto sample_process = uniDist(fRNG);
+          InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
+          fProcessSequence.SelectInteraction(particle, fStack, sample_process,
+                                             inv_lambda_count);
+        } else {
+          std::cout << "decay" << std::endl;
+          InverseTimeType const actual_decay_time =
+              fProcessSequence.GetTotalInverseLifetime(particle);
+
+          corsika::random::UniformRealDistribution<InverseTimeType> uniDist(
+              actual_decay_time);
+          const auto sample_process = uniDist(fRNG);
+          InverseTimeType inv_decay_count = 0 / second;
+          fProcessSequence.SelectDecay(particle, fStack, sample_process, inv_decay_count);
         }
       }
     }

Framework/Cascade/testCascade.cc

@@ -80,7 +80,7 @@ public:
   }
 
   template <typename Particle, typename T, typename Stack>
-  void DoContinuous(Particle& p, T&, Stack& s) const {
+  EProcessReturn DoContinuous(Particle& p, T&, Stack& s) const {
     EnergyType E = p.GetEnergy();
     if (E < 85_MeV) {
       p.Delete();
@@ -95,6 +95,7 @@ public:
       pnew.SetPosition(p.GetPosition());
       pnew.SetMomentum(p.GetMomentum());
     }
+    return EProcessReturn::eOk;
   }
 
   void Init() { fCount = 0; }

Framework/ProcessSequence/ContinuousProcess.h

@@ -35,7 +35,7 @@ namespace corsika::process {
     // -> enforce derived to implement DoContinuous...
     template <typename Particle, typename Track, typename Stack>
     EProcessReturn DoContinuous(Particle&, Track&, Stack&) const;
-    
+
     // -> enforce derived to implement MaxStepLength...
     template <typename Particle, typename Track>
     corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const;

Framework/ProcessSequence/ProcessReturn.h

@@ -20,13 +20,25 @@ namespace corsika::process {
      that can be accumulated easily with "|="
    */
 
-  enum class EProcessReturn {
-    eOk = 1,
-    eParticleAbsorbed = 2,
-    eInteracted = 3,
-    eDecayed = 4,
+  enum class EProcessReturn : int {
+    eOk = (1 << 0),
+    eParticleAbsorbed = (1 << 2),
+    eInteracted = (1 << 3),
+    eDecayed = (1 << 4),
   };
 
+  inline EProcessReturn operator|(EProcessReturn a, EProcessReturn b) {
+    return static_cast<EProcessReturn>(static_cast<int>(a) | static_cast<int>(b));
+  }
+
+  inline EProcessReturn& operator|=(EProcessReturn& a, EProcessReturn b) {
+    return a = a | b;
+  }
+
+  inline bool operator==(EProcessReturn a, EProcessReturn b) {
+    return (static_cast<int>(a) & static_cast<int>(b)) != 0;
+  }
+
 } // namespace corsika::process
 
 #endif

Framework/ProcessSequence/ProcessSequence.h

@@ -60,20 +60,21 @@ namespace corsika::process {
       EProcessReturn ret = EProcessReturn::eOk;
       if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
                     is_process_sequence<T1>::value) {
-        A.DoContinuous(p, t, s);
+        ret |= A.DoContinuous(p, t, s);
       }
       if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
                     is_process_sequence<T2>::value) {
-        B.DoContinuous(p, t, s);
+        ret |= B.DoContinuous(p, t, s);
       }
       return ret;
     }
 
     template <typename Particle, typename Track>
     corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const {
-      corsika::units::si::LengthType max_length = // if no other process in the sequence implements it
+      corsika::units::si::LengthType
+          max_length = // if no other process in the sequence implements it
           std::numeric_limits<double>::infinity() * corsika::units::si::meter;
-          
+
       if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
                     is_process_sequence<T1>::value) {
         corsika::units::si::LengthType const len = A.MaxStepLength(p, track);

Framework/Units/PhysicalUnits.h

@@ -12,7 +12,6 @@
   on. This breaks ADL (argument-dependent lookup). Here we "fix" this:
  */
 namespace phys::units {
-  // using namespace phys::units::io;
   using phys::units::io::operator<<;
 } // namespace phys::units
 

Processes/Sibyll/Decay.h

@@ -15,8 +15,11 @@ namespace corsika::process {
   namespace sibyll {
 
     class Decay : public corsika::process::DecayProcess<Decay> {
+      mutable int fCount = 0;
+
     public:
       Decay() {}
+      ~Decay() { std::cout << "Sibyll::Decay n=" << fCount << std::endl; }
       void Init() {
         setHadronsUnstable();
         setTrackedParticlesStable();
@@ -142,6 +145,7 @@ namespace corsika::process {
       void DoDecay(Particle& p, Stack& s) const {
         using corsika::geometry::Point;
         using namespace corsika::units::si;
+        fCount++;
         SibStack ss;
         ss.Clear();
         // copy particle to sibyll stack

Processes/Sibyll/Interaction.h

@@ -15,9 +15,11 @@ namespace corsika::process::sibyll {
 
   class Interaction : public corsika::process::InteractionProcess<Interaction> {
 
+    mutable int fCount = 0;
+
   public:
     Interaction() {}
-    ~Interaction() {}
+    ~Interaction() { std::cout << "Sibyll::Interaction n=" << fCount << std::endl; }
 
     void Init() {
 
@@ -73,8 +75,8 @@ namespace corsika::process::sibyll {
       const hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
       const double Ecm = sqs / 1_GeV;
 
-      std::cout << "Interaction: "
-                << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
+      std::cout << "Interaction: LambdaInt: \n"
+                << " input energy: " << p.GetEnergy() / 1_GeV << endl
                 << " beam can interact:" << kBeam << endl
                 << " beam XS code:" << kBeam << endl
                 << " beam pid:" << p.GetPID() << endl
@@ -102,7 +104,6 @@ namespace corsika::process::sibyll {
                   << "nucleon mass " << nucleon_mass << std::endl;
         // calculate interaction length in medium
         GrammageType int_length = kTarget * nucleon_mass / sig;
-        // pick random step lenth
         std::cout << "Interaction: "
                   << "interaction length (g/cm2): " << int_length << std::endl;
 
@@ -110,17 +111,6 @@ namespace corsika::process::sibyll {
       }
 
       return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
-
-      /*
-        what are the units of the output? slant depth or 3space length?
-
-      */
-      //
-      // int a = 0;
-      // const double next_step = -int_length * log(s_rndm_(a));
-      // std::cout << "Interaction: "
-      //        << "next step (g/cm2): " << next_step << std::endl;
-      // return next_step;
     }
 
     template <typename Particle, typename Stack>
@@ -137,9 +127,7 @@ namespace corsika::process::sibyll {
            << "DoInteraction: " << p.GetPID() << " interaction? "
            << process::sibyll::CanInteract(p.GetPID()) << endl;
       if (process::sibyll::CanInteract(p.GetPID())) {
-        cout << "defining coordinates" << endl;
-        // coordinate system, get global frame of reference
-        CoordinateSystem& rootCS =
+        const CoordinateSystem& rootCS =
             RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
         Point pOrig = p.GetPosition();
@@ -157,14 +145,12 @@ namespace corsika::process::sibyll {
         const int kTarget = corsika::particles::Oxygen::
             GetNucleusA(); // env.GetTargetParticle().GetPID();
 
-        cout << "defining target momentum.." << endl;
         // FOR NOW: target is always at rest
-        const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
+        const EnergyType Etarget = 0_GeV + corsika::particles::Proton::GetMass();
         const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
-        cout << "target momentum (GeV/c): "
-             << pTarget.GetComponents() / 1_GeV * constants::c << endl;
-        cout << "beam momentum (GeV/c): "
-             << p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
+        cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
+        cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
+             << endl;
         cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
         cout << "time: " << tOrig << endl;
 
@@ -182,8 +168,7 @@ namespace corsika::process::sibyll {
         // total momentum
         MomentumVector Ptot = p.GetMomentum();
         // invariant mass, i.e. cm. energy
-        EnergyType Ecm =
-            sqrt(Etot * Etot - Ptot.squaredNorm()); // sqrt( 2. * E * 0.93827_GeV );
+        EnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
         /*
          get transformation between Stack-frame and SibStack-frame
          for EAS Stack-frame is lab. frame, could be different for CRMC-mode
@@ -207,8 +192,9 @@ namespace corsika::process::sibyll {
                     << " DoInteraction: should have dropped particle.." << std::endl;
           // p.Delete(); delete later... different process
         } else {
+          fCount++;
           // Sibyll does not know about units..
-          double sqs = Ecm / 1_GeV;
+          const double sqs = Ecm / 1_GeV;
           // running sibyll, filling stack
           sibyll_(kBeam, kTarget, sqs);
           // running decays
@@ -228,7 +214,8 @@ namespace corsika::process::sibyll {
           // SibStack does not know about momentum yet so we need counter to access
           // momentum array in Sibyll
           int i = -1;
-          MomentumVector Ptot_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+          MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+          EnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
           for (auto& psib : ss) {
             ++i;
             // skip particles that have decayed in Sibyll
@@ -263,10 +250,14 @@ namespace corsika::process::sibyll {
             pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
             pnew.SetPosition(pOrig);
             pnew.SetTime(tOrig);
-            Ptot_final += pnew.GetMomentum();
+            Plab_final += pnew.GetMomentum();
+            E_final += en_lab;
+            Ecm_final += psib.GetEnergy();
           }
-          // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
-          // * constants::c << endl;
+          std::cout << "conservation (all GeV): E_final=" << E_final / 1_GeV
+                    << ", Ecm_final=" << Ecm_final / 1_GeV
+                    << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
+                    << std::endl;
         }
       }
       return process::EProcessReturn::eOk;

Processes/Sibyll/ParticleConversion.h

@@ -35,7 +35,7 @@ namespace corsika::process::sibyll {
 
   SibyllCode constexpr ConvertToSibyll(corsika::particles::Code pCode) {
     return static_cast<SibyllCode>(
-      corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
+        corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
   }
 
   corsika::particles::Code constexpr ConvertFromSibyll(SibyllCode pCode) {

Processes/Sibyll/sibyll2.3c.cc

@@ -6,7 +6,7 @@
 double s_rndm_(int&) {
   static corsika::random::RNG& rng =
       corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-  
+
   std::uniform_real_distribution<double> dist;
   return dist(rng);
 }

Processes/TrackingLine/TrackingLine.h

@@ -12,6 +12,7 @@
 #define _include_corsika_processes_TrackingLine_h_
 
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/QuantityVector.h>
 #include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Vector.h>
 
@@ -75,12 +76,19 @@ namespace corsika::process {
       void Init() {}
 
       auto GetTrack(Particle const& p) {
+        using std::cout;
+        using std::endl;
         using namespace corsika::units::si;
+        using namespace corsika::geometry;
         geometry::Vector<SpeedType::dimension_type> const velocity =
             p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
 
         auto const currentPosition = p.GetPosition();
 
+        std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
+                  << std::endl;
+        std::cout << "TrackingLine   v: " << velocity.GetComponents() << std::endl;
+
         // to do: include effect of magnetic field
         geometry::Line line(currentPosition, velocity);
 
@@ -125,6 +133,8 @@ namespace corsika::process {
           min = *minIter;
         }
 
+        std::cout << " t-intersect: " << min << std::endl;
+
         return geometry::Trajectory<corsika::geometry::Line>(line, min);
       }
     };

Stack/SuperStupidStack/SuperStupidStack.h

@@ -120,7 +120,7 @@ namespace corsika::stack {
       void SetPosition(const int i, const corsika::geometry::Point& v) {
         fPosition[i] = v;
       }
-      
+
       void SetTime(const int i, const corsika::units::si::TimeType& v) { fTime[i] = v; }
 
       corsika::particles::Code GetPID(const int i) const { return fDataPID[i]; }

Stack/SuperStupidStack/testSuperStupidStack.cc

@@ -45,10 +45,9 @@ TEST_CASE("SuperStupidStack", "[stack]") {
     auto pout = s.GetNextParticle();
     REQUIRE(pout.GetPID() == particles::Code::Electron);
     REQUIRE(pout.GetEnergy() == 1.5_GeV);
-#warning Fix the next two lines:
-    // REQUIRE(pout.GetMomentum() == MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 *
-    // joule})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 * meter, 1 *
-    // meter}));
+    // REQUIRE(pout.GetMomentum() == stack::super_stupid::MomentumVector(dummyCS, {1_GeV,
+    // 1_GeV, 1_GeV})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 *
+    // meter, 1 * meter}));
     REQUIRE(pout.GetTime() == 100_s);
   }