style guide

Framework/Cascade/CMakeLists.txt

@@ -42,8 +42,9 @@ target_link_libraries (
   CORSIKArandom
   ProcessSibyll
   CORSIKAcascade
-  ProcessStackInspector
+#  ProcessStackInspector
   ProcessTrackingLine
+  ProcessNullModel
   CORSIKAstackinterface
   CORSIKAprocesses
   CORSIKAparticles

Framework/Cascade/Cascade.h

@@ -16,9 +16,19 @@
 #include <corsika/process/ProcessReturn.h>
 #include <corsika/random/RNGManager.h>
 #include <corsika/random/UniformRealDistribution.h>
-#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/stack/SecondaryView.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/setup/SetupTrajectory.h>
+
+/*  see Issue 161, we need to include SetupStack only because we need
+    to globally define StackView. This is clearly not nice and should
+    be changed, when possible. It might be that StackView needs to be
+    templated in Cascade, but this would be even worse... so we don't
+    do that until it is really needed.
+ */
+#include <corsika/setup/SetupStack.h>
+
 #include <cmath>
 #include <iostream>
 
@@ -35,25 +45,24 @@ namespace corsika::cascade {
    * it very versatile. Via the template arguments physics models are
    * plugged into the cascade simulation.
    *
-   * <b>Tracking</b> must be a class according to the
+   * <b>TTracking</b> must be a class according to the
    * TrackingInterface providing the functions: <code>void
    * Init();</code> and <code>auto GetTrack(Particle const& p)</auto>,
    * where the latter has a return type of <code>
    * geometry::Trajectory<corsika::geometry::Line or Helix> </code>
    *
-   * <b>ProcessList</b> must be a ProcessSequence.
+   * <b>TProcessList</b> must be a ProcessSequence.
    *            TimeOfIntersection(corsika::geometry::Line const& line,
    *
    * <b>Stack</b> is the storage object for particle data, i.e. with
    * Particle class type <code>Stack::ParticleType</code>
    *
    *
-
    */
 
-  template <typename Tracking, typename ProcessList, typename Stack>
+  template <typename TTracking, typename TProcessList, typename TStack>
   class Cascade {
-    using Particle = typename Stack::ParticleType;
+    using Particle = typename TStack::ParticleType;
 
     // we only want fully configured objects
     Cascade() = delete;
@@ -63,8 +72,8 @@ namespace corsika::cascade {
      * Cascade class cannot be default constructed, but needs a valid
      * list of physics processes for configuration at construct time.
      */
-    Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl,
-            Stack& stack)
+    Cascade(corsika::environment::Environment const& env, TTracking& tr, TProcessList& pl,
+            TStack& stack)
         : fEnvironment(env)
         , fTracking(tr)
         , fProcessSequence(pl)
@@ -107,15 +116,16 @@ namespace corsika::cascade {
      * New particles produced in one step are subject to further
      * processing, e.g. thinning, etc.
      */
-    void Step(Particle& particle) {
+    void Step(Particle& vParticle) {
+      using namespace corsika;
       using namespace corsika::units::si;
 
       // determine geometric tracking
-      corsika::setup::Trajectory step = fTracking.GetTrack(particle);
+      setup::Trajectory step = fTracking.GetTrack(vParticle);
 
       // determine combined total interaction length (inverse)
       InverseGrammageType const total_inv_lambda =
-          fProcessSequence.GetTotalInverseInteractionLength(particle, step);
+          fProcessSequence.GetTotalInverseInteractionLength(vParticle, step);
 
       // sample random exponential step length in grammage
       std::exponential_distribution expDist(total_inv_lambda * (1_g / (1_m * 1_m)));
@@ -126,7 +136,7 @@ namespace corsika::cascade {
 
       // convert next_step from grammage to length
       auto const* currentNode =
-          fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+          fEnvironment.GetUniverse()->GetContainingNode(vParticle.GetPosition());
 
       if (currentNode == &*fEnvironment.GetUniverse()) {
         throw std::runtime_error("particle entered void universe");
@@ -136,12 +146,12 @@ namespace corsika::cascade {
           currentNode->GetModelProperties().ArclengthFromGrammage(step, next_interact);
 
       // determine the maximum geometric step length
-      LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
+      LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
       std::cout << "distance_max=" << distance_max << std::endl;
 
       // determine combined total inverse decay time
       InverseTimeType const total_inv_lifetime =
-          fProcessSequence.GetTotalInverseLifetime(particle);
+          fProcessSequence.GetTotalInverseLifetime(vParticle);
 
       // sample random exponential decay time
       std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s);
@@ -150,9 +160,8 @@ namespace corsika::cascade {
                 << ", next_decay=" << next_decay << std::endl;
 
       // convert next_decay from time to length [m]
-      LengthType const distance_decay = next_decay * particle.GetMomentum().norm() /
-                                        particle.GetEnergy() *
-                                        corsika::units::constants::c;
+      LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
+                                        vParticle.GetEnergy() * units::constants::c;
 
       // take minimum of geometry, interaction, decay for next step
       auto const min_distance =
@@ -161,33 +170,45 @@ namespace corsika::cascade {
       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));
+      // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
+      vParticle.SetPosition(step.PositionFromArclength(min_distance));
       // .... also update time, momentum, direction, ...
+      vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
 
       step.LimitEndTo(min_distance);
 
       // particle.GetNode(); // previous VolumeNode
-      particle.SetNode(
+      vParticle.SetNode(
           currentNode); // NOTE @Max : here we need to distinguish: IF particle step is
       // limited by tracking (via fTracking.GetTrack()), THEN we need
       // to check/update VolumeNodes. In all other cases it is
       // guaranteed that we are still in the same volume
 
       // apply all continuous processes on particle + track
-      corsika::process::EProcessReturn status =
-          fProcessSequence.DoContinuous(particle, step, fStack);
+      process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
 
-      if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
-        std::cout << "Cascade: delete absorbed particle " << particle.GetPID() << " "
-                  << particle.GetEnergy() / 1_GeV << "GeV" << std::endl;
-        particle.Delete();
+      if (status == process::EProcessReturn::eParticleAbsorbed) {
+        std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
+                  << vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl;
+        vParticle.Delete();
         return;
       }
 
       std::cout << "sth. happening before geometric limit ? "
                 << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
 
+      /*
+        Create SecondaryView object on Stack. The data container
+        remains untouched and identical, and 'projectil' is identical
+        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,
+        and Decay!
+       */
+      setup::StackView secondaries(vParticle);
+      [[maybe_unused]] auto projectile = secondaries.GetProjectile();
+
       if (min_distance < distance_max) { // interaction to happen within geometric limit
         // check whether decay or interaction limits this step
 
@@ -195,33 +216,35 @@ namespace corsika::cascade {
           std::cout << "collide" << std::endl;
 
           InverseGrammageType const actual_inv_length =
-              fProcessSequence.GetTotalInverseInteractionLength(particle, step);
+              fProcessSequence.GetTotalInverseInteractionLength(vParticle, step);
 
-          corsika::random::UniformRealDistribution<InverseGrammageType> uniDist(
-              actual_inv_length);
+          random::UniformRealDistribution<InverseGrammageType> uniDist(actual_inv_length);
           const auto sample_process = uniDist(fRNG);
           InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
-          fProcessSequence.SelectInteraction(particle, step, fStack, sample_process,
+          fProcessSequence.SelectInteraction(vParticle, projectile, step, sample_process,
                                              inv_lambda_count);
         } else {
           std::cout << "decay" << std::endl;
           InverseTimeType const actual_decay_time =
-              fProcessSequence.GetTotalInverseLifetime(particle);
+              fProcessSequence.GetTotalInverseLifetime(vParticle);
 
-          corsika::random::UniformRealDistribution<InverseTimeType> uniDist(
-              actual_decay_time);
+          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);
+          fProcessSequence.SelectDecay(vParticle, projectile, sample_process,
+                                       inv_decay_count);
         }
+
+        fProcessSequence.DoSecondaries(secondaries);
+        vParticle.Delete(); // last thing in Step function
       }
     }
 
   private:
-    corsika::environment::Environment const& fEnvironment;
-    Tracking& fTracking;
-    ProcessList& fProcessSequence;
-    Stack& fStack;
+    environment::Environment const& fEnvironment;
+    TTracking& fTracking;
+    TProcessList& fProcessSequence;
+    TStack& fStack;
     corsika::random::RNG& fRNG =
         corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
   };

Framework/Cascade/testCascade.cc

@@ -16,6 +16,7 @@
 #include <corsika/cascade/Cascade.h>
 
 #include <corsika/process/ProcessSequence.h>
+#include <corsika/process/null_model/NullModel.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/process/tracking_line/TrackingLine.h>
 
@@ -75,13 +76,13 @@ public:
   ProcessSplit(HEPEnergyType e)
       : fEcrit(e) {}
 
-  template <typename Particle, typename T>
-  LengthType MaxStepLength(Particle&, T&) const {
+  template <typename Particle, typename Track>
+  LengthType MaxStepLength(Particle&, Track&) const {
     return 1_m;
   }
 
-  template <typename Particle, typename T, typename Stack>
-  EProcessReturn DoContinuous(Particle& p, T&, Stack&) {
+  template <typename Particle, typename Track>
+  EProcessReturn DoContinuous(Particle& p, Track&) {
     fCalls++;
     HEPEnergyType E = p.GetEnergy();
     if (E < fEcrit) {
@@ -115,11 +116,12 @@ TEST_CASE("Cascade", "[Cascade]") {
   auto env = MakeDummyEnv();
   tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
 
-  stack_inspector::StackInspector<setup::Stack> p0(true);
+  // stack_inspector::StackInspector<setup::Stack> stackInspect(true);
+  null_model::NullModel nullModel;
 
   const HEPEnergyType Ecrit = 85_MeV;
   ProcessSplit p1(Ecrit);
-  auto sequence = p0 << p1;
+  auto sequence = nullModel << p1;
   setup::Stack stack;
 
   cascade::Cascade EAS(env, tracking, sequence, stack);