improved ProcessSequence. Should be stable in all conditions now.

Documentation/Examples/cascade_example.cc

@@ -29,6 +29,9 @@
 
 #include <corsika/random/RNGManager.h>
 
+#include <boost/type_index.hpp>
+using boost::typeindex::type_id_with_cvr;
+
 #include <iostream>
 #include <limits>
 #include <typeinfo>
@@ -221,7 +224,9 @@ int main() {
   ProcessCut cut(8_GeV);
 
   // assemble all processes into an ordered process list
-  const auto sequence = /*p0 <<*/ sibyll << decay << cut;
+  auto sequence = p0 << sibyll << decay << cut;
+
+  //cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name() << "\n";
 
   // setup particle stack, and add primary particle
   setup::Stack stack;
@@ -230,13 +235,14 @@ int main() {
   {
     auto particle = stack.NewParticle();
     particle.SetPID(Code::Proton);
-    hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
+    hep::MomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass());
     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);
+    cout << particle.GetEnergy() / 1_GeV << endl;
   }
 
   // define air shower object, run simulation

Documentation/Examples/staticsequence_example.cc

@@ -86,7 +86,7 @@ void modular() {
   Process3 m3;
   Process4 m4(0.9);
 
-  const auto sequence = m1 << m2 << m3 << m4;
+  auto sequence = m1 << m2 << m3 << m4;
 
   DummyData p;
   DummyStack s;

Framework/Cascade/Cascade.h

@@ -63,6 +63,8 @@ namespace corsika::cascade {
       using std::endl;
       using std::log;
 
+      cout << particle.GetEnergy() / 1_GeV << endl;
+
       // determine geometric tracking
       corsika::setup::Trajectory step = fTracking.GetTrack(particle);
 

Framework/Cascade/testCascade.cc

@@ -114,7 +114,7 @@ TEST_CASE("Cascade", "[Cascade]") {
 
   stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;
-  const auto sequence = p0 << p1;
+  auto sequence = p0 << p1;
   setup::Stack stack;
 
   corsika::cascade::Cascade EAS(env, tracking, sequence, stack);

Framework/Geometry/BaseVector.h

@@ -31,10 +31,8 @@ namespace corsika::geometry {
     BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
         : qVector(pQVector)
         , cs(&pCS) {}
-        
-    auto const& GetCoordinateSystem() const {
-        return *cs;
-    }
+
+    auto const& GetCoordinateSystem() const { return *cs; }
   };
 
 } // namespace corsika::geometry

Framework/ProcessSequence/BaseProcess.h

@@ -25,12 +25,21 @@ namespace corsika::process {
 
    */
 
-  template <typename derived>
+  template <typename Derived>
   struct BaseProcess {
-    derived& GetRef() { return static_cast<derived&>(*this); }
-    const derived& GetRef() const { return static_cast<const derived&>(*this); }
+  private:
+    BaseProcess() {}
+    friend Derived;
+
+  public:
+    Derived& GetRef() { return static_cast<Derived&>(*this); }
+    const Derived& GetRef() const { return static_cast<const Derived&>(*this); }
   };
 
+  // overwrite the default trait class, to mark BaseProcess<T> as useful process
+  template <class T>
+  std::true_type is_process_impl(const BaseProcess<T>* impl);
+
 } // namespace corsika::process
 
 #endif

Framework/ProcessSequence/ContinuousProcess.h

@@ -41,6 +41,10 @@ namespace corsika::process {
     corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const;
   };
 
+  // overwrite the default trait class, to mark BaseProcess<T> as useful process
+  template <class T>
+  std::true_type is_process_impl(const ContinuousProcess<T>* impl);
+
 } // namespace corsika::process
 
 #endif

Framework/ProcessSequence/DecayProcess.h

@@ -35,17 +35,21 @@ namespace corsika::process {
     /// here starts the interface-definition part
     // -> enforce derived to implement DoDecay...
     template <typename Particle, typename Stack>
-    EProcessReturn DoDecay(Particle&, Stack&) const;
+    EProcessReturn DoDecay(Particle&, Stack&);
 
     template <typename Particle>
-    corsika::units::si::TimeType GetLifetime(Particle& p) const;
+    corsika::units::si::TimeType GetLifetime(Particle& p);
 
     template <typename Particle>
-    corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) const {
+    corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) {
       return 1. / GetRef().GetLifetime(p);
     }
   };
 
+  // overwrite the default trait class, to mark DecayProcess<T> as useful process
+  template <class T>
+  std::true_type is_process_impl(const DecayProcess<T>* impl);
+
 } // namespace corsika::process
 
 #endif

Framework/ProcessSequence/InteractionProcess.h

@@ -36,18 +36,22 @@ namespace corsika::process {
     /// here starts the interface-definition part
     // -> enforce derived to implement DoInteraction...
     template <typename P, typename S>
-    inline EProcessReturn DoInteraction(P&, S&) const;
+    inline EProcessReturn DoInteraction(P&, S&);
 
     template <typename Particle, typename Track>
-    corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t) const;
+    corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t);
 
     template <typename Particle, typename Track>
     corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p,
-                                                                        Track& t) const {
+                                                                        Track& t) {
       return 1. / GetRef().GetInteractionLength(p, t);
     }
   };
 
+  // overwrite the default trait class, to mark BaseProcess<T> as useful process
+  template <class T>
+  std::true_type is_process_impl(const InteractionProcess<T>* impl);
+
 } // namespace corsika::process
 
 #endif

Framework/ProcessSequence/ProcessSequence.h

@@ -21,6 +21,7 @@
 
 #include <cmath>
 #include <limits>
+#include <type_traits>
 
 namespace corsika::process {
 
@@ -35,20 +36,38 @@ namespace corsika::process {
      https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
    */
 
+  // define a marker (trait class) to tag any class that qualifies as "Process" for the
+  // "ProcessSequence"
+  std::false_type is_process_impl(...);
+  template <class T>
+  using is_process = decltype(is_process_impl(std::declval<T*>()));
+
   // this is a marker to track which BaseProcess is also a ProcessSequence
   template <typename T>
   struct is_process_sequence {
     static const bool value = false;
   };
 
+  /**
+     T1 and T2 are both references if possible (lvalue), otherwise
+     (rvalue) they are just classes. This allows us to handle both,
+     rvalue as well as lvalue Processes in the ProcessSequence.
+   */
   template <typename T1, typename T2>
   class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2> > {
 
+    using T1type = typename std::decay<T1>::type;
+    using T2type = typename std::decay<T2>::type;
+
   public:
-    const T1& A;
-    const T2& B;
+    T1 A; // this is a reference, if possible
+    T2 B; // this is a reference, if possible
+
+    // ProcessSequence(ProcessSequence<T1,T2>&& v) : A(v.A), B(v.B) {}
+    // ProcessSequence<T1,T2>& operator=(ProcessSequence<T1,T2>&& v) { A=v.A; B=v.B;
+    // return *this; }
 
-    ProcessSequence(const T1& in_A, const T2& in_B)
+    ProcessSequence(T1 in_A, T2 in_B)
         : A(in_A)
         , B(in_B) {}
 
@@ -56,13 +75,13 @@ namespace corsika::process {
     // void Hello() const  { detail::CallHello<T1,T2>::Call(A, B); }
 
     template <typename Particle, typename Track, typename Stack>
-    EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) const {
+    EProcessReturn DoContinuous(Particle& p, Track& t, Stack& s) {
       EProcessReturn ret = EProcessReturn::eOk;
-      if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
+      if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value ||
                     is_process_sequence<T1>::value) {
         ret |= A.DoContinuous(p, t, s);
       }
-      if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
+      if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value ||
                     is_process_sequence<T2>::value) {
         ret |= B.DoContinuous(p, t, s);
       }
@@ -70,17 +89,17 @@ namespace corsika::process {
     }
 
     template <typename Particle, typename Track>
-    corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const {
+    corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) {
       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 ||
+      if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value ||
                     is_process_sequence<T1>::value) {
         corsika::units::si::LengthType const len = A.MaxStepLength(p, track);
         max_length = std::min(max_length, len);
       }
-      if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
+      if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value ||
                     is_process_sequence<T2>::value) {
         corsika::units::si::LengthType const len = B.MaxStepLength(p, track);
         max_length = std::min(max_length, len);
@@ -89,29 +108,28 @@ namespace corsika::process {
     }
 
     template <typename Particle, typename Track>
-    corsika::units::si::GrammageType GetTotalInteractionLength(Particle& p,
-                                                               Track& t) const {
+    corsika::units::si::GrammageType GetTotalInteractionLength(Particle& p, Track& t) {
       return 1. / GetInverseInteractionLength(p, t);
     }
 
     template <typename Particle, typename Track>
-    corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
-        Particle& p, Track& t) const {
+    corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(Particle& p,
+                                                                             Track& t) {
       return GetInverseInteractionLength(p, t);
     }
 
     template <typename Particle, typename Track>
     corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p,
-                                                                        Track& t) const {
+                                                                        Track& t) {
       using namespace corsika::units::si;
 
       InverseGrammageType tot = 0 * meter * meter / gram;
 
-      if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value ||
+      if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value ||
                     is_process_sequence<T1>::value) {
         tot += A.GetInverseInteractionLength(p, t);
       }
-      if constexpr (std::is_base_of<InteractionProcess<T2>, T2>::value ||
+      if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value ||
                     is_process_sequence<T2>::value) {
         tot += B.GetInverseInteractionLength(p, t);
       }
@@ -121,15 +139,15 @@ namespace corsika::process {
     template <typename Particle, typename Stack>
     EProcessReturn SelectInteraction(
         Particle& p, Stack& s, corsika::units::si::InverseGrammageType lambda_select,
-        corsika::units::si::InverseGrammageType& lambda_inv_count) const {
+        corsika::units::si::InverseGrammageType& lambda_inv_count) {
 
-      if constexpr (is_process_sequence<T1>::value) {
+      if constexpr (is_process_sequence<T1type>::value) {
         // if A is a process sequence --> check inside
         const EProcessReturn ret =
             A.SelectInteraction(p, s, lambda_select, lambda_inv_count);
         // if A did succeed, stop routine
         if (ret != EProcessReturn::eOk) { return ret; }
-      } else if constexpr (std::is_base_of<InteractionProcess<T1>, T1>::value) {
+      } else if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value) {
         // if this is not a ContinuousProcess --> evaluate probability
         lambda_inv_count += A.GetInverseInteractionLength(p, s);
         // check if we should execute THIS process and then EXIT
@@ -145,7 +163,7 @@ namespace corsika::process {
             B.SelectInteraction(p, s, lambda_select, lambda_inv_count);
         // if A did succeed, stop routine
         if (ret != EProcessReturn::eOk) { return ret; }
-      } else if constexpr (std::is_base_of<InteractionProcess<T2>, T2>::value) {
+      } else if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value) {
         // if this is not a ContinuousProcess --> evaluate probability
         lambda_inv_count += B.GetInverseInteractionLength(p, s);
         // check if we should execute THIS process and then EXIT
@@ -158,26 +176,26 @@ namespace corsika::process {
     }
 
     template <typename Particle>
-    corsika::units::si::TimeType GetTotalLifetime(Particle& p) const {
+    corsika::units::si::TimeType GetTotalLifetime(Particle& p) {
       return 1. / GetInverseLifetime(p);
     }
 
     template <typename Particle>
-    corsika::units::si::InverseTimeType GetTotalInverseLifetime(Particle& p) const {
+    corsika::units::si::InverseTimeType GetTotalInverseLifetime(Particle& p) {
       return GetInverseLifetime(p);
     }
 
     template <typename Particle>
-    corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) const {
+    corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) {
       using namespace corsika::units::si;
 
       corsika::units::si::InverseTimeType tot = 0 / second;
 
-      if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value ||
+      if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value ||
                     is_process_sequence<T1>::value) {
         tot += A.GetInverseLifetime(p);
       }
-      if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value ||
+      if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value ||
                     is_process_sequence<T2>::value) {
         tot += B.GetInverseLifetime(p);
       }
@@ -186,15 +204,15 @@ namespace corsika::process {
 
     // select decay process
     template <typename Particle, typename Stack>
-    EProcessReturn SelectDecay(
-        Particle& p, Stack& s, corsika::units::si::InverseTimeType decay_select,
-        corsika::units::si::InverseTimeType& decay_inv_count) const {
+    EProcessReturn SelectDecay(Particle& p, Stack& s,
+                               corsika::units::si::InverseTimeType decay_select,
+                               corsika::units::si::InverseTimeType& decay_inv_count) {
       if constexpr (is_process_sequence<T1>::value) {
         // if A is a process sequence --> check inside
         const EProcessReturn ret = A.SelectDecay(p, s, decay_select, decay_inv_count);
         // if A did succeed, stop routine
         if (ret != EProcessReturn::eOk) { return ret; }
-      } else if constexpr (std::is_base_of<DecayProcess<T1>, T1>::value) {
+      } else if constexpr (std::is_base_of<DecayProcess<T1type>, T1type>::value) {
         // if this is not a ContinuousProcess --> evaluate probability
         decay_inv_count += A.GetInverseLifetime(p);
         // check if we should execute THIS process and then EXIT
@@ -210,7 +228,7 @@ namespace corsika::process {
         const EProcessReturn ret = B.SelectDecay(p, s, decay_select, decay_inv_count);
         // if A did succeed, stop routine
         if (ret != EProcessReturn::eOk) { return ret; }
-      } else if constexpr (std::is_base_of<DecayProcess<T2>, T2>::value) {
+      } else if constexpr (std::is_base_of<DecayProcess<T2type>, T2type>::value) {
         // if this is not a ContinuousProcess --> evaluate probability
         decay_inv_count += B.GetInverseLifetime(p);
         // check if we should execute THIS process and then EXIT
@@ -222,10 +240,9 @@ namespace corsika::process {
       return EProcessReturn::eOk;
     }
 
-    /// TODO the const_cast is not nice, think about the constness here
-    void Init() const {
-      const_cast<T1*>(&A)->Init();
-      const_cast<T2*>(&B)->Init();
+    void Init() {
+      A.Init();
+      B.Init();
     }
   };
 
@@ -234,28 +251,46 @@ namespace corsika::process {
   /// must be allowed, this is why we define a macro to define all
   /// combinations here:
 
-#define OPSEQ(C1, C2)                                                                 \
-  template <typename T1, typename T2>                                                 \
-  inline const ProcessSequence<T1, T2> operator<<(const C1<T1>& A, const C2<T2>& B) { \
-    return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef());                           \
+  template <
+      typename P1, typename P2,
+      typename std::enable_if<is_process<typename std::decay<P1>::type>::value &&
+                              is_process<typename std::decay<P2>::type>::value>::type...>
+  inline auto operator<<(P1&& A, P2&& B) -> ProcessSequence<P1, P2> {
+    return ProcessSequence<P1, P2>(A.GetRef(), B.GetRef());
   }
 
-  OPSEQ(BaseProcess, BaseProcess)
-  OPSEQ(BaseProcess, InteractionProcess)
-  OPSEQ(BaseProcess, ContinuousProcess)
-  OPSEQ(BaseProcess, DecayProcess)
-  OPSEQ(ContinuousProcess, BaseProcess)
-  OPSEQ(ContinuousProcess, InteractionProcess)
-  OPSEQ(ContinuousProcess, ContinuousProcess)
-  OPSEQ(ContinuousProcess, DecayProcess)
-  OPSEQ(InteractionProcess, BaseProcess)
-  OPSEQ(InteractionProcess, InteractionProcess)
-  OPSEQ(InteractionProcess, ContinuousProcess)
-  OPSEQ(InteractionProcess, DecayProcess)
-  OPSEQ(DecayProcess, BaseProcess)
-  OPSEQ(DecayProcess, InteractionProcess)
-  OPSEQ(DecayProcess, ContinuousProcess)
-  OPSEQ(DecayProcess, DecayProcess)
+  /* #define OPSEQ(C1, C2) \ */
+  /*   template < \ */
+  /*       typename P1, typename P2, \ */
+  /*       typename std::enable_if<is_process<typename std::decay<P1>::type>::value &&
+   * \ */
+  /*                               is_process<typename
+   * std::decay<P2>::type>::value>::type...> \ */
+  /*   inline auto operator+(P1&& A, P2&& B)->ProcessSequence<P1, P2> { \ */
+  /*     return ProcessSequence<P1, P2>(A.GetRef(), B.GetRef()); \ */
+  /*   } */
+
+  /*   /\*template <typename T1, typename T2>				\ */
+  /* inline ProcessSequence<T1, T2> operator%(C1<T1>& A, C2<T2>& B) {	\ */
+  /* return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef());	\ */
+  /* }*\/ */
+
+  /*   OPSEQ(BaseProcess, BaseProcess) */
+  /*   OPSEQ(BaseProcess, InteractionProcess) */
+  /*   OPSEQ(BaseProcess, ContinuousProcess) */
+  /*   OPSEQ(BaseProcess, DecayProcess) */
+  /*   OPSEQ(ContinuousProcess, BaseProcess) */
+  /*   OPSEQ(ContinuousProcess, InteractionProcess) */
+  /*   OPSEQ(ContinuousProcess, ContinuousProcess) */
+  /*   OPSEQ(ContinuousProcess, DecayProcess) */
+  /*   OPSEQ(InteractionProcess, BaseProcess) */
+  /*   OPSEQ(InteractionProcess, InteractionProcess) */
+  /*   OPSEQ(InteractionProcess, ContinuousProcess) */
+  /*   OPSEQ(InteractionProcess, DecayProcess) */
+  /*   OPSEQ(DecayProcess, BaseProcess) */
+  /*   OPSEQ(DecayProcess, InteractionProcess) */
+  /*   OPSEQ(DecayProcess, ContinuousProcess) */
+  /*   OPSEQ(DecayProcess, DecayProcess) */
 
   /// marker to identify objectas ProcessSequence
   template <typename A, typename B>

Framework/ProcessSequence/testProcessSequence.cc

@@ -188,7 +188,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
     Process3 m3(2);
     Process4 m4(3);
 
-    const auto sequence = m1 << m2 << m3 << m4;
+    auto sequence = m1 << m2 << m3 << m4;
 
     globalCount = 0;
     sequence.Init();
@@ -208,7 +208,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
     DummyStack s;
     DummyTrajectory t;
 
-    const auto sequence2 = cp1 << m2 << m3;
+    auto sequence2 = cp1 << m2 << m3;
     GrammageType const tot = sequence2.GetTotalInteractionLength(s, t);
     InverseGrammageType const tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
     cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
@@ -222,7 +222,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
 
     DummyStack s;
 
-    const auto sequence2 = cp1 << m2 << m3 << d3;
+    auto sequence2 = cp1 << m2 << m3 << d3;
     TimeType const tot = sequence2.GetTotalLifetime(s);
     InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(s);
     cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
@@ -235,7 +235,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
     Process2 m2(1);
     Process3 m3(2);
 
-    const auto sequence2 = cp1 << m2 << m3 << cp2;
+    auto sequence2 = cp1 << m2 << m3 << cp2;
 
     DummyData p;
     DummyStack s;

Framework/Utilities/testCOMBoost.cc

@@ -55,7 +55,7 @@ TEST_CASE("boosts") {
   // boost projecticle
   auto const [eProjectileCoM, pProjectileCoM] =
       boost.toCoM(eProjectileLab, pProjectileLab);
-      
+
   // boost target
   auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
 

Processes/Sibyll/Decay.h

@@ -25,11 +25,11 @@ namespace corsika::process {
         setTrackedParticlesStable();
       }
 
-      void setTrackedParticlesStable() const {
+      void setTrackedParticlesStable() {
         /*
-          Sibyll is hadronic generator
-          only hadrons decay
-        */
+           Sibyll is hadronic generator
+           only hadrons decay
+         */
         // set particles unstable
         setHadronsUnstable();
         // make tracked particles stable
@@ -45,17 +45,17 @@ namespace corsika::process {
         }
       }
 
-      void setUnstable(const corsika::particles::Code pCode) const {
+      void setUnstable(const corsika::particles::Code pCode) {
         int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
         s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
       }
 
-      void setStable(const corsika::particles::Code pCode) const {
+      void setStable(const corsika::particles::Code pCode) {
         int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
         s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
       }
 
-      void setAllStable() const {
+      void setAllStable() {
         // name? also makes EM particles stable
 
         using std::cout;
@@ -75,7 +75,7 @@ namespace corsika::process {
         }
       }
 
-      void setHadronsUnstable() const {
+      void setHadronsUnstable() {
 
         using std::cout;
         using std::endl;
@@ -118,7 +118,7 @@ namespace corsika::process {
       }
 
       template <typename Particle>
-      corsika::units::si::TimeType GetLifetime(Particle& p) const {
+      corsika::units::si::TimeType GetLifetime(Particle& p) {
         using std::cout;
         using std::endl;
         using namespace corsika::units::si;
@@ -142,7 +142,7 @@ namespace corsika::process {
       }
 
       template <typename Particle, typename Stack>
-      void DoDecay(Particle& p, Stack& s) const {
+      void DoDecay(Particle& p, Stack& s) {
         using corsika::geometry::Point;
         using namespace corsika::units::si;
         fCount++;

Processes/Sibyll/Interaction.h

@@ -35,7 +35,7 @@ namespace corsika::process::sibyll {
     }
 
     template <typename Particle, typename Track>
-    corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) const {
+    corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) {
 
       using namespace corsika::units;
       using namespace corsika::units::hep;
@@ -114,7 +114,7 @@ namespace corsika::process::sibyll {
     }
 
     template <typename Particle, typename Stack>
-    corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const {
+    corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
 
       using namespace corsika::units;
       using namespace corsika::units::hep;

Processes/StackInspector/StackInspector.cc

@@ -37,7 +37,7 @@ StackInspector<Stack>::~StackInspector() {}
 
 template <typename Stack>
 process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
-                                                            Stack& s) const {
+                                                            Stack& s) {
 
   if (!fReport) return process::EProcessReturn::eOk;
   [[maybe_unused]] int i = 0;
@@ -60,8 +60,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
 }
 
 template <typename Stack>
-corsika::units::si::LengthType StackInspector<Stack>::MaxStepLength(
-    Particle&, setup::Trajectory&) const {
+corsika::units::si::LengthType StackInspector<Stack>::MaxStepLength(Particle&,
+                                                                    setup::Trajectory&) {
   return std::numeric_limits<double>::infinity() * meter;
 }
 

Processes/StackInspector/StackInspector.h

@@ -31,14 +31,14 @@ namespace corsika::process {
       ~StackInspector();
 
       void Init();
-      EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
+      EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s);
       //~ template <typename Particle>
       corsika::units::si::LengthType MaxStepLength(Particle&,
-                                                   corsika::setup::Trajectory&) const;
+                                                   corsika::setup::Trajectory&);
 
     private:
       bool fReport;
-      mutable int fCountStep = 0;
+      int fCountStep = 0;
     };
 
   } // namespace stack_inspector

Processes/TrackingLine/TrackingLine.h

@@ -87,6 +87,9 @@ namespace corsika::process {
 
         std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
                   << std::endl;
+        std::cout << "TrackingLine   E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+        std::cout << "TrackingLine   p: " << p.GetMomentum().GetComponents() / 1_GeV
+                  << " GeV " << std::endl;
         std::cout << "TrackingLine   v: " << velocity.GetComponents() << std::endl;
 
         // to do: include effect of magnetic field