improvements in Cascade

Framework/Cascade/Cascade.h

@@ -14,6 +14,7 @@
 
 #include <corsika/environment/Environment.h>
 #include <corsika/process/ProcessReturn.h>
+#include <corsika/random/ExponentialDistribution.h>
 #include <corsika/random/RNGManager.h>
 #include <corsika/random/UniformRealDistribution.h>
 #include <corsika/setup/SetupTrajectory.h>
@@ -21,6 +22,7 @@
 
 #include <cmath>
 #include <iostream>
+#include <limits>
 
 /**
  * The cascade namespace assembles all objects needed to simulate full particles cascades.
@@ -118,8 +120,8 @@ namespace corsika::cascade {
           fProcessSequence.GetTotalInverseInteractionLength(particle, step);
 
       // sample random exponential step length in grammage
-      std::exponential_distribution expDist(total_inv_lambda * (1_g / (1_m * 1_m)));
-      GrammageType const next_interact = (1_g / (1_m * 1_m)) * expDist(fRNG);
+      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;
@@ -133,7 +135,10 @@ namespace corsika::cascade {
       }
 
       LengthType const distance_interact =
-          currentNode->GetModelProperties().ArclengthFromGrammage(step, next_interact);
+          (next_interact == std::numeric_limits<double>::infinity() * 1_g / (1_m * 1_m))
+              ? std::numeric_limits<double>::infinity() * 1_m
+              : currentNode->GetModelProperties().ArclengthFromGrammage(step,
+                                                                        next_interact);
 
       // determine the maximum geometric step length
       LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
@@ -144,8 +149,8 @@ namespace corsika::cascade {
           fProcessSequence.GetTotalInverseLifetime(particle);
 
       // sample random exponential decay time
-      std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s);
-      TimeType const next_decay = 1_s * expDistDecay(fRNG);
+      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;