From b54178905782da0e44a03302ebcbd9e4c56ceda1 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Sat, 23 Feb 2019 15:18:20 +0100
Subject: [PATCH] improvements in Cascade
---
Framework/Cascade/Cascade.h | 15 ++++++++++-----
1 file changed, 10 insertions(+), 5 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index c55ac0b43..1a559a210 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;
--
GitLab