diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h index 7079a22bfa39f4b061e1bfe111553c34e595d8c7..19e42c725c4b5204dc9c6fc166f9be00f44bdb7a 100644 --- a/Environment/FlatExponential.h +++ b/Environment/FlatExponential.h @@ -1,6 +1,5 @@ - /* - * (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; + } } } }; diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc index 5beab6587a20199a24bb8946b96ccedfe4b8d293..206849e7eb17b6d33a23fffddb4abefbb16900e8 100644 --- a/Environment/testEnvironment.cc +++ b/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)); } diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index c55ac0b43b86e1956750618ec65e868a0d1334e0..1a559a210eb16c60f3d8165b0b18de5589dbdbd3 100644 --- a/Framework/Cascade/Cascade.h +++ b/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; diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc index 06d301e3a9f4d264e3b6b56c969fd8d3742269cc..b16d48332a76aaa84edd403ed1a3736487a0269e 100644 --- a/Processes/Sibyll/Interaction.cc +++ b/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++;