Merge branch 'master' into 'master'

Improvements in Cascade and FlatExponential

See merge request !93

Environment/FlatExponential.h

-
 /*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
  *
  * See file AUTHORS for a list of contributors.
  *
@@ -21,6 +20,7 @@
 #include <corsika/units/PhysicalUnits.h>
 
 #include <cassert>
+#include <limits>
 
 /**
  *
@@ -76,7 +76,14 @@ namespace corsika::environment {
       if (vDotA == 0) {
         return pGrammage / GetMassDensity(line.GetR0());
       } else {
-        return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
+        auto const logArg = pGrammage * vDotA / (fRho0 * fLambda) + 1;
+        if (logArg > 0) {
+          return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
+        } else {
+          return std::numeric_limits<typename decltype(
+                     pGrammage)::value_type>::infinity() *
+                 corsika::units::si::meter;
+        }
       }
     }
   };

Environment/testEnvironment.cc

@@ -75,7 +75,19 @@ TEST_CASE("FlatExponential") {
     REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
   }
 
-  SECTION("vertical") {
+  SECTION("escape grammage") {
+    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+                                 gCS, {0_m / second, 0_m / second, -5_m / second}));
+    Trajectory<Line> const trajectory(line, tEnd);
+
+    GrammageType const escapeGrammage = rho0 * lambda;
+
+    REQUIRE(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
+    REQUIRE(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
+            std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
+  }
+
+  SECTION("inclined") {
     Line const line(gOrigin, Vector<SpeedType::dimension_type>(
                                  gCS, {0_m / second, 5_m / second, 5_m / second}));
     Trajectory<Line> const trajectory(line, tEnd);
@@ -83,7 +95,6 @@ TEST_CASE("FlatExponential") {
     LengthType const length = 2 * lambda;
     GrammageType const exact =
         rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
-
     REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
     REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
   }

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;
 

Processes/Sibyll/Interaction.cc

@@ -171,7 +171,7 @@ namespace corsika::process::sibyll {
       int i = -1;
       si::CrossSectionType weightedProdCrossSection = 0_mbarn;
       // get weights of components from environment/medium
-      const auto w = mediumComposition.GetFractions();
+      const auto& w = mediumComposition.GetFractions();
       // loop over components in medium
       for (auto const targetId : mediumComposition.GetComponents()) {
         i++;