Added SecondariesProcess, some wording -> style

Framework/ProcessSequence/CMakeLists.txt

@@ -12,6 +12,7 @@ set (
   CORSIKAprocesssequence_HEADERS
   BaseProcess.h
   ContinuousProcess.h
+  SecondariesProcess.h
   InteractionProcess.h
   DecayProcess.h
   ProcessSequence.h

Framework/ProcessSequence/ContinuousProcess.h

@@ -33,8 +33,8 @@ namespace corsika::process {
 
     // here starts the interface part
     // -> enforce derived to implement DoContinuous...
-    template <typename Particle, typename Track, typename Stack>
-    EProcessReturn DoContinuous(Particle&, Track&, Stack&) const;
+    template <typename Particle, typename Track>
+    EProcessReturn DoContinuous(Particle&, Track&) const;
 
     // -> enforce derived to implement MaxStepLength...
     template <typename Particle, typename Track>

Framework/ProcessSequence/DecayProcess.h

@@ -34,15 +34,15 @@ namespace corsika::process {
 
     /// here starts the interface-definition part
     // -> enforce derived to implement DoDecay...
-    template <typename Particle, typename Stack>
-    EProcessReturn DoDecay(Particle&, Stack&);
+    template <typename Particle>
+    EProcessReturn DoDecay(Particle&);
 
     template <typename Particle>
     corsika::units::si::TimeType GetLifetime(Particle& p);
 
     template <typename Particle>
-    corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) {
-      return 1. / GetRef().GetLifetime(p);
+    corsika::units::si::InverseTimeType GetInverseLifetime(Particle& vP) {
+      return 1. / GetRef().GetLifetime(vP);
     }
   };
 

Framework/ProcessSequence/InteractionProcess.h

@@ -35,8 +35,8 @@ namespace corsika::process {
 
     /// here starts the interface-definition part
     // -> enforce derived to implement DoInteraction...
-    template <typename P, typename S>
-    inline EProcessReturn DoInteraction(P&, S&);
+    template <typename Particle>
+    inline EProcessReturn DoInteraction(Particle&);
 
     template <typename Particle, typename Track>
     corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track& t);

Framework/ProcessSequence/ProcessSequence.h

@@ -17,6 +17,7 @@
 #include <corsika/process/DecayProcess.h>
 #include <corsika/process/InteractionProcess.h>
 #include <corsika/process/ProcessReturn.h>
+#include <corsika/process/SecondariesProcess.h>
 #include <corsika/units/PhysicalUnits.h>
 
 #include <cmath>
@@ -70,78 +71,93 @@ namespace corsika::process {
     // example for a trait-based call:
     // 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) {
+    template <typename TParticle, typename TTrack>
+    EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) {
       EProcessReturn ret = EProcessReturn::eOk;
       if constexpr (std::is_base_of<ContinuousProcess<T1type>, T1type>::value ||
                     is_process_sequence<T1>::value) {
-        ret |= A.DoContinuous(p, t, s);
+        ret |= A.DoContinuous(vP, vT);
       }
       if constexpr (std::is_base_of<ContinuousProcess<T2type>, T2type>::value ||
                     is_process_sequence<T2>::value) {
-        ret |= B.DoContinuous(p, t, s);
+        ret |= B.DoContinuous(vP, vT);
       }
       return ret;
     }
 
-    template <typename Particle, typename Track>
-    corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) {
+    template <typename TSecondaries>
+    EProcessReturn DoSecondaries(TSecondaries& vS) {
+      EProcessReturn ret = EProcessReturn::eOk;
+      if constexpr (std::is_base_of<SecondariesProcess<T1type>, T1type>::value ||
+                    is_process_sequence<T1>::value) {
+        ret |= A.DoSecondaries(vS);
+      }
+      if constexpr (std::is_base_of<SecondariesProcess<T2type>, T2type>::value ||
+                    is_process_sequence<T2>::value) {
+        ret |= B.DoSecondaries(vS);
+      }
+      return ret;
+    }
+
+    template <typename TParticle, typename TTrack>
+    corsika::units::si::LengthType MaxStepLength(TParticle& vP, TTrack& vTrack) {
       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<T1type>, T1type>::value ||
                     is_process_sequence<T1>::value) {
-        corsika::units::si::LengthType const len = A.MaxStepLength(p, track);
+        corsika::units::si::LengthType const len = A.MaxStepLength(vP, vTrack);
         max_length = std::min(max_length, len);
       }
       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);
+        corsika::units::si::LengthType const len = B.MaxStepLength(vP, vTrack);
         max_length = std::min(max_length, len);
       }
       return max_length;
     }
 
-    template <typename Particle, typename Track>
-    corsika::units::si::GrammageType GetTotalInteractionLength(Particle& p, Track& t) {
-      return 1. / GetInverseInteractionLength(p, t);
+    template <typename TParticle, typename TTrack>
+    corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP,
+                                                               TTrack& vT) {
+      return 1. / GetInverseInteractionLength(vP, vT);
     }
 
-    template <typename Particle, typename Track>
-    corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(Particle& p,
-                                                                             Track& t) {
-      return GetInverseInteractionLength(p, t);
+    template <typename TParticle, typename TTrack>
+    corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
+        TParticle& vP, TTrack& vT) {
+      return GetInverseInteractionLength(vP, vT);
     }
 
-    template <typename Particle, typename Track>
-    corsika::units::si::InverseGrammageType GetInverseInteractionLength(Particle& p,
-                                                                        Track& t) {
+    template <typename TParticle, typename TTrack>
+    corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP,
+                                                                        TTrack& vT) {
       using namespace corsika::units::si;
 
       InverseGrammageType tot = 0 * meter * meter / gram;
 
       if constexpr (std::is_base_of<InteractionProcess<T1type>, T1type>::value ||
                     is_process_sequence<T1>::value) {
-        tot += A.GetInverseInteractionLength(p, t);
+        tot += A.GetInverseInteractionLength(vP, vT);
       }
       if constexpr (std::is_base_of<InteractionProcess<T2type>, T2type>::value ||
                     is_process_sequence<T2>::value) {
-        tot += B.GetInverseInteractionLength(p, t);
+        tot += B.GetInverseInteractionLength(vP, vT);
       }
       return tot;
     }
 
-    template <typename Particle, typename Track, typename Stack>
+    template <typename TParticle, typename TSecondaries, typename TTrack>
     EProcessReturn SelectInteraction(
-        Particle& vP, Track& vT, Stack& vS,
+        TParticle& vP, TSecondaries& vS, TTrack& vT,
         [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
         corsika::units::si::InverseGrammageType& lambda_inv_count) {
 
       if constexpr (is_process_sequence<T1type>::value) {
         // if A is a process sequence --> check inside
         const EProcessReturn ret =
-            A.SelectInteraction(vP, vT, vS, lambda_select, lambda_inv_count);
+            A.SelectInteraction(vP, vS, vT, 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<T1type>, T1type>::value) {
@@ -149,7 +165,7 @@ namespace corsika::process {
         lambda_inv_count += A.GetInverseInteractionLength(vP, vT);
         // check if we should execute THIS process and then EXIT
         if (lambda_select < lambda_inv_count) {
-          A.DoInteraction(vP, vS);
+          A.DoInteraction(vS);
           return EProcessReturn::eInteracted;
         }
       } // end branch A
@@ -157,7 +173,7 @@ namespace corsika::process {
       if constexpr (is_process_sequence<T2>::value) {
         // if A is a process sequence --> check inside
         const EProcessReturn ret =
-            B.SelectInteraction(vP, vT, vS, lambda_select, lambda_inv_count);
+            B.SelectInteraction(vP, vS, vT, 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<T2type>, T2type>::value) {
@@ -165,25 +181,25 @@ namespace corsika::process {
         lambda_inv_count += B.GetInverseInteractionLength(vP, vT);
         // check if we should execute THIS process and then EXIT
         if (lambda_select < lambda_inv_count) {
-          B.DoInteraction(vP, vS);
+          B.DoInteraction(vS);
           return EProcessReturn::eInteracted;
         }
       } // end branch A
       return EProcessReturn::eOk;
     }
 
-    template <typename Particle>
-    corsika::units::si::TimeType GetTotalLifetime(Particle& p) {
+    template <typename TParticle>
+    corsika::units::si::TimeType GetTotalLifetime(TParticle& p) {
       return 1. / GetInverseLifetime(p);
     }
 
-    template <typename Particle>
-    corsika::units::si::InverseTimeType GetTotalInverseLifetime(Particle& p) {
+    template <typename TParticle>
+    corsika::units::si::InverseTimeType GetTotalInverseLifetime(TParticle& p) {
       return GetInverseLifetime(p);
     }
 
-    template <typename Particle>
-    corsika::units::si::InverseTimeType GetInverseLifetime(Particle& p) {
+    template <typename TParticle>
+    corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& p) {
       using namespace corsika::units::si;
 
       corsika::units::si::InverseTimeType tot = 0 / second;
@@ -200,38 +216,38 @@ namespace corsika::process {
     }
 
     // select decay process
-    template <typename Particle, typename Stack>
+    template <typename TParticle, typename TSecondaries>
     EProcessReturn SelectDecay(
-        Particle& p, Stack& s,
+        TParticle& vP, TSecondaries& vS,
         [[maybe_unused]] 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);
+        const EProcessReturn ret = A.SelectDecay(vP, vS, 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<T1type>, T1type>::value) {
         // if this is not a ContinuousProcess --> evaluate probability
-        decay_inv_count += A.GetInverseLifetime(p);
+        decay_inv_count += A.GetInverseLifetime(vP);
         // check if we should execute THIS process and then EXIT
         if (decay_select < decay_inv_count) { // more pedagogical: rndm_select <
                                               // decay_inv_count / decay_inv_tot
-          A.DoDecay(p, s);
+          A.DoDecay(vS);
           return EProcessReturn::eDecayed;
         }
       } // end branch A
 
       if constexpr (is_process_sequence<T2>::value) {
         // if A is a process sequence --> check inside
-        const EProcessReturn ret = B.SelectDecay(p, s, decay_select, decay_inv_count);
+        const EProcessReturn ret = B.SelectDecay(vP, vS, 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<T2type>, T2type>::value) {
         // if this is not a ContinuousProcess --> evaluate probability
-        decay_inv_count += B.GetInverseLifetime(p);
+        decay_inv_count += B.GetInverseLifetime(vP);
         // check if we should execute THIS process and then EXIT
         if (decay_select < decay_inv_count) {
-          B.DoDecay(p, s);
+          B.DoDecay(vS);
           return EProcessReturn::eDecayed;
         }
       } // end branch B
@@ -253,8 +269,8 @@ namespace corsika::process {
       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());
+  inline auto operator<<(P1&& vA, P2&& vB) -> ProcessSequence<P1, P2> {
+    return ProcessSequence<P1, P2>(vA.GetRef(), vB.GetRef());
   }
 
   /// marker to identify objectas ProcessSequence

Framework/ProcessSequence/SecondaryProcess.h

+
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _include_corsika_secondariesprocess_h_
+#define _include_corsika_secondariesprocess_h_
+
+#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::process {
+
+  /**
+     \class SecondariesProcess
+
+     The structural base type of a process object in a
+     ProcessSequence. Both, the ProcessSequence and all its elements
+     are of type SecondariesProcess<T>
+
+   */
+
+  template <typename derived>
+  struct SecondariesProcess {
+
+    derived& GetRef() { return static_cast<derived&>(*this); }
+    const derived& GetRef() const { return static_cast<const derived&>(*this); }
+
+    /// here starts the interface-definition part
+    // -> enforce derived to implement DoSecondaries...
+    template <typename Particle>
+    inline EProcessReturn DoSecondaries(Particle&);
+  };
+
+  // overwrite the default trait class, to mark BaseProcess<T> as useful process
+  template <class T>
+  std::true_type is_process_impl(const SecondariesProcess<T>* impl);
+
+} // namespace corsika::process
+
+#endif

Framework/ProcessSequence/testProcessSequence.cc

@@ -39,8 +39,8 @@ public:
     assert(globalCount == fV);
     globalCount++;
   }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+  template <typename D, typename T>
+  inline EProcessReturn DoContinuous(D& d, T&) const {
     cout << "ContinuousProcess1::DoContinuous" << endl;
     for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
@@ -58,8 +58,8 @@ public:
     assert(globalCount == fV);
     globalCount++;
   }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+  template <typename D, typename T>
+  inline EProcessReturn DoContinuous(D& d, T&) const {
     cout << "ContinuousProcess2::DoContinuous" << endl;
     for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
     return EProcessReturn::eOk;
@@ -95,8 +95,8 @@ public:
     assert(globalCount == fV);
     globalCount++;
   }
-  template <typename Particle, typename Stack>
-  inline EProcessReturn DoInteraction(Particle&, Stack&) const {
+  template <typename Particle>
+  inline EProcessReturn DoInteraction(Particle&) const {
     cout << "Process2::DoInteraction" << endl;
     return EProcessReturn::eOk;
   }
@@ -118,8 +118,8 @@ public:
     assert(globalCount == fV);
     globalCount++;
   }
-  template <typename Particle, typename Stack>
-  inline EProcessReturn DoInteraction(Particle&, Stack&) const {
+  template <typename Particle>
+  inline EProcessReturn DoInteraction(Particle&) const {
     cout << "Process3::DoInteraction" << endl;
     return EProcessReturn::eOk;
   }
@@ -141,14 +141,14 @@ public:
     assert(globalCount == fV);
     globalCount++;
   }
-  template <typename D, typename T, typename S>
-  inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+  template <typename D, typename T>
+  inline EProcessReturn DoContinuous(D& d, T&) const {
     for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; }
     return EProcessReturn::eOk;
   }
   // inline double MinStepLength(D& d) {
-  template <typename Particle, typename Stack>
-  EProcessReturn DoInteraction(Particle&, Stack&) const {
+  template <typename Particle>
+  EProcessReturn DoInteraction(Particle&) const {
     return EProcessReturn::eOk;
   }
 };
@@ -168,8 +168,8 @@ public:
   TimeType GetLifetime(Particle&) const {
     return 1_s;
   }
-  template <typename Particle, typename Stack>
-  EProcessReturn DoDecay(Particle&, Stack&) const {
+  template <typename Particle>
+  EProcessReturn DoDecay(Particle&) const {
     return EProcessReturn::eOk;
   }
 };
@@ -177,7 +177,6 @@ public:
 struct DummyData {
   double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 };
-struct DummyStack {};
 struct DummyTrajectory {};
 
 TEST_CASE("Process Sequence", "[Process Sequence]") {
@@ -205,12 +204,13 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
     Process2 m2(1);
     Process3 m3(2);
 
-    DummyStack s;
-    DummyTrajectory t;
+    DummyData particle;
+    DummyTrajectory track;
 
     auto sequence2 = cp1 << m2 << m3;
-    GrammageType const tot = sequence2.GetTotalInteractionLength(s, t);
-    InverseGrammageType const tot_inv = sequence2.GetTotalInverseInteractionLength(s, t);
+    GrammageType const tot = sequence2.GetTotalInteractionLength(particle, track);
+    InverseGrammageType const tot_inv =
+        sequence2.GetTotalInverseInteractionLength(particle, track);
     cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
   }
 
@@ -220,11 +220,11 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
     Process3 m3(2);
     Decay1 d3(2);
 
-    DummyStack s;
+    DummyData particle;
 
     auto sequence2 = cp1 << m2 << m3 << d3;
-    TimeType const tot = sequence2.GetTotalLifetime(s);
-    InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(s);
+    TimeType const tot = sequence2.GetTotalLifetime(particle);
+    InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(particle);
     cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
   }
 
@@ -237,27 +237,24 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
 
     auto sequence2 = cp1 << m2 << m3 << cp2;
 
-    DummyData p;
-    DummyStack s;
-    DummyTrajectory t;
+    DummyData particle;
+    DummyTrajectory track;
 
     cout << "-->init sequence2" << endl;
     globalCount = 0;
     sequence2.Init();
     cout << "-->docont" << endl;
 
-    sequence2.DoContinuous(p, t, s);
+    sequence2.DoContinuous(particle, track);
     cout << "-->dodisc" << endl;
-    // sequence2.DoInteraction(p, s);
     cout << "-->done" << endl;
 
     const int nLoop = 5;
     cout << "Running loop with n=" << nLoop << endl;
-    for (int i = 0; i < nLoop; ++i) {
-      sequence2.DoContinuous(p, t, s);
-      // sequence2.DoInteraction(p, s);
+    for (int i = 0; i < nLoop; ++i) { sequence2.DoContinuous(particle, track); }
+    for (int i = 0; i < nData; i++) {
+      cout << "data[" << i << "]=" << particle.p[i] << endl;
     }
-    for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; }
     cout << "done" << endl;
   }
 }