first improve log output

corsika/detail/modules/ObservationPlane.inl

@@ -7,6 +7,8 @@
  */
 
 #include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
 
 #include <fstream>
 
@@ -27,115 +29,64 @@ namespace corsika {
 
   ProcessReturn ObservationPlane::doContinuous(
       corsika::setup::Stack::particle_type& particle,
-      corsika::setup::Trajectory& trajectory) {
-    TimeType const timeOfIntersection =
-        (plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
-        trajectory.getVelocity(0).dot(plane_.getNormal());
+      corsika::setup::Trajectory& trajectory,
+      bool const stepLimit) {
 
-    if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
+      auto const* volumeNode = particle.getNode();
 
-    if (plane_.isAbove(trajectory.getPosition(0)) ==
-        plane_.isAbove(trajectory.getPosition(1))) {
-      return ProcessReturn::Ok;
-    }
-
-    const auto energy = particle.getEnergy();
-    auto const displacement = trajectory.getPosition(1) - plane_.getCenter();
-
-    outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
-                  << ' ' << displacement.dot(xAxis_) / 1_m << ' '
-                  << displacement.dot(yAxis_) / 1_m
-                  << (trajectory.getPosition(1) - plane_.getCenter()).getNorm() / 1_m
-                  << std::endl;
+      if (!stepLimit) {
+	return ProcessReturn::Ok;
+      }
 
-    if (deleteOnHit_) {
-      count_ground_++;
-      energy_ground_ += energy;
-      particle.erase();
-      return ProcessReturn::ParticleAbsorbed;
-    } else {
-      return ProcessReturn::Ok;
+      HEPEnergyType const energy = particle.getEnergy();
+      TimeType const timeOfIntersection = particle.getTime();
+      Point const pointOfIntersection = particle.getPosition();
+      Vector const displacement = pointOfIntersection - plane_.getCenter();
+
+      outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
+		    << ' ' << displacement.dot(xAxis_) / 1_m << ' '
+		    << displacement.dot(yAxis_) / 1_m
+		    << (pointOfIntersection - plane_.getCenter()).getNorm() / 1_m
+		    << '\n';
+      
+      if (deleteOnHit_) {
+	count_ground_++;
+	energy_ground_ += energy;
+	particle.erase();
+	return ProcessReturn::ParticleAbsorbed;
+      } else {
+	return ProcessReturn::Ok;
     }
   }
 
   LengthType ObservationPlane::getMaxStepLength(
-      corsika::setup::Stack::particle_type const& vParticle,
+      corsika::setup::Stack::particle_type const& particle,
       corsika::setup::Trajectory const& trajectory) {
 
-    int chargeNumber;
-    if (is_nucleus(vParticle.getPID())) {
-      chargeNumber = vParticle.getNuclearZ();
-    } else {
-      chargeNumber = get_charge_number(vParticle.getPID());
-    }
-    auto const* currentLogicalVolumeNode = vParticle.getNode();
-    auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
-        vParticle.getPosition());
-    auto direction = trajectory.getVelocity(0).normalized();
-
-    if (chargeNumber != 0 &&
-        abs(plane_.getNormal().dot(
-            trajectory.getLine().getVelocity().cross(magneticfield))) *
-                1_s / 1_m / 1_T >
-            1e-6) {
-      auto const* currentLogicalVolumeNode = vParticle.getNode();
-      auto magneticfield =
-          currentLogicalVolumeNode->getModelProperties().getMagneticField(
-              vParticle.getPosition());
-      auto k =
-          chargeNumber * constants::c * 1_eV / (vParticle.getMomentum().getNorm() * 1_V);
-
-      if (direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
-              (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
-               plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
-          0) {
-        return std::numeric_limits<double>::infinity() * 1_m;
-      }
+    auto const& volumeNode = particle.getNode();
 
-      LengthType MaxStepLength1 =
-          (sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
-                (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
-                 plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
-           direction.dot(plane_.getNormal()) / direction.getNorm()) /
-          (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
-      LengthType MaxStepLength2 =
-          (-sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
-                 (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
-                  plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
-           direction.dot(plane_.getNormal()) / direction.getNorm()) /
-          (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
-      if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
-        return std::numeric_limits<double>::infinity() * 1_m;
-      } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
-        std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
-        return MaxStepLength2 *
-               (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
-                   .getNorm() *
-               1.001;
-      } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
-        std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
-        return MaxStepLength1 *
-               (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
-                   .getNorm() *
-               1.001;
-      }
-    }
+    typedef typename std::remove_const_t<
+      std::remove_reference_t<decltype(volumeNode->getModelProperties())>>
+      medium_type;
+    
+    Intersections const intersection = setup::Tracking::intersect<corsika::setup::Stack::particle_type, medium_type>(particle, plane_,
+															    volumeNode->getModelProperties());
 
-    TimeType const timeOfIntersection =
-        (plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
-        trajectory.getVelocity(0).dot(plane_.getNormal());
 
-    if (timeOfIntersection < TimeType::zero()) {
+    TimeType const timeOfIntersection = intersection.getEntry();
+    
+    if (timeOfIntersection < TimeType::zero()) { 
+      return std::numeric_limits<double>::infinity() * 1_m;
+    }    
+    if (timeOfIntersection > trajectory.getDuration()) {
       return std::numeric_limits<double>::infinity() * 1_m;
     }
 
     double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
 
     auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
-    auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm() * 1.0001;
-    CORSIKA_LOG_TRACE("ObservationPlane w/o magnetic field: getMaxStepLength l={} m",
+    auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm() * 1.001; // make max. step length a bit longer, to assure crossing
+    CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m",
                       dist / 1_m);
     return dist;
   }

corsika/framework/core/Logging.hpp

@@ -35,9 +35,9 @@
 // if this is a Debug build, include debug messages in objects
 #ifdef DEBUG
 // trace is the highest level of logging (ALL messages will be printed)
-#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_DEBUG
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_TRACE
 #else // otherwise, remove everything but "error" and worse messages
-#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_INFO
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_DEBUG
 #endif
 
 #include <spdlog/fmt/ostr.h> // will output whenerver a streaming operator is found

corsika/framework/process/ContinuousProcess.hpp

@@ -30,7 +30,7 @@ namespace corsika {
     // here starts the interface part
     // -> enforce TDerived to implement DoContinuous...
     template <typename TParticle, typename TTrack>
-    ProcessReturn doContinuous(TParticle&, TTrack const&) const;
+    ProcessReturn doContinuous(TParticle&, TTrack const&, bool const stepLimit) const;
 
     // -> enforce TDerived to implement MaxStepLength...
     template <typename TParticle, typename TTrack>

corsika/framework/process/ProcessSequence.hpp

@@ -26,6 +26,21 @@
 
 namespace corsika {
 
+  namespace detail {
+    template <typename TProcess, int N=0, typename Enable=void>
+    struct CountContinuous {
+      enum { count = N };
+    }; 
+
+    template <typename TProcess, int N=0>              
+    struct CountContinuous<TProcess, N, 
+                           typename std::enable_if_t<std::is_base_of_v<ContinuousProcess<typename std::decay_t<TProcess>>>> {
+      enum { count = N+1 };
+    };
+  }
+
+    
+
   /**
    *
    *  Definition of a static process list/sequence
@@ -39,13 +54,12 @@ namespace corsika {
    *  they are just classes. This allows us to handle both, rvalue as
    *  well as lvalue Processes in the ProcessSequence.
    *
-   *  The sequence, and the processes use CRTP.
+   *  (The sequence, and the processes use the CRTP, curiously recurring template pattern).
    *
-   *  \todo There are several FIXME's in the ProcessSequence.inl due to
-   *  outstanding migration of SecondaryView::parent()
    **/
 
-  template <typename TProcess1, typename TProcess2 = NullModel>
+  template <typename TProcess1, typename TProcess2 = NullModel, 
+            int NContinuous = detail::CountContinuous<TProcess1, detail::CountContinous<TProcess1>::count>::count>
   class ProcessSequence : public BaseProcess<ProcessSequence<TProcess1, TProcess2>> {
 
     using process1_type = typename std::decay_t<TProcess1>;
@@ -57,6 +71,9 @@ namespace corsika {
     static bool constexpr t1SwitchProcSeq = is_switch_process_sequence_v<process1_type>;
     static bool constexpr t2SwitchProcSeq = is_switch_process_sequence_v<process2_type>;
 
+    // we want to count continuous processes, to each one has a unique index
+    if constexpr (std::is_base_of_v<ContinuousProcess<typename std::decay_t<TProcess1>>   )
+
     // make sure only BaseProcess types TProcess1/2 are passed
     static_assert(std::is_base_of_v<BaseProcess<process1_type>, process1_type>,
                   "can only use process derived from BaseProcess in "
@@ -65,9 +82,6 @@ namespace corsika {
                   "can only use process derived from BaseProcess in "
                   "ProcessSequence, for Process 2");
 
-    TProcess1 A_; /// process/list A, this is a reference, if possible
-    TProcess2 B_; /// process/list B, this is a reference, if possible
-
   public:
     // resource management
     ProcessSequence() = delete; // only initialized objects
@@ -118,7 +132,7 @@ namespace corsika {
     inline void doStack(TStack& stack);
 
     template <typename TParticle, typename TTrack>
-    inline LengthType getMaxStepLength(TParticle& particle, TTrack& vTrack);
+      inline std::pair<LengthType,void*> getMaxStepLength(TParticle& particle, TTrack& vTrack);
 
     template <typename TParticle>
     inline GrammageType getInteractionLength(TParticle&& particle) {
@@ -147,6 +161,12 @@ namespace corsika {
     inline ProcessReturn selectDecay(
         TSecondaryView& view, [[maybe_unused]] InverseTimeType decay_inv_select,
         [[maybe_unused]] InverseTimeType decay_inv_sum = InverseTimeType::zero());
+
+  private:
+    int countContinuous_ = Ncontinuous;
+    TProcess1 A_; /// process/list A, this is a reference, if possible
+    TProcess2 B_; /// process/list B, this is a reference, if possible
+
   };
 
   /**

corsika/modules/ObservationPlane.hpp

@@ -31,7 +31,8 @@ namespace corsika {
                      bool = true);
 
     ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
-                               corsika::setup::Trajectory& vTrajectory);
+                               corsika::setup::Trajectory& vTrajectory,
+                               bool const stepLimit);
 
     LengthType getMaxStepLength(corsika::setup::Stack::particle_type const&,
                                 corsika::setup::Trajectory const& vTrajectory);

examples/vertical_EAS.cpp

@@ -91,8 +91,8 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
 
 int main(int argc, char** argv) {
 
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-  logging::set_level(logging::level::info);
+  corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
+  logging::set_level(logging::level::trace);
 
   CORSIKA_LOG_INFO("vertical_EAS");
 

tests/framework/testProcessSequence.cpp

@@ -247,10 +247,10 @@ struct DummyView {
   DummyData& parent() { return p_; }
 };
 
-TEST_CASE("Process Sequence", "[Process Sequence]") {
+TEST_CASE("Process Sequence General", "ProcessSequence") {
 
   logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+  corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
 
   SECTION("Check construction") {
     globalCount = 0;
@@ -373,10 +373,10 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
   }
 }
 
-TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
+TEST_CASE("Switch Process Sequence", "ProcessSequence") {
 
   logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+  corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
 
   SECTION("Check construction") {
 
@@ -479,3 +479,13 @@ TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
     CHECK(checkSec == 0);
   }
 }
+
+TEST_CASE("Continuous Process Indexing", "ProcessSequence") {
+
+  logging::set_level(logging::level::info);
+  corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
+
+  SECTION("Check construction") {
+
+  }
+}