History

We need this MR to discuss details and technicalities of this branch.

The code approval procedure is described in the wiki: Code approval procedure wiki

• The MR is without WIP (work in progress) status
• Make sure the most recent CI jobs (config, quality, build_test_example) all run fine with no failures
  • if "check clang-format" failed: the code contributor has to run ./do-clang-format.py --apply eventually with the --all option
  • if "check copyright" failed the code contributor has to run./do-copyright.py --add=20xy
• Make sure also the jobs with label ready for code review succeed. This includes the optional jobs, in particular 'coverage', 'release-clang-8", "release-u-18.04" and make sure no problems occur
• Check in the "coverage" job output that the coverage did not decreases. It should always stay, or increase. If it decreased --> ask contributor to add further needed unit tests, and check coverage report.
• On the MR page, open the "Open in Web IDE" tool
  • Check if the provided solution solves the Issue, discuss on gitlab
  • Check that all changes are actually related to the issue
  • There are no debug statements left, not even commented out
  • Check all changes for coding rules and guidelines
• When all above is done:
  • Add MR label "Code Review Finished"

Framework/Cascade/Cascade.h

@@ -16,6 +16,8 @@
 #include <corsika/random/UniformRealDistribution.h>
 #include <corsika/stack/SecondaryView.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/stack/history/EventType.hpp>
+#include <corsika/stack/history/HistorySecondaryProducer.hpp>
 
 #include <corsika/setup/SetupTrajectory.h>
 
@@ -29,9 +31,7 @@
 
 #include <cassert>
 #include <cmath>
-#include <iostream>
 #include <limits>
-#include <type_traits>
 
 /**
  * The cascade namespace assembles all objects needed to simulate full particles cascades.
@@ -61,8 +61,9 @@ namespace corsika::cascade {
 
   template <typename TTracking, typename TProcessList, typename TStack,
             /*
-              TStackView is needed as template parameter because of issue 161 and the
-              inability of clang to understand "MakeView" so far.
+              TStackView is needed as explicit template parameter because
+              of issue 161 and the
+              inability of clang to understand "stack::MakeView" so far.
              */
             typename TStackView = corsika::setup::StackView>
   class Cascade {
@@ -71,6 +72,17 @@ namespace corsika::cascade {
         std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
     using MediumInterface = typename VolumeTreeNode::IModelProperties;
 
+  private:
+    // Data members
+    corsika::environment::Environment<MediumInterface> const& fEnvironment;
+    TTracking& fTracking;
+    TProcessList& fProcessSequence;
+    TStack& fStack;
+    corsika::random::RNG& fRNG =
+        corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
+    unsigned int count_ = 0;
+
+  private:
     // we only want fully configured objects
     Cascade() = delete;
 
@@ -85,9 +97,12 @@ namespace corsika::cascade {
         , fTracking(tr)
         , fProcessSequence(pl)
         , fStack(stack)
-        , energy_cut_(0 * corsika::units::si::electronvolt) {}
-
-    corsika::units::si::HEPEnergyType GetEnergyCut() const { return energy_cut_; }
+        , count_(0) {
+      C8LOG_INFO(c8_ascii_);
+#ifdef WITH_HISTORY
+      C8LOG_INFO(" - With full cascade HISTORY.");
+#endif
+    }
 
     /**
      * set the nodes for all particles on the stack according to their numerical
@@ -110,10 +125,15 @@ namespace corsika::cascade {
 
       while (!fStack.IsEmpty()) {
         while (!fStack.IsEmpty()) {
+          C8LOG_TRACE(fmt::format("Stack: {}", fStack.as_string()));
+          count_++;
           auto pNext = fStack.GetNextParticle();
-          std::cout << "========= next: " << pNext.GetPID() << std::endl;
+          C8LOG_DEBUG(fmt::format(
+              "============== next particle : count={}, pid={}, "
+              ", stack entries={}"
+              ", stack deleted={}",
+              count_, pNext.GetPID(), fStack.getEntries(), fStack.getDeleted()));
           Step(pNext);
-          std::cout << "========= stack ============" << std::endl;
           fProcessSequence.DoStack(fStack);
         }
         // do cascade equations, which can put new particles on Stack,
@@ -128,11 +148,10 @@ namespace corsika::cascade {
      * want to call forceInteraction() for the primary interaction.
      */
     void forceInteraction() {
-      std::cout << "forced interaction!" << std::endl;
+      C8LOG_DEBUG("forced interaction!");
       auto vParticle = fStack.GetNextParticle();
       TStackView secondaries(vParticle);
-      auto projectile = secondaries.GetProjectile();
-      interaction(vParticle, projectile);
+      interaction(vParticle, secondaries);
       fProcessSequence.DoSecondaries(secondaries);
       vParticle.Delete(); // todo: this should be reviewed, see below
     }
@@ -164,8 +183,11 @@ namespace corsika::cascade {
       corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
       GrammageType const next_interact = expDist(fRNG);
 
-      std::cout << "total_inv_lambda=" << total_inv_lambda
-                << ", next_interact=" << next_interact << std::endl;
+      C8LOG_DEBUG(
+          "total_lambda={} g/cm2, "
+          ", next_interact={} g/cm2",
+          double((1. / total_inv_lambda) / 1_g * 1_cm * 1_cm),
+          double(next_interact / 1_g * 1_cm * 1_cm));
 
       auto const* currentLogicalNode = vParticle.GetNode();
 
@@ -181,7 +203,7 @@ namespace corsika::cascade {
 
       // determine the maximum geometric step length from continuous processes
       LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
-      std::cout << "distance_max=" << distance_max << std::endl;
+      C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
 
       // determine combined total inverse decay time
       InverseTimeType const total_inv_lifetime =
@@ -190,8 +212,10 @@ namespace corsika::cascade {
       // sample random exponential decay time
       corsika::random::ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
       TimeType const next_decay = expDistDecay(fRNG);
-      std::cout << "total_inv_lifetime=" << total_inv_lifetime
-                << ", next_decay=" << next_decay << std::endl;
+      C8LOG_DEBUG(
+          "total_lifetime={} s"
+          ", next_decay={} s",
+          (1 / total_inv_lifetime) / 1_s, next_decay / 1_s);
 
       // convert next_decay from time to length [m]
       LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
@@ -201,7 +225,7 @@ namespace corsika::cascade {
       auto const min_distance =
           std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
 
-      std::cout << " move particle by : " << min_distance << std::endl;
+      C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
 
       // here the particle is actually moved along the trajectory to new position:
       // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
@@ -215,15 +239,14 @@ namespace corsika::cascade {
       process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
 
       if (status == process::EProcessReturn::eParticleAbsorbed) {
-        std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
-                  << vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl;
-        energy_cut_ += vParticle.GetEnergy();
+        C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
+                    vParticle.GetPID(), vParticle.GetEnergy() / 1_GeV);
         vParticle.Delete();
         return;
       }
 
-      std::cout << "sth. happening before geometric limit ? "
-                << ((min_distance < geomMaxLength) ? "yes" : "no") << std::endl;
+      C8LOG_DEBUG("sth. happening before geometric limit ? {}",
+                  ((min_distance < geomMaxLength) ? "yes" : "no"));
 
       if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
 
@@ -241,19 +264,19 @@ namespace corsika::cascade {
             to 'vParticle' above this line. However,
             projectil.AddSecondaries populate the SecondaryView, which can
             then be used afterwards for further processing. Thus: it is
-            important to use projectle (and not vParticle) for Interaction,
+            important to use projectle/view (and not vParticle) for Interaction,
             and Decay!
           */
 
           [[maybe_unused]] auto projectile = secondaries.GetProjectile();
 
           if (min_distance == distance_interact) {
-            interaction(vParticle, projectile);
+            interaction(vParticle, secondaries);
           } else {
             assert(min_distance == distance_decay);
-            decay(vParticle, projectile);
+            decay(vParticle, secondaries);
             // make sure particle actually did decay if it should have done so
-            if (secondaries.GetSize() == 1 &&
+            if (secondaries.getSize() == 1 &&
                 projectile.GetPID() == secondaries.GetNextParticle().GetPID())
               throw std::runtime_error(
                   fmt::format("Cascade: {} decayed into itself!",
@@ -261,15 +284,10 @@ namespace corsika::cascade {
           }
 
           fProcessSequence.DoSecondaries(secondaries);
-          vParticle.Delete(); // todo: this should be reviewed. Where
-                              // exactly are particles best deleted, and
-                              // where they should NOT be
-                              // deleted... maybe Delete function should
-                              // be "protected" and not accessible to physics
+          vParticle.Delete();
 
         } else { // step-length limitation within volume
-
-          std::cout << "step-length limitation" << std::endl;
+          C8LOG_DEBUG("step-length limitation");
           // no extra physics happens here. just proceed to next step.
         }
 
@@ -283,10 +301,7 @@ namespace corsika::cascade {
         };
 
         assert(assertion()); // numerical and logical nodes don't match
-
-      } else { // boundary crossing, step is limited by volume boundary
-
-        std::cout << "boundary crossing! next node = " << nextVol << std::endl;
+      } else {               // boundary crossing, step is limited by volume boundary
         vParticle.SetNode(nextVol);
         /*
           DoBoundary may delete the particle (or not)
@@ -299,9 +314,8 @@ namespace corsika::cascade {
       }
     }
 
-    auto decay(Particle& particle,
-               decltype(std::declval<TStackView>().GetProjectile()) projectile) {
-      std::cout << "decay" << std::endl;
+    auto decay(Particle& particle, TStackView& view) {
+      C8LOG_DEBUG("decay");
       units::si::InverseTimeType const actual_decay_time =
           fProcessSequence.GetTotalInverseLifetime(particle);
 
@@ -309,13 +323,14 @@ namespace corsika::cascade {
           actual_decay_time);
       const auto sample_process = uniDist(fRNG);
       units::si::InverseTimeType inv_decay_count = units::si::InverseTimeType::zero();
-      return fProcessSequence.SelectDecay(particle, projectile, sample_process,
-                                          inv_decay_count);
+      auto const returnCode =
+          fProcessSequence.SelectDecay(particle, view, sample_process, inv_decay_count);
+      SetEventType(view, history::EventType::Decay);
+      return returnCode;
     }
 
-    auto interaction(Particle& particle,
-                     decltype(std::declval<TStackView>().GetProjectile()) projectile) {
-      std::cout << "collide" << std::endl;
+    auto interaction(Particle& particle, TStackView& view) {
+      C8LOG_DEBUG("collide");
 
       units::si::InverseGrammageType const current_inv_length =
           fProcessSequence.GetTotalInverseInteractionLength(particle);
@@ -324,20 +339,30 @@ namespace corsika::cascade {
           current_inv_length);
       const auto sample_process = uniDist(fRNG);
       auto inv_lambda_count = units::si::InverseGrammageType::zero();
-      return fProcessSequence.SelectInteraction(particle, projectile, sample_process,
-                                                inv_lambda_count);
+      auto const returnCode = fProcessSequence.SelectInteraction(
+          particle, view, sample_process, inv_lambda_count);
+      SetEventType(view, history::EventType::Interaction);
+      return returnCode;
     }
 
-  private:
-    corsika::environment::Environment<MediumInterface> const& fEnvironment;
-    TTracking& fTracking;
-    TProcessList& fProcessSequence;
-    TStack& fStack;
-    corsika::random::RNG& fRNG =
-        corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
-
-    corsika::units::si::HEPEnergyType energy_cut_;
+    void SetEventType(TStackView& view, [[maybe_unused]] history::EventType eventType) {
+      if constexpr (TStackView::has_event) {
+        for (auto&& sec : view) { sec.GetEvent()->setEventType(eventType); }
+      }
+    }
 
-  }; // namespace corsika::cascade
+    // but this here temporarily. Should go into dedicated file later:
+    const char* c8_ascii_ =
+        R"V0G0N(
+  ,ad8888ba,     ,ad8888ba,    88888888ba    ad88888ba   88  88      a8P          db              ad88888ba   
+ d8"'    `"8b   d8"'    `"8b   88      "8b  d8"     "8b  88  88    ,88'          d88b            d8"     "8b  
+d8'            d8'        `8b  88      ,8P  Y8,          88  88  ,88"           d8'`8b           Y8a     a8P  
+88             88          88  88aaaaaa8P'  `Y8aaaaa,    88  88,d88'           d8'  `8b           "Y8aaa8P"   
+88             88          88  88""""88'      `"""""8b,  88  8888"88,         d8YaaaaY8b          ,d8"""8b,   
+Y8,            Y8,        ,8P  88    `8b            `8b  88  88P   Y8b       d8""""""""8b        d8"     "8b  
+ Y8a.    .a8P   Y8a.    .a8P   88     `8b   Y8a     a8P  88  88     "88,    d8'        `8b       Y8a     a8P  
+  `"Y8888Y"'     `"Y8888Y"'    88      `8b   "Y88888P"   88  88       Y8b  d8'          `8b       "Y88888P"
+	)V0G0N";
+  };
 
 } // namespace corsika::cascade