diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index aa52dd3f38b47442e0c1d3697735b2306dfc0d22..a52219a52896ab6382d6ab72148da5288d041124 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -16,6 +16,7 @@ variables:
# _alternatively_ corsika-data can be downloaded as submodule:
GIT_SUBMODULE_STRATEGY: normal # none: we get the submodules in before_script,
# normal: get submodules automatically
+ CTEST_OUTPUT_ON_FAILURE: 1
#
@@ -251,6 +252,8 @@ test-clang-8:
artifacts:
when: always
expire_in: 3 days
+ paths:
+ - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
reports:
junit:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
@@ -311,6 +314,7 @@ build_test-clang-8:
- if: $CI_COMMIT_TAG
when: manual
allow_failure: true
+ allow_failure: true
artifacts:
when: always
expire_in: 3 days
@@ -319,6 +323,7 @@ build_test-clang-8:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
paths:
- ${CI_PROJECT_DIR}/build/build_examples/examples.log.gz
+ - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
cache:
paths:
- ${CI_PROJECT_DIR}/build/
@@ -440,6 +445,7 @@ install-clang-8:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
paths:
- ${CI_PROJECT_DIR}/build/build_examples/examples.log.gz
+ - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
# release for gcc
release-full-u-18_04:
diff --git a/CMakeModules/FindCONEX.cmake b/cmake/FindCONEX.cmake
similarity index 100%
rename from CMakeModules/FindCONEX.cmake
rename to cmake/FindCONEX.cmake
diff --git a/conanfile.txt b/conanfile.txt
index 85c1238b40b290ca181660b0a8bda665b1c88c5b..78cda273ada5bf3c944cb45f6f22df53d0567a8c 100644
--- a/conanfile.txt
+++ b/conanfile.txt
@@ -1,9 +1,28 @@
[requires]
-spdlog/1.8.0
+readline/8.0 # this is only needed to fix a missing dependency in "bison" which is pulled-in by arrow
+spdlog/1.8.5
catch2/2.13.3
eigen/3.3.8
-boost/1.74.0
+boost/1.76.0
zlib/1.2.11
+proposal/7.0.2
+yaml-cpp/0.6.3
+arrow/2.0.0
[generators]
cmake
+
+[options]
+readline:shared=True
+arrow:shared=False
+arrow:parquet=True
+arrow:fPIC=False
+arrow:with_re2=False
+arrow:with_protobuf=False
+arrow:with_openssl=False
+arrow:with_gflags=False
+arrow:with_glog=False
+arrow:with_grpc=False
+arrow:with_utf8proc=False
+arrow:with_zstd=False
+arrow:with_bz2=False
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 28306c590bb7da30d19d70f7ddfa19917ce525a0..32127908c714d79dfdd7e22e74615aed449b1b7c 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -60,7 +60,7 @@ namespace corsika {
setNodes();
auto vParticle = stack_.getNextParticle();
TStackView secondaries(vParticle);
- interaction(secondaries);
+ interaction(secondaries, sequence_.getInverseInteractionLength(vParticle));
sequence_.doSecondaries(secondaries);
vParticle.erase(); // primary particle is done
}
@@ -243,15 +243,15 @@ namespace corsika {
[[maybe_unused]] auto projectile = secondaries.getProjectile();
if (distance_interact < distance_decay) {
- interaction(secondaries);
+ interaction(secondaries, total_inv_lambda);
} else {
- decay(secondaries);
- // make sure particle actually did decay if it should have done so
- // \todo this should go to a validation code and not be included here
- if (secondaries.getSize() == 1 &&
- projectile.getPID() == secondaries.getNextParticle().getPID())
- throw std::runtime_error(fmt::format("Particle {} decays into itself!",
- get_name(projectile.getPID())));
+ if (decay(secondaries, total_inv_lifetime) == ProcessReturn::Decayed) {
+ if (secondaries.getSize() == 1 &&
+ projectile.getPID() == secondaries.getNextParticle().getPID()) {
+ throw std::runtime_error(fmt::format("Particle {} decays into itself!",
+ get_name(projectile.getPID())));
+ }
+ }
}
sequence_.doSecondaries(secondaries);
@@ -261,16 +261,29 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
inline ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::decay(
- TStackView& view) {
+ TStackView& view, InverseTimeType initial_inv_decay_time) {
CORSIKA_LOG_DEBUG("decay");
+
+#ifdef DEBUG
InverseTimeType const actual_decay_time = sequence_.getInverseLifetime(view.parent());
+ if (actual_decay_time * 0.99 > initial_inv_decay_time) {
+ CORSIKA_LOG_WARN(
+ "Decay time decreased during step! This leads to un-physical step length. "
+ "delta_inverse_decay_time={}",
+ 1 / initial_inv_decay_time - 1 / actual_decay_time);
+ }
+#endif
- UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
+ // one option is that decay_time is now larger (less
+ // probability for decay) than it was before the step, thus,
+ // no decay might actually occur and is allowed
+
+ UniformRealDistribution<InverseTimeType> uniDist(initial_inv_decay_time);
const auto sample_process = uniDist(rng_);
auto const returnCode = sequence_.selectDecay(view, sample_process);
if (returnCode != ProcessReturn::Decayed) {
- CORSIKA_LOG_WARN("Particle did not decay!");
+ CORSIKA_LOG_DEBUG("Particle did not decay!");
}
setEventType(view, history::EventType::Decay);
return returnCode;
@@ -279,18 +292,31 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
inline ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::interaction(
- TStackView& view) {
+ TStackView& view, InverseGrammageType initial_inv_int_length) {
+
CORSIKA_LOG_DEBUG("collide");
- InverseGrammageType const current_inv_length =
- sequence_.getInverseInteractionLength(view.parent());
+#ifdef DEBUG
+ InverseGrammageType const actual_inv_length = sequence_.getInverseInteractionLength(
+ view.parent()); // 1/lambda_int after step, -dE/dX etc.
+
+ if (actual_inv_length * 0.99 > initial_inv_int_length) {
+ CORSIKA_LOG_WARN(
+ "Interaction length decreased during step! This leads to un-physical step "
+ "length. delta_innverse_interaction_length={}",
+ 1 / initial_inv_int_length - 1 / actual_inv_length);
+ }
+#endif
- UniformRealDistribution<InverseGrammageType> uniDist(current_inv_length);
+ // one option is that interaction_length is now larger (less
+ // probability for collision) than it was before the step, thus,
+ // no interaction might actually occur and is allowed
+ UniformRealDistribution<InverseGrammageType> uniDist(initial_inv_int_length);
const auto sample_process = uniDist(rng_);
auto const returnCode = sequence_.selectInteraction(view, sample_process);
if (returnCode != ProcessReturn::Interacted) {
- CORSIKA_LOG_WARN("Particle did not interace!");
+ CORSIKA_LOG_DEBUG("Particle did not interact!");
}
setEventType(view, history::EventType::Interaction);
return returnCode;
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 6c9bd789289a6fa24204bfe10df34e5413b5b6fa..5eb3d20ae785b7645101b362cb4697282a29ecdd 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -221,7 +221,7 @@ namespace corsika {
IndexProcess2>::doStack(TStack& stack) {
if constexpr (is_process_v<process1_type>) { // to protect from further compiler
// errors if process1_type is invalid
- if constexpr (std::is_base_of_v<StackProcess<process1_type>, process1_type> ||
+ if constexpr (is_stack_process_v<process1_type> ||
(process1_type::is_process_sequence &&
!process1_type::is_switch_process_sequence)) {
@@ -237,7 +237,7 @@ namespace corsika {
}
if constexpr (is_process_v<process2_type>) { // to protect from further compiler
// errors if process2_type is invalid
- if constexpr (std::is_base_of_v<StackProcess<process2_type>, process2_type> ||
+ if constexpr (is_stack_process_v<process2_type> ||
(process2_type::is_process_sequence &&
!process2_type::is_switch_process_sequence)) {
@@ -316,14 +316,14 @@ namespace corsika {
if constexpr (is_process_v<process1_type>) { // to protect from further compiler
// errors if process1_type is invalid
- if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, process1_type> ||
+ if constexpr (is_interaction_process_v<process1_type> ||
process1_type::is_process_sequence) {
tot += A_.getInverseInteractionLength(particle);
}
}
if constexpr (is_process_v<process2_type>) { // to protect from further compiler
// errors if process2_type is invalid
- if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, process2_type> ||
+ if constexpr (is_interaction_process_v<process2_type> ||
process2_type::is_process_sequence) {
tot += B_.getInverseInteractionLength(particle);
}
@@ -350,9 +350,7 @@ namespace corsika {
A_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
// if A_ did succeed, stop routine. Not checking other static branch B_.
if (ret != ProcessReturn::Ok) { return ret; }
- } else if constexpr (is_process_v<process1_type> &&
- std::is_base_of_v<InteractionProcess<process1_type>,
- process1_type>) {
+ } else if constexpr (is_interaction_process_v<process1_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
@@ -376,9 +374,7 @@ namespace corsika {
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
- } else if constexpr (is_process_v<process2_type> &&
- std::is_base_of_v<InteractionProcess<process2_type>,
- process2_type>) {
+ } else if constexpr (is_interaction_process_v<process2_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
// soon as SecondaryView::parent() is migrated!
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index 2f87ade6d329ab249f5a61f21b79c72d07d315fd..634395f55e0a59524263df2f9d72ebe80cac8be6 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -36,13 +36,11 @@ namespace corsika {
typename TParticle::node_type const& from,
typename TParticle::node_type const& to) {
if (select_(particle)) {
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
- process1_type> ||
+ if constexpr (is_boundary_process_v<process1_type> ||
process1_type::is_process_sequence) {
// interface checking on TSequence
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
- process1_type>) {
+ if constexpr (is_boundary_process_v<process1_type>) {
static_assert(
has_method_doBoundaryCrossing_v<TSequence, ProcessReturn, TParticle&>,
@@ -56,13 +54,11 @@ namespace corsika {
}
} else {
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>,
- process2_type> ||
+ if constexpr (is_boundary_process_v<process2_type> ||
process2_type::is_process_sequence) {
// interface checking on USequence
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>,
- process2_type>) {
+ if constexpr (is_boundary_process_v<process2_type>) {
static_assert(
has_method_doBoundaryCrossing_v<USequence, ProcessReturn, TParticle>,
@@ -136,7 +132,7 @@ namespace corsika {
IndexProcess2>::doSecondaries(TSecondaries& vS) {
const auto& particle = vS.parent();
if (select_(particle)) {
- if constexpr (std::is_base_of_v<SecondariesProcess<process1_type>, process1_type> ||
+ if constexpr (is_secondaries_process_v<process1_type> ||
process1_type::is_process_sequence) {
// interface checking on TSequence
@@ -150,7 +146,7 @@ namespace corsika {
A_.doSecondaries(vS);
}
} else {
- if constexpr (std::is_base_of_v<SecondariesProcess<process2_type>, process2_type> ||
+ if constexpr (is_secondaries_process_v<process2_type> ||
process2_type::is_process_sequence) {
// interface checking on USequence
@@ -219,14 +215,14 @@ namespace corsika {
IndexProcess2>::getInverseInteractionLength(TParticle&& particle) {
if (select_(particle)) {
- if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, process1_type> ||
+ if constexpr (is_interaction_process_v<process1_type> ||
process1_type::is_process_sequence) {
return A_.getInverseInteractionLength(particle);
}
} else {
- if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, process2_type> ||
+ if constexpr (is_interaction_process_v<process2_type> ||
process2_type::is_process_sequence) {
return B_.getInverseInteractionLength(particle);
}
@@ -249,8 +245,7 @@ namespace corsika {
A_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
// if A_ did succeed, stop routine. Not checking other static branch B_.
if (ret != ProcessReturn::Ok) { return ret; }
- } else if constexpr (std::is_base_of_v<InteractionProcess<process1_type>,
- process1_type>) {
+ } else if constexpr (is_interaction_process_v<process1_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
@@ -272,8 +267,7 @@ namespace corsika {
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
- } else if constexpr (std::is_base_of_v<InteractionProcess<process2_type>,
- process2_type>) {
+ } else if constexpr (is_interaction_process_v<process2_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
@@ -301,14 +295,14 @@ namespace corsika {
IndexProcess2>::getInverseLifetime(TParticle&& particle) {
if (select_(particle)) {
- if constexpr (std::is_base_of_v<DecayProcess<process1_type>, process1_type> ||
+ if constexpr (is_decay_process_v<process1_type> ||
process1_type::is_process_sequence) {
return A_.getInverseLifetime(particle);
}
} else {
- if constexpr (std::is_base_of_v<DecayProcess<process2_type>, process2_type> ||
+ if constexpr (is_decay_process_v<process2_type> ||
process2_type::is_process_sequence) {
return B_.getInverseLifetime(particle);
}
@@ -331,8 +325,7 @@ namespace corsika {
ProcessReturn const ret = A_.selectDecay(view, decay_inv_select, decay_inv_sum);
// if A_ did succeed, stop routine here (not checking other static branch B_)
if (ret != ProcessReturn::Ok) { return ret; }
- } else if constexpr (std::is_base_of_v<DecayProcess<process1_type>,
- process1_type>) {
+ } else if constexpr (is_decay_process_v<process1_type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A_.getInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
@@ -355,8 +348,7 @@ namespace corsika {
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
return B_.selectDecay(view, decay_inv_select, decay_inv_sum);
- } else if constexpr (std::is_base_of_v<DecayProcess<process2_type>,
- process2_type>) {
+ } else if constexpr (is_decay_process_v<process2_type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B_.getInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
new file mode 100644
index 0000000000000000000000000000000000000000..3a68a3d6ebcf39290807c3402c09fcfccfe8ab6e
--- /dev/null
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -0,0 +1,298 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <cmath>
+
+namespace corsika {
+
+ namespace andre {
+
+ //---------------------------------------------------------------------------
+ // solve cubic equation A x^3 + B*x^2 + C*x + D = 0
+ inline std::vector<double> solve_cubic_real_analytic(long double A, long double B,
+ long double C, long double D,
+ double const epsilon) {
+
+ if (std::abs(A) < epsilon) { return solve_quadratic_real(B, C, epsilon); }
+
+ long double a = B / A;
+ long double b = C / A;
+ long double c = D / A;
+
+ long double a2 = a * a;
+ long double q = (a2 - 3 * b) / 9;
+ long double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
+ long double q3 = q * q * q;
+
+ // disc = q**3 + r**2
+ // w:e = r*r exactly
+ long double w = r * r;
+ long double e = std::fma(r, r, -w);
+ // s:t = q*q exactly
+ long double s = -q * q;
+ long double t = std::fma(-q, q, -s);
+ // s:t * q + w:e = s*q + w + t*q +e = s*q+w + u:v + e = f + u:v + e
+ long double f = std::fma(s, q, w);
+ long double u = t * q;
+ long double v = std::fma(t, q, -u);
+ // error-free sum f+u. See Knuth, TAOCP vol. 2
+ long double uf = u + f;
+ long double au = uf - u;
+ long double ab = uf - au;
+ au = f - au;
+ ab = u - ab;
+ // sum all terms into final result
+ long double const disc = -(((e + uf) + au) + ab) + v;
+
+ if (disc >= 0) {
+ long double t = r / std::sqrt(q3);
+ if (t < -1) t = -1;
+ if (t > 1) t = 1;
+ t = std::acos(t);
+ a /= 3;
+ q = -2 * std::sqrt(q);
+ return {double(q * std::cos(t / 3) - a),
+ double(q * std::cos((t + 2 * M_PI) / 3) - a),
+ double(q * std::cos((t - 2 * M_PI) / 3) - a)};
+ } else {
+ long double term1 = -cbrt(std::fabs(r) + std::sqrt(-disc));
+ if (r < 0) term1 = -term1;
+ long double term2 = (0 == term1 ? 0 : q / term1);
+
+ a /= 3;
+ long double test = 0.5 * std::sqrt(3.) * (term1 - term2);
+ if (std::fabs(test) < epsilon) {
+ return {double((term1 + term2) - 1), double(-0.5 * (term1 + term2) - a)};
+ }
+
+ return {double((term1 + term2) - a)};
+ }
+ }
+ } // namespace andre
+
+ inline std::vector<double> solve_cubic_depressed_disciminant_real(
+ long double p, long double q, long double const disc, double const epsilon) {
+
+ CORSIKA_LOG_TRACE("p={}, q={}, disc={}", p, q, disc);
+
+ if (std::abs(disc) < epsilon) { // disc==0 multiple roots !
+ if (std::abs(p) < epsilon) { // tripple root
+ return {0};
+ }
+ // double root, single root
+ CORSIKA_LOG_TRACE("cubic double root");
+ return {double(-3 * q / (2 * p)), double(3 * q / p)};
+ }
+
+ if (std::abs(p) < epsilon) { // p==0 --> x^3 + q = 0
+ return {double(std::cbrt(-q))};
+ }
+
+ if (disc > 0) { // three real roots
+ CORSIKA_LOG_TRACE("cubic three roots");
+ long double const p_third = p / 3;
+ long double const sqrt_minus_p_third = std::sqrt(-p_third);
+
+ long double const arg = std::acos(q / (2 * p_third) / sqrt_minus_p_third) / 3;
+
+ long double constexpr pi = M_PI;
+ return {double(2 * sqrt_minus_p_third * std::cos(arg)),
+ double(2 * sqrt_minus_p_third * std::cos(arg - 2 * pi / 3)),
+ double(2 * sqrt_minus_p_third * std::cos(arg - 4 * pi / 3))};
+ }
+
+ // thus disc < 0 -> one real root
+ if (p < 0) {
+ CORSIKA_LOG_TRACE("cubic cosh");
+ long double const abs_q = std::abs(q);
+ long double const p_third = p / 3;
+ long double const sqrt_minus_p_third = std::sqrt(-p_third);
+ CORSIKA_LOG_TRACE("sqrt_minus_p_third={}, arcosh({})={}", sqrt_minus_p_third,
+ -abs_q / (2 * p_third) / sqrt_minus_p_third,
+ std::acosh(-abs_q / (2 * p_third) / sqrt_minus_p_third));
+ CORSIKA_LOG_TRACE(
+ "complex: {}",
+ -2 * abs_q / q * sqrt_minus_p_third *
+ std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3));
+ double const z =
+ -2 * abs_q / q * sqrt_minus_p_third *
+ std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3);
+ return {z};
+ } else { // p > 0
+ CORSIKA_LOG_TRACE("cubic sinh");
+ long double const p_third = p / 3;
+ long double const sqrt_p_third = std::sqrt(p_third);
+ return {double(-2 * sqrt_p_third *
+ std::sinh(std::asinh(q / (2 * p_third * sqrt_p_third)) / 3))};
+ }
+ }
+
+ inline std::vector<double> solve_cubic_depressed_real(long double p, long double q,
+ double const epsilon) {
+
+ // thanks!
+ // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic
+ // long double const disc = -(4 * p * p * p + 27 * q * q);
+ // disc = (p/3)**3 + (q/2)**2
+ long double p_third = p / 3;
+ long double q_half = q / 2;
+ // w:e = (q/2)*(q/2) exactly
+ long double w = q_half * q_half;
+ long double e = std::fma(q_half, q_half, -w);
+ // s:t = (p/3)*(p/3) exactly
+ long double s = p_third * p_third;
+ long double t = std::fma(p_third, p_third, -s);
+ // s:t * p + w:e = s*p + w + t*p +e = s*p+w + u:v + e = f + u:v + e
+ long double f = std::fma(s, p_third, w);
+ long double u = t * p_third;
+ long double v = std::fma(t, p_third, -u);
+ // error-free sum f+u. See Knuth, TAOCP vol. 2
+ long double a = u + f;
+ long double b = a - u;
+ long double c = a - b;
+ b = f - b;
+ c = u - c;
+ // sum all terms into final result
+ long double const disc = -(((e + a) + b) + c) + v;
+ return solve_cubic_depressed_disciminant_real(p, q, disc, epsilon);
+ }
+
+ /**
+ * Analytical approach. Not very stable in all conditions.
+ */
+ inline std::vector<double> solve_cubic_real_analytic(long double a, long double b,
+ long double c, long double d,
+ double const epsilon) {
+
+ CORSIKA_LOG_TRACE("cubic: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b, c,
+ d, epsilon, (std::abs(a - 1) < epsilon), (std::abs(b) < epsilon));
+
+ if (std::abs(a) < epsilon) { // this is just a quadratic
+ return solve_quadratic_real(b, c, d, epsilon);
+ }
+
+ if ((std::abs(a - 1) < epsilon) &&
+ (std::abs(b) < epsilon)) { // this is a depressed cubic
+ return solve_cubic_depressed_real(c, d, epsilon);
+ }
+
+ // p = (3ac - b^2) / (3a^2) = 3 * ( c/(3*a) - b^2/(9*a^2) )
+ long double b_over_a = b / a;
+ long double const p_third = std::fma(-b_over_a, b_over_a / 9, c / (a * 3));
+
+ // p = std::fma(a * 3, c, -b2) / (3 * a2);
+ // q = (2*b^3 - 9*abc + 27*a^2*d) / (27a^3) = 2 * ( b^3/(27*a^3) - bc/(6*a^2) +
+ // d/(2*a) )
+ long double const q_half_term1 = std::fma(b_over_a, b_over_a / 27, -c / (a * 6));
+ long double const q_half = std::fma(b_over_a, q_half_term1, d / (a * 2));
+
+ std::vector<double> zs = solve_cubic_depressed_real(3 * p_third, 2 * q_half, epsilon);
+ CORSIKA_LOG_TRACE("cubic: solve_depressed={}, b/3a={}", fmt::join(zs, ", "),
+ b / (a * 3));
+ for (auto& z : zs) {
+ z -= b / (a * 3);
+ double const b1 = z + b / a;
+ double const b0 = b1 * z + c / a;
+ std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
+ CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
+ static_pow<3>(z) * a + static_pow<2>(z) * b + z * c + d);
+ }
+ CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(zs, ", "));
+ return zs;
+ }
+
+ template <typename T> // T must be floating point type
+ inline T cubic_function(T x, T a, T b, T c, T d) {
+ T x2 = x * x;
+ return x2 * x * a + x2 * b + x * c + d;
+ }
+
+ template <typename T> // T must be floating point type
+ inline T cubic_function_dfdx(T x, T a, T b, T c) {
+ T x2 = x * x;
+ return x2 * a * 3 + x * b * 2 + c;
+ }
+
+ template <typename T> // T must be floating point type
+ inline T cubic_function_d2fd2x(T x, T a, T b) {
+ return x * a * 6 + b * 2;
+ }
+
+ /**
+ * Iterative approach.
+ */
+
+ inline std::vector<double> solve_cubic_real(long double a, long double b, long double c,
+ long double d, double const epsilon) {
+
+ CORSIKA_LOG_TRACE("cubic_2: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b,
+ c, d, epsilon, (std::abs(a - 1) < epsilon),
+ (std::abs(b) < epsilon));
+
+ if (std::abs(a) < epsilon) { // this is just a quadratic
+ return solve_quadratic_real(b, c, d, epsilon);
+ }
+
+ long double const dist = std::fma(b / a, b / a, -3 * c / a);
+ long double const xinfl = -b / (a * 3);
+
+ long double x1 = xinfl;
+ long double f_x1 = cubic_function(xinfl, a, b, c, d);
+
+ CORSIKA_LOG_TRACE("dist={} xinfl={} f_x1={} {} {}", dist, xinfl, f_x1,
+ cubic_function(xinfl - 2 / 3 * std::sqrt(dist), a, b, c, d),
+ cubic_function(xinfl + 2 / 3 * std::sqrt(dist), a, b, c, d));
+
+ if (std::abs(f_x1) > epsilon) {
+ if (std::abs(dist) < epsilon) {
+ x1 = xinfl - std::cbrt(f_x1);
+ } else if (dist > 0) {
+ if (f_x1 > 0)
+ x1 = xinfl - 2 / 3 * std::sqrt(dist);
+ else
+ x1 = xinfl + 2 / 3 * std::sqrt(dist);
+ }
+
+ int niter = 0;
+ const int maxiter = 100;
+ do {
+ long double const f_prime_x1 = cubic_function_dfdx(x1, a, b, c);
+ long double const f_prime2_x1 = cubic_function_d2fd2x(x1, a, b);
+ // if (potential) saddle point... avoid
+ if (std::abs(f_prime_x1) < epsilon) {
+ x1 -= std::cbrt(f_x1);
+ } else {
+ x1 -=
+ f_x1 * f_prime_x1 / (static_pow<2>(f_prime_x1) - 0.5 * f_x1 * f_prime2_x1);
+ }
+ f_x1 = cubic_function(x1, a, b, c, d);
+ CORSIKA_LOG_TRACE("niter={} x1={} f_x1={} f_prime={} f_prime2={} eps={}", niter,
+ x1, f_x1, f_prime_x1, f_prime2_x1, epsilon);
+ } while ((++niter < maxiter) && (std::abs(f_x1) > epsilon));
+
+ CORSIKA_LOG_TRACE("niter={}", niter);
+ }
+
+ CORSIKA_LOG_TRACE("x1={} f_x1={}", x1, f_x1);
+
+ double const b1 = x1 + b / a;
+ double const b0 = b1 * x1 + c / a;
+ std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
+ CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
+ cubic_function(x1, a, b, c, d));
+
+ quad_check.push_back(x1);
+ CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(quad_check, ", "));
+ return quad_check;
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/framework/utility/LinearSolver.inl b/corsika/detail/framework/utility/LinearSolver.inl
new file mode 100644
index 0000000000000000000000000000000000000000..6eccee91ce4614d5cf95cd0948c6497824741213
--- /dev/null
+++ b/corsika/detail/framework/utility/LinearSolver.inl
@@ -0,0 +1,35 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <complex>
+
+namespace corsika {
+
+ inline std::vector<double> solve_linear_real(double a, double b) {
+
+ if (a == 0) {
+ return {}; // no (b!=0), or infinite number (b==0) of solutions....
+ }
+
+ return {-b / a};
+ }
+
+ inline std::vector<std::complex<double>> solve_linear(double a, double b) {
+
+ if (std::abs(a) == 0) {
+ return {}; // no (b!=0), or infinite number (b==0) of solutions....
+ }
+
+ return {std::complex<double>(-b / a, 0)};
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/framework/utility/QuadraticSolver.inl b/corsika/detail/framework/utility/QuadraticSolver.inl
new file mode 100644
index 0000000000000000000000000000000000000000..9c9759e975ae6d322dee7a478f935c095bdf0ff9
--- /dev/null
+++ b/corsika/detail/framework/utility/QuadraticSolver.inl
@@ -0,0 +1,83 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+ inline std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
+ long double c,
+ double const epsilon) {
+ if (std::abs(a) < epsilon) { return solve_linear(b, c); }
+
+ if (std::abs(c) < epsilon) {
+ std::vector<std::complex<double>> lin_result = solve_linear(a, b);
+ lin_result.push_back({0.});
+ return lin_result;
+ }
+
+ long double const radicant = static_pow<2>(b) - a * c * 4;
+
+ if (radicant < -epsilon) { // two complex solutions
+ double const rpart = -b / 2 * a;
+ double const ipart = std::sqrt(-radicant);
+ return {{rpart, ipart}, {rpart, -ipart}};
+ }
+
+ if (radicant < epsilon) { // just one real solution
+ return {{double(-b / 2 * a), 0}};
+ }
+
+ // two real solutions, use Citardauq formula and actively avoid cancellation in the
+ // denominator
+
+ const long double x1 =
+ 2 * c / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant));
+
+ return {{double(x1), 0}, {double(c / (a * x1)), 0}};
+ }
+
+ inline std::vector<double> solve_quadratic_real(long double a, long double b,
+ long double c, double const epsilon) {
+
+ CORSIKA_LOG_TRACE("quadratic: a={} b={} c={}", a, b, c);
+
+ if (std::abs(a) < epsilon) { return solve_linear_real(b, c); }
+ if (std::abs(c) < epsilon) {
+ std::vector<double> lin_result = solve_linear_real(a, b);
+ lin_result.push_back(0.);
+ return lin_result;
+ }
+
+ // long double const radicant = std::pow(b, 2) - a * c * 4;
+ // Thanks!
+ // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic
+ long double w = a * 4 * c;
+ long double e = std::fma(a * 4, c, -w);
+ long double f = std::fma(b, b, -w);
+ long double radicant = f + e;
+
+ CORSIKA_LOG_TRACE(" radicant={} {} ", radicant, b * b - a * c * 4);
+
+ if (std::abs(radicant) < epsilon) { // just one real solution
+ return {double(-b / (2 * a))};
+ }
+
+ if (radicant < 0) { return {}; }
+
+ // two real solutions, use Citardauq formula and actively avoid cancellation in the
+ // denominator
+
+ long double const x1 =
+ c * 2 / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant));
+
+ CORSIKA_LOG_TRACE("quadratic x1={} x2={}", double(x1), double(c / (a * x1)));
+
+ return {double(x1), double(c / (a * x1))};
+ }
+} // namespace corsika
diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
index a82abf7ac022c3c99f6e342f88eb35f2ddaa6425..988bbe80b315d59344c8c71f8960ac5833e270aa 100644
--- a/corsika/detail/framework/utility/QuarticSolver.inl
+++ b/corsika/detail/framework/utility/QuarticSolver.inl
@@ -8,128 +8,154 @@
#pragma once
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CubicSolver.hpp>
#include <cmath>
-namespace corsika::quartic_solver {
-
- //---------------------------------------------------------------------------
- // solve cubic equation x^3 + a*x^2 + b*x + c = 0
- // x - array of size 3
- // In case 3 real roots: => x[0], x[1], x[2], return 3
- // 2 real roots: x[0], x[1], return 2
- // 1 real root : x[0], x[1] ± i*x[2], return 1
- inline unsigned int solveP3(double* x, double a, double b, double c) {
- double a2 = a * a;
- double q = (a2 - 3 * b) / 9;
- double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
- double r2 = r * r;
- double q3 = q * q * q;
- double A, B;
- if (r2 < q3) {
- double t = r / sqrt(q3);
- if (t < -1) t = -1;
- if (t > 1) t = 1;
- t = acos(t);
- a /= 3;
- q = -2 * sqrt(q);
- x[0] = q * cos(t / 3) - a;
- x[1] = q * cos((t + M_2PI) / 3) - a;
- x[2] = q * cos((t - M_2PI) / 3) - a;
- return 3;
- } else {
- A = -pow(fabs(r) + sqrt(r2 - q3), 1. / 3);
- if (r < 0) A = -A;
- B = (0 == A ? 0 : q / A);
-
- a /= 3;
- x[0] = (A + B) - a;
- x[1] = -0.5 * (A + B) - a;
- x[2] = 0.5 * sqrt(3.) * (A - B);
- if (fabs(x[2]) < eps) {
- x[2] = x[1];
- return 2;
+namespace corsika {
+
+ namespace andre {
+
+ inline std::vector<double> solve_quartic_real(long double a, long double b,
+ long double c, long double d,
+ long double e, double const epsilon) {
+
+ if (std::abs(a) < epsilon) { return solve_cubic_real(b, c, d, e, epsilon); }
+
+ b /= a;
+ c /= a;
+ d /= a;
+ e /= a;
+
+ long double a3 = -c;
+ long double b3 = b * d - 4. * e;
+ long double c3 = -b * b * e - d * d + 4. * c * e;
+
+ // cubic resolvent
+ // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+
+ std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
+ long double y = x3[0]; // there is always at least one solution
+ // The essence - choosing Y with maximal absolute value.
+ if (x3.size() == 3) {
+ if (fabs(x3[1]) > fabs(y)) y = x3[1];
+ if (fabs(x3[2]) > fabs(y)) y = x3[2];
+ }
+
+ long double q1, q2, p1, p2;
+ // h1+h2 = y && h1*h2 = e <=> h^2 -y*h + e = 0 (h === q)
+
+ long double Det = y * y - 4 * e;
+ CORSIKA_LOG_TRACE("Det={}", Det);
+ if (fabs(Det) < epsilon) // in other words - D==0
+ {
+ q1 = q2 = y * 0.5;
+ // g1+g2 = b && g1+g2 = c-y <=> g^2 - b*g + c-y = 0 (p === g)
+ Det = b * b - 4 * (c - y);
+ if (fabs(Det) < epsilon) { // in other words - D==0
+ p1 = p2 = b * 0.5;
+ } else {
+ if (Det < 0) return {};
+ long double sqDet = sqrt(Det);
+ p1 = (b + sqDet) * 0.5;
+ p2 = (b - sqDet) * 0.5;
+ }
+ } else {
+ if (Det < 0) return {};
+ long double sqDet1 = sqrt(Det);
+ q1 = (y + sqDet1) * 0.5;
+ q2 = (y - sqDet1) * 0.5;
+ // g1+g2 = b && g1*h2 + g2*h1 = c ( && g === p ) Krammer
+ p1 = (b * q1 - d) / (q1 - q2);
+ p2 = (d - b * q2) / (q1 - q2);
}
- return 1;
+ // solving quadratic eqs. x^2 + p1*x + q1 = 0
+ // x^2 + p2*x + q2 = 0
+
+ std::vector<double> quad1 = solve_quadratic_real(1, p1, q1);
+ std::vector<double> quad2 = solve_quadratic_real(1, p2, q2);
+ if (quad2.size() > 0) {
+ for (auto val : quad2) quad1.push_back(val);
+ }
+ return quad1;
}
- }
+ } // namespace andre
- //---------------------------------------------------------------------------
- // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
- // Attention - this function returns dynamically allocated array. It has to be released
- // afterwards.
- inline DComplex* solve_quartic(double a, double b, double c, double d) {
- double a3 = -b;
- double b3 = a * c - 4. * d;
- double c3 = -a * a * d - c * c + 4. * b * d;
+ inline std::vector<double> solve_quartic_depressed_real(long double p, long double q,
+ long double r,
+ double const epsilon) {
- // cubic resolvent
- // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+ CORSIKA_LOG_TRACE("quartic-depressed: p={:f}, q={:f}, r={:f}, epsilon={}", p, q, r,
+ epsilon);
- double x3[3];
- unsigned int iZeroes = solveP3(x3, a3, b3, c3);
+ long double const p2 = static_pow<2>(p);
+ long double const q2 = static_pow<2>(q);
- double q1, q2, p1, p2, D, sqD, y;
+ std::vector<double> const resolve_cubic =
+ solve_cubic_real(1, p, p2 / 4 - r, -q2 / 8, epsilon);
- y = x3[0];
- // The essence - choosing Y with maximal absolute value.
- if (iZeroes != 1) {
- if (fabs(x3[1]) > fabs(y)) y = x3[1];
- if (fabs(x3[2]) > fabs(y)) y = x3[2];
- }
+ CORSIKA_LOG_TRACE("resolve_cubic: N={}, m=[{}]", resolve_cubic.size(),
+ fmt::join(resolve_cubic, ", "));
- // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q)
-
- D = y * y - 4 * d;
- if (fabs(D) < eps) // in other words - D==0
- {
- q1 = q2 = y * 0.5;
- // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g)
- D = a * a - 4 * (b - y);
- if (fabs(D) < eps) // in other words - D==0
- p1 = p2 = a * 0.5;
-
- else {
- sqD = sqrt(D);
- p1 = (a + sqD) * 0.5;
- p2 = (a - sqD) * 0.5;
- }
- } else {
- sqD = sqrt(D);
- q1 = (y + sqD) * 0.5;
- q2 = (y - sqD) * 0.5;
- // g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer
- p1 = (a * q1 - c) / (q1 - q2);
- p2 = (c - a * q2) / (q1 - q2);
+ if (!resolve_cubic.size()) return {};
+
+ long double m = 0;
+ for (auto const& v : resolve_cubic) {
+ CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
+ v * (p2 / 4 - r) - q2 / 8));
+ if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
}
+ CORSIKA_LOG_TRACE("check m={}", m);
+ if (m == 0) { return {0}; }
+
+ CORSIKA_LOG_TRACE("check m={}", m);
+
+ long double const quad_term1 = p / 2 + m;
+ long double const quad_term2 = std::sqrt(2 * m);
+ long double const quad_term3 = q / (2 * quad_term2);
- DComplex* retval = new DComplex[4];
-
- // solving quadratic eq. - x^2 + p1*x + q1 = 0
- D = p1 * p1 - 4 * q1;
- if (D < 0.0) {
- retval[0].real(-p1 * 0.5);
- retval[0].imag(sqrt(-D) * 0.5);
- retval[1] = std::conj(retval[0]);
- } else {
- sqD = sqrt(D);
- retval[0].real((-p1 + sqD) * 0.5);
- retval[1].real((-p1 - sqD) * 0.5);
+ std::vector<double> z_quad1 =
+ solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, epsilon);
+ std::vector<double> z_quad2 =
+ solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, epsilon);
+ for (auto const& z : z_quad2) z_quad1.push_back(z);
+ return z_quad1;
+ }
+
+ inline std::vector<double> solve_quartic_real(long double a, long double b,
+ long double c, long double d,
+ long double e, double const epsilon) {
+
+ CORSIKA_LOG_TRACE("quartic: a={:f}, b={:f}, c={:f}, d={:f}, e={:f}, epsilon={}", a, b,
+ c, d, e, epsilon);
+
+ if (std::abs(a) < epsilon) { // this is just a quadratic
+ return solve_cubic_real(b, c, d, e, epsilon);
}
- // solving quadratic eq. - x^2 + p2*x + q2 = 0
- D = p2 * p2 - 4 * q2;
- if (D < 0.0) {
- retval[2].real(-p2 * 0.5);
- retval[2].imag(sqrt(-D) * 0.5);
- retval[3] = std::conj(retval[2]);
- } else {
- sqD = sqrt(D);
- retval[2].real((-p2 + sqD) * 0.5);
- retval[3].real((-p2 - sqD) * 0.5);
+ if ((std::abs(a - 1) < epsilon) &&
+ (std::abs(b) < epsilon)) { // this is a depressed quartic
+ return solve_quartic_depressed_real(c, d, e, epsilon);
}
- return retval;
+ long double const b2 = static_pow<2>(b);
+ long double const b3 = static_pow<3>(b);
+ long double const b4 = static_pow<4>(b);
+ long double const a2 = static_pow<2>(a);
+ long double const a3 = static_pow<3>(a);
+ long double const a4 = static_pow<4>(a);
+
+ long double const p = (c * a * 8 - b2 * 3) / (a4 * 8);
+ long double const q = (b3 - b * c * a * 4 + d * a2 * 8) / (a4 * 8);
+ long double const r =
+ (-b4 * 3 + e * a3 * 256 - b * d * a2 * 64 + b2 * c * a * 16) / (a4 * 256);
+
+ std::vector<double> zs = solve_quartic_depressed_real(p, q, r, epsilon);
+ CORSIKA_LOG_TRACE("quartic: solve_depressed={}, b/4a={}", fmt::join(zs, ", "),
+ b / (4 * a));
+ for (auto& z : zs) { z -= b / (4 * a); }
+ CORSIKA_LOG_TRACE("quartic: solve_quartic_real returns={}", fmt::join(zs, ", "));
+ return zs;
}
-
-} // namespace corsika::quartic_solver
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/detail/modules/LongitudinalProfile.inl b/corsika/detail/modules/LongitudinalProfile.inl
index 48deafe6e79a474f06ecf664d9b70640b8e7756b..5b1f8438bc0438c45940e34c9502cd31ec637196 100644
--- a/corsika/detail/modules/LongitudinalProfile.inl
+++ b/corsika/detail/modules/LongitudinalProfile.inl
@@ -35,15 +35,15 @@ namespace corsika {
GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
- CORSIKA_LOG_TRACE("pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2",
+ CORSIKA_LOG_DEBUG("longprof: pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2",
vTrack.getPosition(0).getCoordinates() / 1_m,
vTrack.getPosition(1).getCoordinates() / 1_m,
grammageStart / 1_g * square(1_cm),
grammageEnd / 1_g * square(1_cm));
// Note: particle may go also "upward", thus, grammageEnd<grammageStart
- const int binStart = std::ceil(grammageStart / dX_);
- const int binEnd = std::floor(grammageEnd / dX_);
+ int const binStart = std::ceil(grammageStart / dX_);
+ int const binEnd = std::floor(grammageEnd / dX_);
for (int b = binStart; b <= binEnd; ++b) {
if (pid == Code::Gamma) {
@@ -66,6 +66,7 @@ namespace corsika {
inline void LongitudinalProfile::save(std::string const& filename, const int width,
const int precision) {
+ CORSIKA_LOG_DEBUG("Write longprof to {}", filename);
std::ofstream f{filename};
f << "# X / g·cm¯², gamma, e+, e-, mu+, mu-, all hadrons" << std::endl;
for (size_t b = 0; b < profiles_.size(); ++b) {
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index c4da50e940f88e3b8a3d3e9f36ac0018eeb09466..e5be30bedf2f81b65d68faeca314669c58c31ac6 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -24,6 +24,7 @@ namespace corsika {
, count_ground_(0)
, xAxis_(x_axis.normalized())
, yAxis_(obsPlane.getNormal().cross(xAxis_)) {
+ CORSIKA_LOG_DEBUG("Plane height: {}", obsPlane.getCenter());
outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m"
<< std::endl;
}
@@ -35,7 +36,18 @@ namespace corsika {
/*
The current step did not yet reach the ObservationPlane, do nothing now and wait:
*/
- if (!stepLimit) { return ProcessReturn::Ok; }
+ if (!stepLimit) {
+#ifdef DEBUG
+ if (deleteOnHit_) {
+ LengthType const check =
+ (particle.getPosition() - plane_.getCenter()).dot(plane_.getNormal());
+ if (check < 0_m) {
+ CORSIKA_LOG_DEBUG("PARTICLE AVOIDED OBSERVATIONPLANE {}", check);
+ }
+ }
+#endif
+ return ProcessReturn::Ok;
+ }
HEPEnergyType const energy = particle.getEnergy();
Point const pointOfIntersection = particle.getPosition();
@@ -45,6 +57,8 @@ namespace corsika {
<< ' ' << displacement.dot(xAxis_) / 1_m << ' '
<< displacement.dot(yAxis_) / 1_m << '\n';
+ CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
+
if (deleteOnHit_) {
count_ground_++;
energy_ground_ += energy;
@@ -58,13 +72,14 @@ namespace corsika {
corsika::setup::Stack::particle_type const& particle,
corsika::setup::Trajectory const& trajectory) {
+ CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}", particle.asString(),
+ particle.getPosition(), particle.getDirection(), plane_.asString());
+
Intersections const intersection =
setup::Tracking::intersect<corsika::setup::Stack::particle_type>(particle,
plane_);
TimeType const timeOfIntersection = intersection.getEntry();
- CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}, timeOfIntersection={}",
- particle.asString(), particle.getPosition(),
- particle.getDirection(), plane_.asString(), timeOfIntersection);
+ CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection);
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
diff --git a/corsika/detail/modules/StackInspector.inl b/corsika/detail/modules/StackInspector.inl
index caae6c59c8cb43b01cd1e9789366b51db2f9ec14..2d9c8238a4a5bb74b6ae80d6b218df0f1ef31221 100644
--- a/corsika/detail/modules/StackInspector.inl
+++ b/corsika/detail/modules/StackInspector.inl
@@ -72,13 +72,13 @@ namespace corsika {
CORSIKA_LOG_INFO(
"StackInspector: "
- " time= {}"
- ", running= {} seconds"
+ " time={}"
+ ", running={} seconds"
" ( {}%)"
- ", nStep= {}"
- ", stackSize= {}"
- ", Estack= {} GeV"
- ", ETA=",
+ ", nStep={}"
+ ", stackSize={}"
+ ", Estack={} GeV"
+ ", ETA={}",
std::put_time(std::localtime(&now_time), "%T"), elapsed_seconds.count(),
int(progress * 100), getStep(), vS.getSize(), Etot / 1_GeV,
std::put_time(std::localtime(&eta_time), "%T"));
diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl
index ebd52339e72eddde23f1f50dc66e709712165963..e55d40cb9cfb64d23aaccd21fc11347e10035226 100644
--- a/corsika/detail/modules/TrackWriter.inl
+++ b/corsika/detail/modules/TrackWriter.inl
@@ -30,9 +30,14 @@ namespace corsika {
<< '\n';
}
+ inline TrackWriter::~TrackWriter() { file_.close(); }
+
template <typename TParticle, typename TTrack>
inline ProcessReturn TrackWriter::doContinuous(TParticle const& vP, TTrack const& vT,
bool const) {
+
+ CORSIKA_LOG_DEBUG("TrackWriter");
+
auto const start = vT.getPosition(0).getCoordinates();
auto const delta = vT.getPosition(1).getCoordinates() - start;
auto const pdg = static_cast<int>(get_PDG(vP.getPID()));
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 73fbf351707c35c894dd8cf48e0d250589909527..310613d24d2d8fcbd75402822773891b8374dfca 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -184,9 +184,9 @@ namespace corsika {
energy * 0.9, // either 10% relative loss max., or
get_energy_threshold(vParticle.getPID()) // energy thresholds globally defined for
// individual particles
- *
- 0.99 // need to go 1% below global e-cut to assure removal in ParticleCut. The
- // 1% does not matter since at cut-time the entire energy is removed.
+ * 0.99999 // need to go slightly below global e-cut to assure removal in
+ // ParticleCut. This does not matter since at cut-time the entire
+ // energy is removed.
);
auto const maxGrammage = (vParticle.getEnergy() - energy_lim) / dEdX;
diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl
index 217fe30b72c3fa27d1ad3794a4a6b1ec8852ffd5..9a86aa880e2333d9247bd2e18d347a83d9184541 100644
--- a/corsika/detail/modules/qgsjetII/Interaction.inl
+++ b/corsika/detail/modules/qgsjetII/Interaction.inl
@@ -27,20 +27,14 @@
namespace corsika::qgsjetII {
- inline Interaction::Interaction(const std::string& dataPath)
- : data_path_(dataPath) {
- if (dataPath == "") {
- if (std::getenv("CORSIKA_DATA")) {
- data_path_ = std::string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/";
- CORSIKA_LOG_DEBUG("Searching for QGSJetII data tables in {}", data_path_);
- }
- }
+ inline Interaction::Interaction(boost::filesystem::path dataPath) {
+ CORSIKA_LOG_DEBUG("Reading QGSJetII data tables from {}", dataPath);
// initialize QgsjetII
static bool initialized = false;
if (!initialized) {
qgset_();
- datadir DIR(data_path_);
+ datadir DIR(dataPath.string() + "/");
qgaini_(DIR.data);
initialized = true;
}
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 5770fe439dbc309588e990da8669c9f85ceb619c..a6553a71dd3646534645bd8a1ee90c9fa29b08ca 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -34,15 +34,14 @@ namespace corsika {
return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
double(0));
} // charge of the particle
- int const chargeNumber = particle.getChargeNumber();
+ ElectricChargeType const charge = particle.getCharge();
auto const* currentLogicalVolumeNode = particle.getNode();
MagneticFieldVector const& magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
- VelocityVector velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
+ VelocityVector velocity = particle.getVelocity();
decltype(meter / (second * volt)) k =
- chargeNumber * constants::cSquared * 1_eV /
+ charge * constants::cSquared * 1_eV /
(velocity.getNorm() * particle.getEnergy() * 1_V);
DirectionVector direction = velocity.normalized();
auto position = particle.getPosition(); // First Movement
@@ -58,23 +57,22 @@ namespace corsika {
template <typename TParticle>
inline auto Tracking::getTrack(TParticle const& particle) {
- VelocityVector const initialVelocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
+ VelocityVector const initialVelocity = particle.getVelocity();
- auto const position = particle.getPosition();
+ auto const& position = particle.getPosition();
CORSIKA_LOG_DEBUG(
- "Tracking pid: {}"
- " , E = {} GeV",
- particle.getPID(), particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking pos: {}", position.getCoordinates());
- CORSIKA_LOG_DEBUG("Tracking E: {} GeV", particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking p: {} GeV",
- particle.getMomentum().getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking v: {} ", initialVelocity.getComponents());
+ "TrackingLeapfrog_Curved pid: {}"
+ " , E = {} GeV \n"
+ "\tTracking pos: {} \n"
+ "\tTracking p: {} GeV \n"
+ "\tTracking v: {}",
+ particle.getPID(), particle.getEnergy() / 1_GeV, position.getCoordinates(),
+ particle.getMomentum().getComponents() / 1_GeV,
+ initialVelocity.getComponents());
typedef
typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
- node_type& volumeNode = *particle.getNode();
+ node_type const& volumeNode = *particle.getNode();
// for the event of magnetic fields and curved trajectories, we need to limit
// maximum step-length since we need to follow curved
@@ -83,25 +81,28 @@ namespace corsika {
MagneticFieldVector const& magneticfield =
volumeNode.getModelProperties().getMagneticField(position);
MagneticFluxType const magnitudeB = magneticfield.getNorm();
- int const chargeNumber = particle.getChargeNumber();
- bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
+ ElectricChargeType const charge = particle.getCharge();
+ bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T;
- if (no_deflection) { return getLinearTrajectory(particle); }
+ if (no_deflection) {
+ CORSIKA_LOG_TRACE("no_deflection");
+ return getLinearTrajectory(particle);
+ }
- HEPMomentumType const pAlongB_delta =
+ HEPMomentumType const p_perp =
(particle.getMomentum() -
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
- if (pAlongB_delta == 0_GeV) {
+ if (p_perp < 1_eV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
- CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
+ CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
return getLinearTrajectory(particle);
}
- LengthType const gyroradius =
- (pAlongB_delta * 1_V / (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
+ constants::c / (abs(charge) * magnitudeB));
double const maxRadians = 0.01;
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
@@ -113,8 +114,15 @@ namespace corsika {
// intersection
auto [minTime, minNode] = nextIntersect(particle, steplimit_time);
- auto const k =
- chargeNumber * constants::cSquared * 1_eV / (particle.getEnergy() * 1_V);
+ auto const p_norm =
+ constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+ particle.getMomentum().getNorm()); // kg *m /s
+ // k = q/|p|
+ decltype(1 / (tesla * second)) const k =
+ charge / p_norm *
+ initialVelocity.getNorm(); // since we use steps in time and not length
+ // units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1/(T s)
+
return std::make_tuple(
LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
minTime), // trajectory
@@ -129,99 +137,104 @@ namespace corsika {
return Intersections();
}
- int const chargeNumber = particle.getChargeNumber();
+ ElectricChargeType const charge = particle.getCharge();
auto const& position = particle.getPosition();
auto const* currentLogicalVolumeNode = particle.getNode();
MagneticFieldVector const& magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
- VelocityVector const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- DirectionVector const directionBefore =
- velocity.normalized(); // determine steplength to next volume
+ VelocityVector const velocity = particle.getVelocity();
+ DirectionVector const directionBefore = velocity.normalized();
auto const projectedDirection = directionBefore.cross(magneticfield);
auto const projectedDirectionSqrNorm = projectedDirection.getSquaredNorm();
bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
- if (chargeNumber == 0 || magneticfield.getNorm() == 0_T || isParallel) {
+ if ((charge == 0 * constants::e) || magneticfield.getNorm() == 0_T || isParallel) {
return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
}
bool const numericallyInside = sphere.contains(particle.getPosition());
CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
- auto const absVelocity = velocity.getNorm();
- auto const energy = particle.getEnergy();
+ SpeedType const absVelocity = velocity.getNorm();
+ auto const p_norm =
+ constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+ particle.getMomentum().getNorm()); // km * m /s
// this is: k = q/|p|
- auto const k =
- chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V);
+ auto const k = charge / p_norm;
- auto const direction_B_perp = directionBefore.cross(magneticfield);
- auto const denom = 4. / (direction_B_perp.getSquaredNorm() * k * k);
+ MagneticFieldVector const direction_x_B = directionBefore.cross(magneticfield);
+ auto const denom = 4. / (direction_x_B.getSquaredNorm() * k * k);
Vector<length_d> const deltaPos = position - sphere.getCenter();
- double const b = (direction_B_perp.dot(deltaPos) * k + 1) * denom / (1_m * 1_m);
+ double const b = (direction_x_B.dot(deltaPos) * k + 1) * denom / (1_m * 1_m);
double const c = directionBefore.dot(deltaPos) * 2 * denom / (1_m * 1_m * 1_m);
LengthType const deltaPosLength = deltaPos.getNorm();
double const d = (deltaPosLength + sphere.getRadius()) *
(deltaPosLength - sphere.getRadius()) * denom /
(1_m * 1_m * 1_m * 1_m);
CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
- std::complex<double> const* solutions = quartic_solver::solve_quartic(0, b, c, d);
+ std::vector<double> solutions = andre::solve_quartic_real(1, 0, b, c, d);
LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0;
- for (int i = 0; i < 4; i++) {
- if (solutions[i].imag() == 0) {
- LengthType const dist = solutions[i].real() * 1_m;
- CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
- if (numericallyInside) {
- // there must be an entry (negative) and exit (positive) solution
- if (dist < -0.0001_m) { // security margin to assure transfer to next
- // logical volume
- if (first_entry == 0) {
- d_enter = dist;
- } else {
- d_enter = std::max(d_enter, dist); // closest negative to zero (-1e-4) m
- }
- first_entry++;
-
- } else { // thus, dist > -0.0001_m
-
- if (first_exit == 0) {
- d_exit = dist;
- } else {
- d_exit = std::min(d_exit, dist); // closest positive to zero (-1e-4) m
- }
- first_exit++;
+ for (auto solution : solutions) {
+ LengthType const dist = solution * 1_m;
+ CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
+ if (numericallyInside) {
+ // there must be an entry (negative) and exit (positive) solution
+ if (dist < 0.0001_m) { // security margin to assure
+ // transfer to next logical volume
+ // (even if dist suggest marginal
+ // entry already, which we
+ // classify as numerical artifact)
+ if (first_entry == 0) {
+ d_enter = dist;
+ } else {
+ d_enter = std::max(d_enter, dist); // closest negative to zero >1e-4 m
}
- first = int(first_exit > 0) + int(first_entry > 0);
+ first_entry++;
- } else { // thus, numericallyInside == false
+ } else { // thus, dist > +0.0001_m
- // both physical solutions (entry, exit) must be positive, and as small as
- // possible
- if (dist < -0.0001_m) { // need small numerical margin, to assure transport
- // into next logical volume
- continue;
+ if (first_exit == 0) {
+ d_exit = dist;
+ } else {
+ d_exit = std::min(d_exit, dist); // closest positive to zero >1e-4 m
}
- if (first == 0) {
+ first_exit++;
+ }
+ first = int(first_exit > 0) + int(first_entry > 0);
+
+ } else { // thus, numericallyInside == false
+
+ // both physical solutions (entry, exit) must be positive, and as small as
+ // possible
+ if (dist < -0.0001_m) { // need small numerical margin, to
+ // assure transport. We consider
+ // begin marginally already inside
+ // next volume (besides
+ // numericallyInside=false) as numerical glitch.
+ // into next logical volume
+ continue;
+ }
+ if (first == 0) {
+ d_enter = dist;
+ } else {
+ if (dist < d_enter) {
+ d_exit = d_enter;
d_enter = dist;
} else {
- if (dist < d_enter) {
- d_exit = d_enter;
- d_enter = dist;
- } else {
- d_exit = dist;
- }
+ d_exit = dist;
}
- first++;
}
- } // loop over solutions
- }
- delete[] solutions;
+ first++;
+ }
+ } // loop over solutions
- if (first != 2) { // entry and exit points found
- CORSIKA_LOG_DEBUG("no intersection! count={}", first);
+ if (first == 0) { // entry and exit points found
+ CORSIKA_LOG_DEBUG(
+ "no intersections found: count={}, first_entry={}, first_exit={}", first,
+ first_entry, first_exit);
return Intersections();
}
return Intersections(d_enter / absVelocity, d_exit / absVelocity);
@@ -231,79 +244,82 @@ namespace corsika {
inline Intersections Tracking::intersect(TParticle const& particle,
Plane const& plane) {
- int chargeNumber;
- if (is_nucleus(particle.getPID())) {
- chargeNumber = particle.getNuclearZ();
- } else {
- chargeNumber = get_charge_number(particle.getPID());
- }
- auto const* currentLogicalVolumeNode = particle.getNode();
- VelocityVector const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- auto const absVelocity = velocity.getNorm();
- DirectionVector const direction =
- velocity.normalized(); // determine steplength to next volume
- Point const position = particle.getPosition();
-
- auto const magneticfield =
- currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
+ CORSIKA_LOG_TRACE("intersection particle with plane");
- if (chargeNumber != 0 && abs(plane.getNormal().dot(velocity.cross(magneticfield))) >
- 1e-6_T * 1_m / 1_s) {
+ ElectricChargeType const charge = particle.getCharge();
+
+ if (charge != 0 * constants::e) {
auto const* currentLogicalVolumeNode = particle.getNode();
+ VelocityVector const velocity = particle.getVelocity();
+ auto const absVelocity = velocity.getNorm();
+ DirectionVector const direction = velocity.normalized();
+ Point const& position = particle.getPosition();
+
auto const magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
- auto const k =
- chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm();
- auto const direction_B_perp = direction.cross(magneticfield);
- auto const denom = plane.getNormal().dot(direction_B_perp) * k;
- auto const sqrtArg =
- direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
- (plane.getNormal().dot(position - plane.getCenter()) * denom * 2);
+ // solve: denom x^2 + p x + q =0 for x = delta-l
+
+ auto const direction_x_B = direction.cross(magneticfield);
+ double const denom = charge *
+ plane.getNormal().dot(direction_x_B) // unit: C*T = kg/s
+ / 1_kg * 1_s;
+
+ CORSIKA_LOG_TRACE("denom={}", denom);
+
+ auto const p_norm =
+ constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+ particle.getMomentum().getNorm()); // unit: kg * m/s
+
+ double const p = (2 * p_norm * direction.dot(plane.getNormal())) // unit: kg*m/s
+ / (1_m * 1_kg) * 1_s;
+ double const q =
+ (2 * p_norm *
+ plane.getNormal().dot(position - plane.getCenter())) // unit: kg*m/s *m
+ / (1_m * 1_m * 1_kg) * 1_s;
+
+ std::vector<double> const deltaLs = solve_quadratic_real(denom, p, q);
- if (sqrtArg < 0) {
+ CORSIKA_LOG_TRACE("deltaLs=[{}]", fmt::join(deltaLs, ", "));
+
+ if (deltaLs.size() == 0) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
}
- double const sqrSqrtArg = sqrt(sqrtArg);
- auto const norm_projected =
- direction.dot(plane.getNormal()) / direction.getNorm();
- LengthType const MaxStepLength1 = (sqrSqrtArg - norm_projected) / denom;
- LengthType const MaxStepLength2 = (-sqrSqrtArg - norm_projected) / denom;
- CORSIKA_LOG_TRACE("MaxStepLength1={}, MaxStepLength2={}", MaxStepLength1,
- MaxStepLength2);
+ // select smallest but positive solution
+ bool first = true;
+ LengthType maxStepLength = 0_m;
+ for (auto const& deltaL : deltaLs) {
+ if (deltaL < 0) continue;
+ if (first) {
+ first = false;
+ maxStepLength = deltaL * meter;
+ } else if (maxStepLength > deltaL * meter) {
+ maxStepLength = deltaL * meter;
+ }
+ }
- // check: both intersections in past
- if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
+ // check: both intersections in past, or no valid intersection
+ if (first) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
-
- // check: next intersection is MaxStepLength2
- } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
- CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength2);
- return Intersections(
- MaxStepLength2 *
- (direction + direction_B_perp * MaxStepLength2 * k / 2).getNorm() /
- absVelocity);
-
- // check: next intersections is MaxStepLength1
- } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
- CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength1);
- return Intersections(
- MaxStepLength1 *
- (direction + direction_B_perp * MaxStepLength1 * k / 2).getNorm() /
- absVelocity);
}
- CORSIKA_LOG_WARN(
- "Particle wasn't tracked with curved trajectory -> straight (is this an "
- "ERROR?)");
+ CORSIKA_LOG_TRACE("maxStepLength={}", maxStepLength);
+
+ // with final length correction
+ auto const corr =
+ (direction + direction_x_B * maxStepLength * charge / (p_norm * 2))
+ .getNorm() /
+ absVelocity; // unit: s/m
+
+ return Intersections(maxStepLength * corr); // unit: s
- } // end if curved-tracking
+ } // no charge
- CORSIKA_LOG_TRACE("straight tracking with chargeNumber={}, B={}", chargeNumber,
- magneticfield);
+ CORSIKA_LOG_TRACE("(plane) straight tracking with charge={}, B={}", charge,
+ particle.getNode()->getModelProperties().getMagneticField(
+ particle.getPosition()));
return tracking_line::Tracking::intersect(particle, plane);
}
@@ -330,7 +346,7 @@ namespace corsika {
straightTrajectory.getLine().getVelocity(),
MagneticFieldVector(particle.getPosition().getCoordinateSystem(), 0_T, 0_T,
0_T),
- square(0_m) / (square(1_s) * 1_V),
+ 0 * square(meter) / (square(second) * volt),
straightTrajectory.getDuration()), // trajectory
minNode); // next volume node
}
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
index 61ad3793d0537c47c0efb1b870d11281c663072c..4e7efed775468a6525a236fe7b5a72b45896e518 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
@@ -29,15 +29,17 @@ namespace corsika {
particle.getMomentum() / particle.getEnergy() * constants::c;
Point const initialPosition = particle.getPosition();
+
CORSIKA_LOG_DEBUG(
- "TrackingB pid: {}"
- " , E = {} GeV",
- particle.getPID(), particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB pos: {}", initialPosition.getCoordinates());
- CORSIKA_LOG_DEBUG("TrackingB E: {} GeV", particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB p: {} GeV",
- particle.getMomentum().getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB v: {} ", initialVelocity.getComponents());
+ "TrackingLeapFrogStraight pid: {}"
+ " , E = {} GeV \n"
+ "\tTracking pos: {} \n"
+ "\tTracking p: {} GeV \n"
+ "\tTracking v: {} ",
+ particle.getPID(), particle.getEnergy() / 1_GeV,
+ initialPosition.getCoordinates(),
+ particle.getMomentum().getComponents() / 1_GeV,
+ initialVelocity.getComponents());
typedef decltype(particle.getNode()) node_type;
node_type const volumeNode = particle.getNode();
@@ -50,14 +52,13 @@ namespace corsika {
}
// charge of the particle, and magnetic field
- const int chargeNumber = particle.getChargeNumber();
+ auto const charge = particle.getCharge();
auto magneticfield =
volumeNode->getModelProperties().getMagneticField(initialPosition);
auto const magnitudeB = magneticfield.getNorm();
- CORSIKA_LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
- magneticfield.getComponents() / 1_uT, chargeNumber,
- magnitudeB / 1_T);
- bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
+ CORSIKA_LOG_DEBUG("field={} uT, charge={}, magnitudeB={} uT",
+ magneticfield.getComponents() / 1_uT, charge, magnitudeB / 1_T);
+ bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T;
// check, where the first halve-step direction has geometric intersections
auto const [initialTrack, initialTrackNextVolume] =
@@ -75,26 +76,26 @@ namespace corsika {
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
- HEPMomentumType const pAlongB_delta =
+ HEPMomentumType const p_perp =
(particle.getMomentum() -
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
- if (pAlongB_delta == 0_GeV) {
+ if (p_perp == 0_GeV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
- CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
+ CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
- LengthType const gyroradius =
- (pAlongB_delta * 1_V / (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
+ constants::c / (abs(charge) * magnitudeB));
// we need to limit maximum step-length since we absolutely
// need to follow strongly curved trajectories segment-wise,
// at least if we don't employ concepts as "Helix
// Trajectories" or similar
- double const maxRadians = 0.01;
+ double const maxRadians = 0.001;
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
@@ -111,9 +112,9 @@ namespace corsika {
firstHalveSteplength, steplimit, initialTrackLength);
// perform the first halve-step
Point const position_mid = initialPosition + direction * firstHalveSteplength;
- auto const k =
- chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm();
- auto const new_direction =
+ auto const k = charge / (constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+ particle.getMomentum().getNorm()));
+ DirectionVector const new_direction =
direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
auto const new_direction_norm = new_direction.getNorm(); // by design this is >1
CORSIKA_LOG_DEBUG(
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 9a4659250286414cc88f837e96c3ac5898952006..cad8ac39d99ebf67af35c19eea7974897806b3f9 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -20,6 +20,7 @@
#include <type_traits>
#include <utility>
+#include <cmath>
namespace corsika::tracking_line {
@@ -28,16 +29,15 @@ namespace corsika::tracking_line {
VelocityVector const initialVelocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
- auto const initialPosition = particle.getPosition();
+ auto const& initialPosition = particle.getPosition();
CORSIKA_LOG_DEBUG(
- "Tracking pid: {}"
- " , E = {} GeV",
- particle.getPID(), particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking pos: {}", initialPosition.getCoordinates());
- CORSIKA_LOG_DEBUG("Tracking E: {} GeV", particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking p: {} GeV",
- particle.getMomentum().getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking v: {} ", initialVelocity.getComponents());
+ "TrackingStraight pid: {}"
+ " , E = {} GeV \n"
+ "\tTracking pos: {} \n"
+ "\tTracking p: {} GeV \n"
+ "\tTracking v: {}",
+ particle.getPID(), particle.getEnergy() / 1_GeV, initialPosition.getCoordinates(),
+ particle.getMomentum().getComponents() / 1_GeV, initialVelocity.getComponents());
// traverse the environment volume tree and find next
// intersection
@@ -51,7 +51,8 @@ namespace corsika::tracking_line {
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Sphere const& sphere) {
- auto const delta = particle.getPosition() - sphere.getCenter();
+ auto const& position = particle.getPosition();
+ auto const delta = position - sphere.getCenter();
auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
auto const vSqNorm = velocity.getSquaredNorm();
auto const R = sphere.getRadius();
@@ -63,6 +64,8 @@ namespace corsika::tracking_line {
if (discriminant.magnitude() > 0) {
auto const sqDisc = sqrt(discriminant);
auto const invDenom = 1 / vSqNorm;
+
+ CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position));
return Intersections((-vDotDelta - sqDisc) * invDenom,
(-vDotDelta + sqDisc) * invDenom);
}
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index ee29642e00bbd408c355a5c2c645abef312564ec..f0409e815995f93aedd0175cdd31be5b2a834ae4 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -29,7 +29,7 @@
namespace corsika::urqmd {
- inline UrQMD::UrQMD(boost::filesystem::path const& xs_file) {
+ inline UrQMD::UrQMD(boost::filesystem::path xs_file) {
readXSFile(xs_file);
::urqmd::iniurqmdc8_();
}
@@ -268,10 +268,8 @@ namespace corsika::urqmd {
2 * ::urqmd::options_.CTParam[30 - 1];
::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV);
- if (projectileCode == Code::K0Long) {
+ if (projectileCode == Code::K0Long || projectileCode == Code::K0Short) {
projectileCode = booleanDist_(RNG_) ? Code::K0 : Code::K0Bar;
- } else if (projectileCode == Code::K0Short) {
- throw std::runtime_error("K0Short should not interact");
}
auto const [ityp, iso3] = convertToUrQMD(projectileCode);
@@ -396,7 +394,7 @@ namespace corsika::urqmd {
return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
}
- inline void UrQMD::readXSFile(boost::filesystem::path const& filename) {
+ inline void UrQMD::readXSFile(boost::filesystem::path const filename) {
boost::filesystem::ifstream file(filename, std::ios::in);
if (!file.is_open()) {
@@ -431,6 +429,7 @@ namespace corsika::urqmd {
std::getline(file, line);
}
}
+ file.close();
}
} // namespace corsika::urqmd
diff --git a/corsika/detail/stack/NuclearStackExtension.inl b/corsika/detail/stack/NuclearStackExtension.inl
index 728647ca0bfc5ea89738946ccdef9c36cda31ce2..5dab28ab2506181580fd510e28164c507b70b878 100644
--- a/corsika/detail/stack/NuclearStackExtension.inl
+++ b/corsika/detail/stack/NuclearStackExtension.inl
@@ -113,6 +113,14 @@ namespace corsika::nuclear_stack {
return super_type::getMass();
}
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline ElectricChargeType NuclearParticleInterface<
+ InnerParticleInterface, StackIteratorInterface>::getCharge() const {
+ if (super_type::getPID() == Code::Nucleus) return getNuclearZ() * constants::e;
+ return super_type::getCharge();
+ }
+
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline int16_t NuclearParticleInterface<
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index 7ee12186ace6f189b2c375b8e1a04e502a3cf06f..b9ee5361bdd00301a1ffe2f2c4ecd75563379ec2 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -130,8 +130,9 @@ namespace corsika {
*/
void step(Particle& vParticle);
- ProcessReturn decay(TStackView& view);
- ProcessReturn interaction(TStackView& view);
+ ProcessReturn decay(TStackView& view, InverseTimeType initial_inv_decay_time);
+ ProcessReturn interaction(TStackView& view,
+ InverseGrammageType initial_inv_int_length);
void setEventType(TStackView& view, history::EventType);
// data members
diff --git a/corsika/framework/core/Logging.hpp b/corsika/framework/core/Logging.hpp
index 67f059c2fe6d1dc5d4b228e551f4ebbbcd1acf9e..d127040b036dab150ce06a40644050f20a3fab2a 100644
--- a/corsika/framework/core/Logging.hpp
+++ b/corsika/framework/core/Logging.hpp
@@ -49,7 +49,8 @@ namespace corsika {
/*
* The default pattern for CORSIKA8 loggers.
*/
- const std::string default_pattern{"[%n:%^%-8l%$] %v"};
+ const std::string minimal_pattern{"[%n:%^%-8l%$] %v"};
+ const std::string default_pattern{"[%n:%^%-8l%$(%s:%#)] %v"};
const std::string source_pattern{"[%n:%^%-8l%$(%s:%!:%#)] %v"};
/**
@@ -91,7 +92,6 @@ namespace corsika {
*/
static inline std::shared_ptr<spdlog::logger> corsika_logger =
get_logger("corsika", true);
- // corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// many of these free functions are special to the logging
// infrastructure so we hide them in the corsika::logging namespace.
diff --git a/corsika/framework/geometry/LeapFrogTrajectory.hpp b/corsika/framework/geometry/LeapFrogTrajectory.hpp
index 6b67bf6345a5c79400d72d9ec16e13afa380e13f..11e3534a19ae6c8ead194077ae515afb194bbd4c 100644
--- a/corsika/framework/geometry/LeapFrogTrajectory.hpp
+++ b/corsika/framework/geometry/LeapFrogTrajectory.hpp
@@ -39,7 +39,7 @@ namespace corsika {
LeapFrogTrajectory(Point const& pos, VelocityVector const& initialVelocity,
MagneticFieldVector const& Bfield,
- decltype(square(meter) / (square(second) * volt)) const k,
+ decltype(1 / (tesla * second)) const k,
TimeType const timeStep) // leap-from total length
: initialPosition_(pos)
, initialVelocity_(initialVelocity)
@@ -74,7 +74,7 @@ namespace corsika {
VelocityVector initialVelocity_;
DirectionVector initialDirection_;
MagneticFieldVector magneticfield_;
- decltype(square(meter) / (square(second) * volt)) k_;
+ decltype(1 / (tesla * second)) k_;
TimeType timeStep_;
};
diff --git a/corsika/framework/process/ContinuousProcess.hpp b/corsika/framework/process/ContinuousProcess.hpp
index c549bc714da8e4c6b655afbb8690568a2f4c4e24..51a7ea7109439401c23e39b722fdd731ff8b1ddb 100644
--- a/corsika/framework/process/ContinuousProcess.hpp
+++ b/corsika/framework/process/ContinuousProcess.hpp
@@ -18,47 +18,46 @@
namespace corsika {
/**
- @ingroup Processes
- @{
-
- Processes with continuous effects along a particle Trajectory
-
- Create a new ContinuousProcess, e.g. for XYModel, via
- @code{.cpp}
- class XYModel : public ContinuousProcess<XYModel> {};
- @endcode
-
- and provide two necessary interface methods:
- @code{.cpp}
- template <typename TParticle, typename TTrack>
- LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const;
- @endcode
-
- which allows any ContinuousProcess to tell to CORSIKA a maximum
- allowed step length. Such step-length limitation, if it turns out
- to be smaller/sooner than any other limit (decay length,
- interaction length, other continuous processes, geometry, etc.)
- will lead to a limited step length.
-
- @code{.cpp}
- template <typename TParticle, typename TTrack>
- ProcessReturn doContinuous(TParticle& p, TTrack const& t, bool const stepLimit)
- const;
- @endcode
-
- which applied any continuous effects on Particle p along
- Trajectory t. The particle in all typical scenarios will be
- altered by a doContinuous. The flag stepLimit will be true if the
- preious evaluation of getMaxStepLength resulted in this
- particular ContinuousProcess to be responsible for the step
- length limit on the current track t. This information can be
- expoited and avoid e.g. any uncessary calculations.
-
- Particle and Track are the valid classes to
- access particles and track (Trajectory) data on the Stack. Those two methods
- do not need to be templated, they could use the types
- e.g. corsika::setup::Stack::particle_type -- but by the cost of
- loosing all flexibility otherwise provided.
+ * @ingroup Processes
+ * @{
+ * Processes with continuous effects along a particle Trajectory.
+ *
+ * Create a new ContinuousProcess, e.g. for XYModel, via:
+ * @code{.cpp}
+ * class XYModel : public ContinuousProcess<XYModel> {};
+ * @endcode
+ *
+ * and provide two necessary interface methods:
+ * @code{.cpp}
+ * template <typename TParticle, typename TTrack>
+ * LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const;
+ * @endcode
+
+ * which allows any ContinuousProcess to tell to CORSIKA a maximum
+ * allowed step length. Such step-length limitation, if it turns out
+ * to be smaller/sooner than any other limit (decay length,
+ * interaction length, other continuous processes, geometry, etc.)
+ * will lead to a limited step length.
+
+ * @code{.cpp}
+ * template <typename TParticle, typename TTrack>
+ * ProcessReturn doContinuous(TParticle& p, TTrack const& t, bool const stepLimit)
+ * const;
+ * @endcode
+
+ * which applied any continuous effects on Particle p along
+ * Trajectory t. The particle in all typical scenarios will be
+ * altered by a doContinuous. The flag stepLimit will be true if the
+ * preious evaluation of getMaxStepLength resulted in this
+ * particular ContinuousProcess to be responsible for the step
+ * length limit on the current track t. This information can be
+ * expoited and avoid e.g. any uncessary calculations.
+
+ * Particle and Track are the valid classes to
+ * access particles and track (Trajectory) data on the Stack. Those two methods
+ * do not need to be templated, they could use the types
+ * e.g. corsika::setup::Stack::particle_type -- but by the cost of
+ * loosing all flexibility otherwise provided.
*/
@@ -68,8 +67,8 @@ namespace corsika {
};
/**
- * ProcessTraits specialization to flag ContinuousProcess objects
- **/
+ * ProcessTraits specialization to flag ContinuousProcess objects.
+ */
template <typename TProcess>
struct is_continuous_process<
TProcess, std::enable_if_t<
diff --git a/corsika/framework/process/ContinuousProcessIndex.hpp b/corsika/framework/process/ContinuousProcessIndex.hpp
index e36bf2812a33c6f03e4058033585f846a871b897..7192b86bf968ffdc74b5420940ea549befa12c73 100644
--- a/corsika/framework/process/ContinuousProcessIndex.hpp
+++ b/corsika/framework/process/ContinuousProcessIndex.hpp
@@ -11,12 +11,10 @@
namespace corsika {
/**
- @ingroup Processes
-
- To index individual processes (continuous processes) inside a
- ProcessSequence.
-
- **/
+ * @ingroup Processes
+ * To index individual processes (continuous processes) inside a
+ * ProcessSequence.
+ */
class ContinuousProcessIndex {
public:
diff --git a/corsika/framework/process/ContinuousProcessStepLength.hpp b/corsika/framework/process/ContinuousProcessStepLength.hpp
index 575be9063d49a32b5853320dc9d25bac40a7b9d8..55c2f417a470b8e4a475de8269c5e19fe8247ff5 100644
--- a/corsika/framework/process/ContinuousProcessStepLength.hpp
+++ b/corsika/framework/process/ContinuousProcessStepLength.hpp
@@ -14,12 +14,10 @@
namespace corsika {
/**
- @ingroup Processes
-
- To store step length in LengthType and unique index in ProcessSequence of shortest
- step ContinuousProcess.
-
- **/
+ * @ingroup Processes
+ * To store step length in LengthType and unique index in ProcessSequence of shortest
+ * step ContinuousProcess.
+ */
class ContinuousProcessStepLength {
public:
diff --git a/corsika/framework/utility/CubicSolver.hpp b/corsika/framework/utility/CubicSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0959a4fa84c3bee6aad856ccb57c887d21e30c9c
--- /dev/null
+++ b/corsika/framework/utility/CubicSolver.hpp
@@ -0,0 +1,93 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <algorithm>
+#include <numeric>
+#include <complex>
+
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+/**
+ * @file CubicSolver.hpp
+ */
+
+namespace corsika {
+
+ /**
+ * @ingroup Utility
+ * @{
+ */
+
+ namespace andre {
+ /**
+ Solve a x^3 + b x^2 + c x + d = 0
+ */
+ std::vector<double> solveP3(long double a, long double b, long double c,
+ long double d, double const epsilon = 1e-12);
+ } // namespace andre
+
+ /**
+ Solve depressed cubic: x^3 + p x + q = 0
+ */
+ inline std::vector<double> solve_cubic_depressed_disciminant_real(
+ long double p, long double q, long double const disc, double const epsilon = 1e-12);
+
+ std::vector<double> solve_cubic_depressed_real(long double p, long double q,
+ double const epsilon = 1e-12);
+
+ /**
+ Solve a x^3 + b x^2 + c x + d = 0
+
+ Analytical approach. Not very stable in all conditions.
+ */
+ std::vector<double> solve_cubic_real_analytic(long double a, long double b,
+ long double c, long double d,
+ double const epsilon = 1e-12);
+
+ /**
+ * Cubic function.
+ *
+ * T must be a floating point type.
+ *
+ * @return a x^3 + b x^2 + c x + d
+ */
+ template <typename T>
+ T cubic_function(T x, T a, T b, T c, T d);
+
+ /**
+ @return 3 a x^2 + 2 b x + c
+ */
+ template <typename T>
+ T cubic_function_dfdx(T x, T a, T b, T c);
+
+ /**
+ @return 6 a x + 2 b
+ */
+ template <typename T>
+ T cubic_function_d2fd2x(T x, T a, T b);
+
+ /**
+ Iterative approach to solve: a x^3 + b x^2 + c x + d = 0
+
+ This is often fastest and most precise.
+ */
+ std::vector<double> solve_cubic_real(long double a, long double b, long double c,
+ long double d, double const epsilon = 1e-12);
+
+ //! @}
+
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/CubicSolver.inl>
diff --git a/corsika/framework/utility/LinearSolver.hpp b/corsika/framework/utility/LinearSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..754039571287a1a816c77ec741b7d850b87c46da
--- /dev/null
+++ b/corsika/framework/utility/LinearSolver.hpp
@@ -0,0 +1,23 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <complex>
+
+namespace corsika {
+
+ std::vector<double> solve_linear_real(double a, double b);
+
+ std::vector<std::complex<double>> solve_linear(double a, double b);
+
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/LinearSolver.inl>
diff --git a/corsika/framework/utility/QuadraticSolver.hpp b/corsika/framework/utility/QuadraticSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..968c01bf13e3680b91b27005084c909b4e473ea7
--- /dev/null
+++ b/corsika/framework/utility/QuadraticSolver.hpp
@@ -0,0 +1,27 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <complex>
+
+#include <corsika/framework/utility/LinearSolver.hpp>
+
+namespace corsika {
+
+ std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
+ long double c,
+ double const epsilon = 1e-12);
+
+ std::vector<double> solve_quadratic_real(long double a, long double b, long double c,
+ double const epsilon = 1e-12);
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/QuadraticSolver.inl>
diff --git a/corsika/framework/utility/QuarticSolver.hpp b/corsika/framework/utility/QuarticSolver.hpp
index ad9484939d10691e38c7e74d9e2632b3d2e60463..d6e20f02fc3d6d594a64475d9ce1117b4dec159c 100644
--- a/corsika/framework/utility/QuarticSolver.hpp
+++ b/corsika/framework/utility/QuarticSolver.hpp
@@ -8,59 +8,46 @@
#pragma once
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+
#include <complex>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
/**
* @file QuarticSolver.hpp
- * @ingroup Utilities
- * @{
- * \todo convert to class
*/
-namespace corsika::quartic_solver {
-
- const double PI = 3.141592653589793238463L;
- const double M_2PI = 2 * PI;
- const double eps = 1e-12;
-
- typedef std::complex<double> DComplex;
-
- //---------------------------------------------------------------------------
- // useful for testing
- inline DComplex polinom_2(DComplex x, double a, double b) {
- // Horner's scheme for x*x + a*x + b
- return x * (x + a) + b;
- }
+namespace corsika {
- //---------------------------------------------------------------------------
- // useful for testing
- inline DComplex polinom_3(DComplex x, double a, double b, double c) {
- // Horner's scheme for x*x*x + a*x*x + b*x + c;
- return x * (x * (x + a) + b) + c;
- }
+ /**
+ * @ingroup Utilities
+ * @{
+ */
- //---------------------------------------------------------------------------
- // useful for testing
- inline DComplex polinom_4(DComplex x, double a, double b, double c, double d) {
- // Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
- return x * (x * (x * (x + a) + b) + c) + d;
- }
+ namespace andre {
- //---------------------------------------------------------------------------
- // x - array of size 3
- // In case 3 real roots: => x[0], x[1], x[2], return 3
- // 2 real roots: x[0], x[1], return 2
- // 1 real root : x[0], x[1] ± i*x[2], return 1
- unsigned int solveP3(double* x, double a, double b, double c);
+ /**
+ solve quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e
+ Attention - this function returns dynamically allocated array. It has to be
+ released afterwards.
+ */
+ std::vector<double> solve_quartic_real(long double a, long double b, long double c,
+ long double d, long double e,
+ double const epsilon = 1e-12);
+ } // namespace andre
- //---------------------------------------------------------------------------
- // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
- // Attention - this function returns dynamically allocated array. It has to be released
- // afterwards.
- DComplex* solve_quartic(double a, double b, double c, double d);
+ /**
+ solve quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e
+ */
+ std::vector<double> solve_quartic_real(long double a, long double b, long double c,
+ long double d, long double e,
+ double const epsilon = 1e-12);
//! @}
-} // namespace corsika::quartic_solver
+} // namespace corsika
#include <corsika/detail/framework/utility/QuarticSolver.inl>
diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index 8276c0d1235f44539bfc17468f20aff5eda902b5..f7ccba999c531333dbe8dedf9b7550c97f09ba19 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -20,9 +20,19 @@
namespace corsika {
/**
- * The ObservationPlane writes PDG codes, energies, and distances of particles to the
- * central point of the plane into its output file. The particles are considered
- * "absorbed" afterwards.
+ @ingroup Modules
+ @{
+
+ The ObservationPlane writes PDG codes, energies, and distances of particles to the
+ central point of the plane into its output file. The particles are considered
+ "absorbed" afterwards.
+
+ **Note/Limitation:** as discussed in
+ https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/397
+ you cannot put two ObservationPlanes exactly on top of each
+ other. Even if one of them is "permeable". You have to put a
+ small gap in between the two plane in such a scenario, or develop
+ another more specialized output class.
*/
class ObservationPlane : public ContinuousProcess<ObservationPlane> {
@@ -50,6 +60,7 @@ namespace corsika {
DirectionVector const xAxis_;
DirectionVector const yAxis_;
};
+ //! @}
} // namespace corsika
#include <corsika/detail/modules/ObservationPlane.inl>
diff --git a/corsika/modules/TrackWriter.hpp b/corsika/modules/TrackWriter.hpp
index 5bf9856b389fe80774ea629b27a4b4dc10f33c85..8b9b0efe7faca74ade9f9c55ea463186bc111e98 100644
--- a/corsika/modules/TrackWriter.hpp
+++ b/corsika/modules/TrackWriter.hpp
@@ -20,6 +20,7 @@ namespace corsika {
public:
TrackWriter(std::string const& filename);
+ ~TrackWriter();
template <typename TParticle, typename TTrack>
ProcessReturn doContinuous(TParticle const&, TTrack const&, bool const limitFlag);
diff --git a/corsika/modules/qgsjetII/Interaction.hpp b/corsika/modules/qgsjetII/Interaction.hpp
index a45282f58ae3073237adf19a3c0c739d6153047c..63c878df5d4cea4feadbe2ff028760a79643e741 100644
--- a/corsika/modules/qgsjetII/Interaction.hpp
+++ b/corsika/modules/qgsjetII/Interaction.hpp
@@ -13,6 +13,10 @@
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
+
+#include <boost/filesystem/path.hpp>
+
#include <qgsjet-II-04.hpp>
#include <string>
@@ -21,14 +25,8 @@ namespace corsika::qgsjetII {
class Interaction : public corsika::InteractionProcess<Interaction> {
- std::string data_path_;
- int count_ = 0;
- bool initialized_ = false;
- QgsjetIIHadronType alternate_ =
- QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
-
public:
- Interaction(const std::string& dataPath = "");
+ Interaction(boost::filesystem::path dataPath = corsika_data("QGSJetII"));
~Interaction();
bool wasInitialized() { return initialized_; }
@@ -53,6 +51,11 @@ namespace corsika::qgsjetII {
void doInteraction(TSecondaryView&);
private:
+ int count_ = 0;
+ bool initialized_ = false;
+ QgsjetIIHadronType alternate_ =
+ QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
+
corsika::default_prng_type& rng_ =
corsika::RNGManager::getInstance().getRandomStream("qgsjet");
const int maxMassNumber_ = 208;
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
index 3f6d3e9936f1fd13fcb831359f8faee6768f2948..69cb4a93750b3389749c279cdc0b6f8f84765ad5 100644
--- a/corsika/modules/urqmd/UrQMD.hpp
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -27,7 +27,7 @@ namespace corsika::urqmd {
class UrQMD : public InteractionProcess<UrQMD> {
public:
- UrQMD(boost::filesystem::path const& path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat"));
+ UrQMD(boost::filesystem::path const path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat"));
template <typename TParticle>
GrammageType getInteractionLength(TParticle const&) const;
@@ -46,7 +46,7 @@ namespace corsika::urqmd {
private:
static CrossSectionType getCrossSection(Code, Code, HEPEnergyType, int);
- void readXSFile(boost::filesystem::path const&);
+ void readXSFile(boost::filesystem::path);
// data members
default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("urqmd");
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index 34bd1377d99a5167b2ac406e506a4c8492173bc3..19e9a44c32aa401147d57585c89a7c1ab3ec6f4e 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -92,6 +92,10 @@ namespace corsika::nuclear_stack {
* Overwrite normal getParticleMass function with nuclear version
*/
HEPMassType getMass() const;
+ /**
+ * Overwrite normal getParticleCharge function with nuclear version
+ */
+ ElectricChargeType getCharge() const;
/**
* Overwirte normal getChargeNumber function with nuclear version
**/
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index adaae996b63a5b69efa85dc8ad8956e253de44b9..e5a765781ae44f38b38da7d4135a6363c697f5d3 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -87,8 +87,14 @@ namespace corsika {
return this->getMomentum() / this->getEnergy();
}
+ VelocityVector getVelocity() const {
+ return this->getMomentum() / this->getEnergy() * constants::c;
+ }
+
HEPMassType getMass() const { return get_mass(this->getPID()); }
+ ElectricChargeType getCharge() const { return get_charge(this->getPID()); }
+
int16_t getChargeNumber() const { return get_charge_number(this->getPID()); }
///@}
};
diff --git a/do-clang-format.py b/do-clang-format.py
index d6ea597535e69da32de80ff31962ba998523305d..fcfe674116832aea4ab94eb0c930b5d5c512fced 100755
--- a/do-clang-format.py
+++ b/do-clang-format.py
@@ -10,6 +10,7 @@ import os
import sys
import re
+debug = False
do_progress = False
try:
from progress.bar import ChargingBar
@@ -77,6 +78,9 @@ else:
filelist_clean.append(f)
filelist = filelist_clean
+if debug:
+ print ("filelist: ", filelist)
+
cmd = "clang-format"
if "CLANG_FORMAT" in os.environ:
cmd = os.environ["CLANG_FORMAT"]
@@ -99,11 +103,18 @@ if do_progress:
bar = ChargingBar('Processing', max=len(filelist))
if args.apply:
+ changed = []
for filename in filelist:
if bar: bar.next()
+ a = open(filename, "rb").read()
subp.check_call(cmd.split() + ["-i", filename])
+ b = open(filename, "rb").read()
+ if a != b:
+ changed.append(filename)
if bar: bar.finish()
-
+ if debug:
+ print ("changed: ", changed)
+
else:
# only print files which need formatting
files_need_formatting = 0
diff --git a/modules/qgsjetII/qgsjet-II-04.cpp b/modules/qgsjetII/qgsjet-II-04.cpp
index eea2724a88c07f1b028aba0703173a1ebc1f0504..8ffd24df4bd7b860116b7792e16671562af051b6 100644
--- a/modules/qgsjetII/qgsjet-II-04.cpp
+++ b/modules/qgsjetII/qgsjet-II-04.cpp
@@ -2,7 +2,7 @@
#include <iostream>
-datadir::datadir(const std::string& dir) {
+datadir::datadir(std::string const& dir) {
if (dir.length() > 130) {
std::cerr << "QGSJetII error, will cut datadir \"" << dir
<< "\" to 130 characters: " << std::endl;
@@ -13,11 +13,9 @@ datadir::datadir(const std::string& dir) {
data[i + 1] = '\0';
}
-
/**
@function qgran
link to random number generation
*/
-double qgran_(int&) { return qgsjetII::rndm_interface(); }
-
+double qgran_(int&) { return qgsjetII::rndm_interface(); }
diff --git a/src/framework/core/CMakeLists.txt b/src/framework/core/CMakeLists.txt
index 558b227f68801ec98c794c832ae5d8e654c1a703..f0edcdebd80fd61962c258b8e8ab43963030e5f6 100644
--- a/src/framework/core/CMakeLists.txt
+++ b/src/framework/core/CMakeLists.txt
@@ -5,7 +5,7 @@ file (MAKE_DIRECTORY ${output_dir})
add_custom_command (
OUTPUT ${output_dir}/GeneratedParticleProperties.inc
- OUTPUT ${output_dir}/GeneratedParticleClasses.inc
+ ${output_dir}/GeneratedParticleClasses.inc
${output_dir}/particle_db.pkl
COMMAND ${input_dir}/pdxml_reader.py ${input_dir}/ParticleData.xml
${input_dir}/NuclearData.xml
@@ -22,7 +22,9 @@ add_custom_command (
add_custom_target (GenParticlesHeaders
DEPENDS ${output_dir}/GeneratedParticleProperties.inc
- ${output_dir}/GeneratedParticleClasses.inc)
+ ${output_dir}/GeneratedParticleClasses.inc
+ ${output_dir}/particle_db.pkl
+ )
add_dependencies (CORSIKA8 GenParticlesHeaders)
install (
diff --git a/src/modules/qgsjetII/CMakeLists.txt b/src/modules/qgsjetII/CMakeLists.txt
index 36ebe9b693b927de9acfca03be3dc016440c4050..9b9777356e511db01238c2a68389c5c476e738a7 100644
--- a/src/modules/qgsjetII/CMakeLists.txt
+++ b/src/modules/qgsjetII/CMakeLists.txt
@@ -10,7 +10,7 @@ add_custom_command (
${input_dir}/qgsjet-II-04-codes.dat
DEPENDS ${input_dir}/code_generator.py
${input_dir}/qgsjet-II-04-codes.dat
- GenParticlesHeaders # for particle_db.pkl
+ ${PROJECT_BINARY_DIR}/corsika/framework/core/particle_db.pkl
WORKING_DIRECTORY
${output_dir}
COMMENT "Generate conversion tables for particle codes QGSJetII <-> CORSIKA"
diff --git a/src/modules/sibyll/CMakeLists.txt b/src/modules/sibyll/CMakeLists.txt
index b3d97a12349a5a5f2cbd9917404a8685894f0096..969d2f1d2e0e1bada1a4c247901ef381c121d368 100644
--- a/src/modules/sibyll/CMakeLists.txt
+++ b/src/modules/sibyll/CMakeLists.txt
@@ -10,7 +10,7 @@ add_custom_command (
${input_dir}/sibyll_codes.dat
DEPENDS ${input_dir}/code_generator.py
${input_dir}/sibyll_codes.dat
- GenParticlesHeaders # for particle_db.pkl
+ ${PROJECT_BINARY_DIR}/corsika/framework/core/particle_db.pkl
WORKING_DIRECTORY
${output_dir}/
COMMENT "Generate conversion tables for particle codes SIBYLL <-> CORSIKA"
diff --git a/tests/framework/CMakeLists.txt b/tests/framework/CMakeLists.txt
index 35f61737b13dfedef24d6746cf92e56a26658ebe..829709779bd7ec19b8b9fb0037422a1f6ca105da 100644
--- a/tests/framework/CMakeLists.txt
+++ b/tests/framework/CMakeLists.txt
@@ -20,6 +20,7 @@ set (test_framework_sources
testNullModel.cpp
testHelix.cpp
testInteractionCounter.cpp
+ testSolver.cpp
)
CORSIKA_ADD_TEST (testFramework SOURCES ${test_framework_sources})
diff --git a/tests/framework/testCOMBoost.cpp b/tests/framework/testCOMBoost.cpp
index caf6f08965bf4468294df23dee13e4f265dd4c47..d092a4b61ea46dfbcb734a4c8f31013550451204 100644
--- a/tests/framework/testCOMBoost.cpp
+++ b/tests/framework/testCOMBoost.cpp
@@ -40,7 +40,6 @@ auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
TEST_CASE("rotation") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
@@ -172,7 +171,6 @@ TEST_CASE("rotation") {
TEST_CASE("boosts") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// define target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
@@ -329,7 +327,6 @@ TEST_CASE("boosts") {
TEST_CASE("rest frame") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType const P0 = 1_TeV;
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index b860e9955d3c3122c4295cdaa6f32875f9aa6021..3e930911b095699bd09a10813c38c167c56b5dcc 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -137,7 +137,6 @@ public:
TEST_CASE("Cascade", "[Cascade]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
HEPEnergyType E0 = 100_GeV;
diff --git a/tests/framework/testClassTimer.cpp b/tests/framework/testClassTimer.cpp
index d51a56b1a254146847f742bfc0283ee9f6e44be2..b88d2fc97a0f2c99e5100268e168b6134598ca63 100644
--- a/tests/framework/testClassTimer.cpp
+++ b/tests/framework/testClassTimer.cpp
@@ -106,7 +106,6 @@ public:
TEST_CASE("ClassTimer", "[Timer]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Measure runtime of a function without arguments") {
diff --git a/tests/framework/testCombinedStack.cpp b/tests/framework/testCombinedStack.cpp
index 1f76fafc21927832063f5271fb8ed26443dffa50..de6defc6ad7a0f94fcbcc289a012de38761112dd 100644
--- a/tests/framework/testCombinedStack.cpp
+++ b/tests/framework/testCombinedStack.cpp
@@ -89,7 +89,6 @@ using StackTest = CombinedStack<TestStackData, TestStackData2, CombinedTestInter
TEST_CASE("Combined Stack", "[stack]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
@@ -395,7 +394,6 @@ using Particle2 = typename StackTest2::particle_type;
TEST_CASE("Combined Stack - secondary view") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("create secondaries via secondaryview") {
diff --git a/tests/framework/testCorsikaFenv.cpp b/tests/framework/testCorsikaFenv.cpp
index 2b6cdb58034d06c82cf336c2826cdb2e4f4248f6..030d33dfac78b92464b677b489f7f62b8fb6fa4c 100644
--- a/tests/framework/testCorsikaFenv.cpp
+++ b/tests/framework/testCorsikaFenv.cpp
@@ -22,7 +22,6 @@ static void handle_fpe(int /*signo*/) { gRESULT = 0; }
TEST_CASE("CorsikaFenv", "[fenv]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Enable all exceptions") { feenableexcept(FE_ALL_EXCEPT); }
diff --git a/tests/framework/testFourVector.cpp b/tests/framework/testFourVector.cpp
index 09b458295ef0b118684ed2707389a0a5f0fcfd9c..489d6261b17d6c68f85a01dae591992cbbba6a90 100644
--- a/tests/framework/testFourVector.cpp
+++ b/tests/framework/testFourVector.cpp
@@ -20,7 +20,6 @@ using namespace corsika;
TEST_CASE("four vectors") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::info);
// this is just needed as a baseline
diff --git a/tests/framework/testFunctionTimer.cpp b/tests/framework/testFunctionTimer.cpp
index 46d97f710f5793a5f27bbc4ed255c49cbdad627b..1328ea7badbc86f13809197a06ca773593381137 100644
--- a/tests/framework/testFunctionTimer.cpp
+++ b/tests/framework/testFunctionTimer.cpp
@@ -32,7 +32,6 @@ public:
TEST_CASE("FunctionTimer", "[Timer]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Measure runtime of a free function") {
diff --git a/tests/framework/testHelix.cpp b/tests/framework/testHelix.cpp
index cbe997d203a02705001f9304bc767cfedce848db..042e4212f122bb538c86e9c78317fb86259b3ba8 100644
--- a/tests/framework/testHelix.cpp
+++ b/tests/framework/testHelix.cpp
@@ -22,7 +22,6 @@ double constexpr absMargin = 1.0e-8;
TEST_CASE("Helix class") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
const CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
diff --git a/tests/framework/testInteractionCounter.cpp b/tests/framework/testInteractionCounter.cpp
index 556d4733b60b82bf295ff47463a24d3775c48553..4b4b7ff25dca0ba6a5cf44018c9efde66a150e82 100644
--- a/tests/framework/testInteractionCounter.cpp
+++ b/tests/framework/testInteractionCounter.cpp
@@ -43,7 +43,6 @@ struct DummyProcess {
TEST_CASE("InteractionCounter", "[process]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::info);
DummyProcess d;
diff --git a/tests/framework/testLogging.cpp b/tests/framework/testLogging.cpp
index b737d7dfa778b0cdf78b235479177e51534311b2..bf92cf8df00e8f12d64c0d5779f59b7976c4e09d 100644
--- a/tests/framework/testLogging.cpp
+++ b/tests/framework/testLogging.cpp
@@ -14,8 +14,6 @@ using namespace corsika;
TEST_CASE("Logging", "[Logging]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
SECTION("top level functions using default corsika logger") {
logging::info("(1) This is an info message!");
logging::warn("(1) This is a warning message!");
diff --git a/tests/framework/testNullModel.cpp b/tests/framework/testNullModel.cpp
index 5b1f351e13638a1ac2358bd5f0b9d2f295563350..872fc34dbbeb685db1806c74c93f4709443db545 100644
--- a/tests/framework/testNullModel.cpp
+++ b/tests/framework/testNullModel.cpp
@@ -21,7 +21,6 @@ using namespace corsika;
TEST_CASE("NullModel", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("interface") {
[[maybe_unused]] NullModel nm;
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index b74d836bc47d5715ec4f2504b30fb056af545af4..d8c71e6e787ce465d367c29dfa5d06b0b2a51faf 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -17,7 +17,6 @@ using namespace corsika;
TEST_CASE("ParticleProperties", "[Particles]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Types") {
CHECK(Electron::code == Code::Electron);
diff --git a/tests/framework/testRandom.cpp b/tests/framework/testRandom.cpp
index 5774575e6844754289cdcd2130313d3b03d95359..2d0c0720ae163f2777c80837d8179b9355aba044 100644
--- a/tests/framework/testRandom.cpp
+++ b/tests/framework/testRandom.cpp
@@ -55,7 +55,6 @@ SCENARIO("random-number streams can be registered and retrieved") {
TEST_CASE("UniformRealDistribution") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::mt19937 rng;
@@ -101,7 +100,6 @@ TEST_CASE("UniformRealDistribution") {
TEST_CASE("ExponentialDistribution") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::mt19937 rng;
diff --git a/tests/framework/testSaveBoostHistogram.cpp b/tests/framework/testSaveBoostHistogram.cpp
index d5c16438a2ab6de1595393205d4eff4320a43d08..89182311a650e76968226e3dc2cf5a97c8b4e29a 100644
--- a/tests/framework/testSaveBoostHistogram.cpp
+++ b/tests/framework/testSaveBoostHistogram.cpp
@@ -16,7 +16,6 @@ namespace bh = boost::histogram;
TEST_CASE("SaveHistogram") {
- corsika::corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
corsika::logging::set_level(corsika::logging::level::info);
std::mt19937 rng;
diff --git a/tests/framework/testSecondaryView.cpp b/tests/framework/testSecondaryView.cpp
index 87cfd91c0c74cb3b9c38cc8bda31af7a96e52c8c..fa113f88b6e111e9837d9ef96d28a432024974dd 100644
--- a/tests/framework/testSecondaryView.cpp
+++ b/tests/framework/testSecondaryView.cpp
@@ -44,7 +44,6 @@ using Particle = typename StackTest::particle_type;
TEST_CASE("SecondaryStack", "[stack]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
diff --git a/tests/framework/testSolver.cpp b/tests/framework/testSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..52f1fa13fc29f59d7b90e4471b76bb51997b0613
--- /dev/null
+++ b/tests/framework/testSolver.cpp
@@ -0,0 +1,278 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <cmath>
+#include <vector>
+#include <complex>
+#include <algorithm>
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+#include <corsika/framework/utility/LinearSolver.hpp>
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <corsika/framework/utility/QuarticSolver.hpp>
+
+using namespace corsika;
+using namespace std;
+
+double pol4(long double x, long double a, long double b, long double c, long double d,
+ long double e) {
+ return std::pow(x, 4) * a + std::pow(x, 3) * b + std::pow(x, 2) * c + x * d + e;
+}
+
+void remove_duplicates(vector<double>& v, double const eps) {
+ std::sort(v.begin(), v.end());
+ v.erase(std::unique(v.begin(), v.end(),
+ [eps](auto v1, auto v2) {
+ return (std::abs(v2) > eps ? std::abs(2 * (v1 - v2) / (v1 + v2)) <
+ eps // relative
+ : std::abs(v1 - v2) < eps); // absolute
+ }),
+ v.end());
+}
+
+TEST_CASE("Solver") {
+
+ logging::set_level(logging::level::trace);
+
+ double epsilon_check = 1e-3; // for catch2 asserts
+
+ SECTION("linear") {
+
+ std::vector<std::pair<double, double>> vals = {{13.33, 8338.3},
+ {-333.99, -8.15},
+ {-58633.3, 2343.54},
+ {87632.87, -982.37},
+ {7e8, -1e-7}};
+
+ for (auto v : vals) {
+
+ {
+ double a = v.first;
+ double b = v.second;
+
+ vector<double> s1 = solve_linear_real(a, b);
+
+ CORSIKA_LOG_INFO("linear: a={}, b={}, N={}, s1[0]={}", a, b, s1.size(), s1[0]);
+
+ double expected = -b / a;
+
+ CHECK(s1.size() == 1);
+ CHECK((s1[0] == Approx(expected).epsilon(epsilon_check)));
+
+ vector<complex<double>> s2 = solve_linear(a, b);
+ CHECK(s2.size() == 1);
+ CHECK((s2[0].real() == Approx(expected).epsilon(epsilon_check)));
+ }
+ }
+
+ CHECK(solve_linear_real(0, 55.).size() == 0);
+ }
+
+ SECTION("quadratic") {
+
+ // tests of type: (x-z1)*(x-z2) = 0 --> x1=z1, x2=z2 and a=1, b=-z2-z1,
+ // c=z1*z2
+
+ std::vector<std::vector<double>> zeros = {
+ {13.33, 8338.3}, {-333.99, -8.15}, {-58633.3, 2343.54}, {87632.87, -982.37},
+ {1.3e-6, 5e7}, {-4.2e-7, 65e6}, {0.1, -13e-5}, {7e8, -1e-7}};
+
+ for (unsigned int idegree = 0; idegree < 2; ++idegree) {
+ CORSIKA_LOG_INFO("------------- quadratic idegree={}", idegree);
+ for (auto z : zeros) {
+
+ {
+ long double const z1 = z[0];
+ long double const z2 = idegree <= 0 ? z1 : z[1];
+
+ long double a = 1;
+ long double b = -z1 - z2;
+ long double c = z1 * z2;
+
+ std::vector<double> s1 = solve_quadratic_real(a, b, c);
+ remove_duplicates(s1, epsilon_check);
+
+ CORSIKA_LOG_INFO("quadratic: z1={}, z2={}, N={}, s1=[{}]", z1, z2, s1.size(),
+ fmt::join(s1, ", "));
+
+ CHECK(s1.size() == idegree + 1);
+ for (auto value : s1) {
+ if (std::abs(value) < epsilon_check) {
+ CHECK(((value == Approx(z1).margin(epsilon_check)) ||
+ (value == Approx(z2).margin(epsilon_check))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
+ (value == Approx(z2).epsilon(epsilon_check))));
+ }
+ }
+ /*
+ std::vector<complex<double>> s2 = solve_quadratic(a, b, c, epsilon);
+ CHECK(s2.size() == idegree + 1);
+ for (auto value : s2) {
+ if (std::abs(value) < epsilon) {
+ CHECK(((value == Approx(z1).margin(epsilon_check)) ||
+ (value == Approx(z2).margin(epsilon_check))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
+ (value == Approx(z2).epsilon(epsilon_check))));
+ }
+ }
+ }*/
+ }
+ }
+ }
+ }
+
+ SECTION("cubic") {
+
+ // tests of type:
+ // (x-z1) * (x-z2) * (x-z3) = 0 --> x1=z1, x2=z2, x3=z3 and
+
+ // (x^2 - x*(z1+z2) + z1*z2) * (x-z3) =
+ // x^3 - x^2*(z1+z2) + x*z1*z2 - x^2 z3 + x*(z1+z2)*z3 - z1*z2*z3 =
+ // x^3 + x^2 (-z3-z1-z2) + x (z1*z2 + (z1+z2)*z3) - z1*z2*z3
+
+ epsilon_check = 1e-2; // for catch2 asserts
+
+ std::vector<std::vector<double>> zeros = {
+ {13.33, 838.3, 44.}, {-333.99, -8.15, -33e8}, {-58633.3, 2343.54, -1e-5},
+ {87632.87, -982.37, 0.}, {1.3e-4, 5e7, 6.6e9}, {-4.2e-7, 65e6, 433.3334},
+ {23e-1, -13e-2, 10.}, {7e8, -1e-7, 1e8}};
+
+ for (unsigned int idegree = 0; idegree < 3; ++idegree) {
+ CORSIKA_LOG_INFO("------------- cubic idegree={}", idegree);
+ for (auto z : zeros) {
+
+ {
+ long double const z1 = z[0];
+ long double const z2 = idegree <= 0 ? z1 : z[1];
+ long double const z3 = idegree <= 1 ? z1 : z[2];
+
+ long double const a = 1;
+ long double const b = -z1 - z2 - z3;
+ long double const c = z1 * z2 + (z1 + z2) * z3;
+ long double const d = -z1 * z2 * z3;
+ //
+ CORSIKA_LOG_INFO(
+ "cubic: z1={}, z2={}, z3={}, "
+ "a={}, b={}, c={}, d={}",
+ z1, z2, z3, a, b, c, d);
+
+ vector<double> s1 = solve_cubic_real(a, b, c, d);
+ remove_duplicates(s1, epsilon_check);
+
+ CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", "));
+
+ CHECK(s1.size() == idegree + 1);
+ for (double value : s1) {
+ CORSIKA_LOG_INFO("value={}, z1={} z2={} z3={} eps_check={}", value, z1, z2,
+ z3, epsilon_check);
+ if (std::abs(value) < epsilon_check) {
+ CHECK(((value == Approx(z1).margin(epsilon_check)) ||
+ (value == Approx(z2).margin(epsilon_check)) ||
+ (value == Approx(z3).margin(epsilon_check))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
+ (value == Approx(z2).epsilon(epsilon_check)) ||
+ (value == Approx(z3).epsilon(epsilon_check))));
+ }
+ }
+ }
+ }
+ }
+ } // cubic
+
+ SECTION("quartic") {
+
+ epsilon_check = 1e-4; // for catch2 asserts
+
+ // **clang-format-off**
+ // tests of type:
+ // (x-z1) * (x-z2) * (x-z3) * (x-z4) = 0 --> x1=z1, x2=z2, x3=z3, x4=z4 and
+
+ // (x^2 - x (z1+z2) + z1*z2) * (x-z3) *
+ // (x-z4) = (x^3 - x^2 (z1+z2) + x z1*z2 - x^2 z3 + x*(z1+z2)*z3 -
+ // z1*z2*z3) * (x-z4) = (x^3 + x^2 (-z1-z2-z3) + x (z1*z2 + (z1+z2)*z3) -
+ // z1*z2*z3) * (x-z4)
+ // =
+ //
+ // x^4 + x^3 (-z1-z2-z3) + x^2 (z1*z2 + (z1+z2)*z3) - x z1*z2*z3
+ // - x^3 z4 - x^2 (-z1-z2-z3)*z4 - x (z1*z2 + (z1+z2)*z3)*z4
+ // + z1*z2*z3*z4
+ //
+ // = x^4
+ // + x^3 (-z1-z2-z3-z4)
+ // + x^2 (z1*z2 + (z1+z2)*z3 + (z1+z2+z3)*z4)
+ // - x (z1*z2*z3 + (z1*z2 + (z1+z2)*z3)*z4))
+ // + z1*z2*z3*z4
+ // **clang-format-on**
+
+ std::vector<std::vector<double>> zeros = {
+ {13.33, 838.3, 44., 2.3}, {-3333.99, -8.15, -33e4, 8.8e3},
+ {-58633.3, 2343.54, -1e-1, 0.}, {87632.87, -982.37, 10., 1e-2},
+ {1.3e2, 5e5, 6.6e5, 1e3}, {-4.9, 65e2, 433.3334, 6e5},
+ {23e-1, -13e-2, 10., 3.4e6}, {7e6, -1e-1, 1e5, 2.55e4}};
+
+ for (unsigned int idegree = 2; idegree < 4;
+ ++idegree) { // idegree==1 is not very stable !!!!!!!
+ CORSIKA_LOG_INFO("------------- quartic idegree={}", idegree);
+ for (auto z : zeros) {
+
+ {
+ long double const z1 = z[0];
+ long double const z2 = idegree <= 0 ? z1 : z[1];
+ long double const z3 = idegree <= 1 ? z1 : z[2];
+ long double const z4 = idegree <= 2 ? z1 : z[3];
+
+ long double const a = 1;
+ long double const b = -z1 - z2 - z3 - z4;
+ long double const c = z1 * z2 + (z1 + z2) * z3 + (z1 + z2 + z3) * z4;
+ long double const d = -z1 * z2 * z3 - (z1 * z2 + (z1 + z2) * z3) * z4;
+ long double const e = z1 * z2 * z3 * z4;
+
+ //
+ CORSIKA_LOG_INFO(
+ "quartic: z1={}, z2={}, z3={}, z4={}, "
+ "a={}, b={}, c={}, d={}, e={}",
+ z1, z2, z3, z4, a, b, c, d, e);
+
+ // make sure the tested cases are all ZERO (printout)
+ CORSIKA_LOG_INFO("quartic trace: {} {} {} {}", pol4(z1, a, b, c, d, e),
+ pol4(z2, a, b, c, d, e), pol4(z3, a, b, c, d, e),
+ pol4(z4, a, b, c, d, e));
+
+ vector<double> s1 = andre::solve_quartic_real(a, b, c, d, e);
+ // vector<double> s1 = solve_quartic_real(a, b, c, d, e, epsilon);
+ remove_duplicates(s1, epsilon_check);
+
+ CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", "));
+
+ CHECK(s1.size() == idegree + 1);
+ for (double value : s1) {
+ CORSIKA_LOG_INFO("value={}, z1={} z2={} z3={} z4={} eps_check={}", value, z1,
+ z2, z3, z4, epsilon_check);
+ if (std::abs(value) < epsilon_check) {
+ CHECK(((value == Approx(z1).margin(epsilon_check)) ||
+ (value == Approx(z2).margin(epsilon_check)) ||
+ (value == Approx(z3).margin(epsilon_check)) ||
+ (value == Approx(z4).margin(epsilon_check))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
+ (value == Approx(z2).epsilon(epsilon_check)) ||
+ (value == Approx(z3).epsilon(epsilon_check)) ||
+ (value == Approx(z4).epsilon(epsilon_check))));
+ }
+ }
+ }
+ }
+ }
+ } // quartic
+}
diff --git a/tests/framework/testStackInterface.cpp b/tests/framework/testStackInterface.cpp
index b37dc495d7b03850ce135dd87942f8227531963d..4171647af5a547baeb55e6a2e3c3701eaeb77b29 100644
--- a/tests/framework/testStackInterface.cpp
+++ b/tests/framework/testStackInterface.cpp
@@ -27,7 +27,6 @@ typedef Stack<TestStackData, TestParticleInterface> StackTest;
TEST_CASE("Stack", "[Stack]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
diff --git a/tests/framework/testUnits.cpp b/tests/framework/testUnits.cpp
index 8d15544aca0d5c2dc1923af5853e6de7cafeb4f4..5c90aa9cd6ad65500a13a178d7dd6b0e52355006 100644
--- a/tests/framework/testUnits.cpp
+++ b/tests/framework/testUnits.cpp
@@ -19,7 +19,6 @@ using namespace corsika;
TEST_CASE("PhysicalUnits", "[Units]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Consistency") {
CHECK(1_m / 1_m == Approx(1));
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index 34718d69bb94e93c7e032e42df8f7e8e9c541f6e..2128e0a47a14141bcede691721a6dfea2da32d68 100644
--- a/tests/modules/testCONEX.cpp
+++ b/tests/modules/testCONEX.cpp
@@ -39,7 +39,6 @@ using MExtraEnvirnoment = MediumPropertyModel<UniformMagneticField<T>>;
TEST_CASE("CONEXSourceCut") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
RNGManager::getInstance().registerRandomStream("cascade");
RNGManager::getInstance().registerRandomStream("sibyll");
@@ -123,7 +122,6 @@ TEST_CASE("CONEXSourceCut") {
TEST_CASE("ConexOutput", "[output validation]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto file = GENERATE(as<std::string>{}, "conex_fit", "conex_output");
diff --git a/tests/modules/testExecTime.cpp b/tests/modules/testExecTime.cpp
index 68a85d81b2af7343ed323777daef39f55286a69e..4ebc03342914dc28e9137e7fd96dde010e8ac544 100644
--- a/tests/modules/testExecTime.cpp
+++ b/tests/modules/testExecTime.cpp
@@ -27,7 +27,6 @@ using namespace corsika::process::example_processors;
TEST_CASE("Timing process", "[proccesses][analytic_processors ExecTime]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
int tmp = 0;
diff --git a/tests/modules/testNullModel.cpp b/tests/modules/testNullModel.cpp
index a24433b080b4272955f01dfb2d116b43e0d5a2d6..6485b7df95a804e11301c0bc01095d705586fdea 100644
--- a/tests/modules/testNullModel.cpp
+++ b/tests/modules/testNullModel.cpp
@@ -14,7 +14,6 @@ using namespace corsika;
TEST_CASE("NullModel", "[processes]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
SECTION("interface") {
diff --git a/tests/modules/testOnShellCheck.cpp b/tests/modules/testOnShellCheck.cpp
index eccf73f287e9b07642ba5ed65e807d15d524f935..e28c5ffcd03a2d9eba64a890f5009ba9c9f18fb3 100644
--- a/tests/modules/testOnShellCheck.cpp
+++ b/tests/modules/testOnShellCheck.cpp
@@ -25,7 +25,6 @@ using namespace corsika;
TEST_CASE("OnShellCheck", "[processes]") {
logging::set_level(logging::level::debug);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
feenableexcept(FE_INVALID);
using EnvType = setup::Environment;
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index 1bc67084bc633c62539ec8a644eb42528a8e4106..0a2388772f960667cf99f7431089e5938ae1f66b 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -26,7 +26,6 @@ using namespace corsika;
TEST_CASE("ParticleCut", "processes") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
feenableexcept(FE_INVALID);
using EnvType = setup::Environment;
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index 8e178affc3f7c26c5a34df8caf2777386dba8d24..a4a93126a6c7456608f0cbd3f126e72910ae876d 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -23,7 +23,6 @@ using namespace corsika;
TEST_CASE("Pythia", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("linking pythia") {
using namespace Pythia8;
@@ -95,7 +94,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("PythiaInterface", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
auto const& cs = *csPtr;
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 4a604fd8bbf32f440db2eb0e073aed684b412fe4..3d68eb0cd4d14cd07c09bfd9bbc8afc0a2116350 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -49,7 +49,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("CORSIKA_DATA", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("check CORSIKA_DATA") {
@@ -68,7 +67,6 @@ TEST_CASE("CORSIKA_DATA", "[processes]") {
TEST_CASE("QgsjetII", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Corsika -> QgsjetII") {
CHECK(corsika::qgsjetII::convertToQgsjetII(PiMinus::code) ==
@@ -133,7 +131,6 @@ TEST_CASE("QgsjetII", "[processes]") {
TEST_CASE("QgsjetIIInterface", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env;
@@ -184,7 +181,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
corsika::qgsjetII::Interaction model;
model.doInteraction(view); // this also should produce some fragments
- CHECK(view.getSize() == Approx(188).margin(2)); // this is not physics validation
+ CHECK(view.getSize() == Approx(228).margin(100)); // this is not physics validation
int countFragments = 0;
for (auto const& sec : view) { countFragments += (sec.getPID() == Code::Nucleus); }
CHECK(countFragments == Approx(2).margin(1)); // this is not physics validation
@@ -219,6 +216,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Electron, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
GrammageType const length = model.getInteractionLength(particle);
@@ -228,7 +226,8 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Pi0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
- setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
CHECK(view.getSize() == Approx(18).margin(2)); // this is not physics validation
@@ -237,28 +236,31 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Rho0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
- setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
- CHECK(view.getSize() == Approx(7).margin(2)); // this is not physics validation
+ CHECK(view.getSize() == Approx(12).margin(8)); // this is not physics validation
}
{ // Lambda is internally converted into neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Lambda0, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
- setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
- CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+ CHECK(view.getSize() == Approx(15).margin(10)); // this is not physics validation
}
{ // AntiLambda is internally converted into anti neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Lambda0Bar, 0, 0, 100_GeV,
+ Code::Lambda0Bar, 0, 0, 1000_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
- setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
- CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+ CHECK(view.getSize() == Approx(40).margin(20)); // this is not physics validation
}
}
}
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index dc237b4a379613dff02be6802ca7047776aa45d4..dc5249a7218e2e2dea96b7a23f791833ef4e8c9a 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -32,7 +32,6 @@ using namespace corsika::sibyll;
TEST_CASE("Sibyll", "[processes]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
SECTION("Sibyll -> Corsika") {
@@ -99,7 +98,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("SibyllInterface", "[processes]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
index 00187ed9addbf6076641540dd494beccbd85758b..0b2a3ab7b9a489df955f00814ce371dfc4a026e9 100644
--- a/tests/modules/testStackInspector.cpp
+++ b/tests/modules/testStackInspector.cpp
@@ -24,7 +24,6 @@ using namespace corsika;
TEST_CASE("StackInspector", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto const& rootCS = get_root_CoordinateSystem();
Point const origin(rootCS, {0_m, 0_m, 0_m});
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index 87cd51e93baf36dd89430c230b0d97e9c5ee9357..b8574124227f477d5b2dfda969b10c4110585e95 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -20,32 +20,52 @@
#include <catch2/catch.hpp>
+#include <boost/type_index.hpp>
+
using namespace corsika;
+struct NonExistingDummyObject : public IVolume {
+ NonExistingDummyObject const& getVolume() const { return *this; }
+ bool contains(Point const&) const { return false; }
+};
+
template <typename T>
int sgn(T val) {
return (T(0) < val) - (val < T(0));
}
/**
- \file testTracking.cpp
+ @file testTracking.cpp
- This is the unified and commond unit test for all Tracking algorithms:
+ This is the unified and common unit test for all Tracking algorithms:
- tracking_leapfrog_curved::Tracking
- tracking_leapfrog_straight::Tracking
- tracking_line::Tracking
+
+ The main part of tests are to inject particles at 10GeV momentum at
+ (-Rg,0,0) in +x direction into a sphere of radius Rg, where Rg is
+ the gyroradius (or 10m for neutral particles). Then it is checked
+ where the particles leaves the sphere for different charges
+ (-1,0,+1) and field strength (-50uT, 0T, +50uT).
+
+ Each test is perfromed once, with the particles starting logically
+ outside of the Rg sphere (thus it first has to enter insides) and a
+ second time with the particle already logically inside the sphere.
+
+ There is a second smaller, internal sphere at +z displacement. Only
+ neutral particles are allowed and expected to hit this.
+
+ All those tests are parameterized, thus, they can be easily extended
+ or applied to new algorithms.
+
*/
-TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
- tracking_leapfrog_curved::Tracking,
+TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
- logging::set_level(logging::level::trace);
const HEPEnergyType P0 = 10_GeV;
@@ -69,17 +89,19 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
Bfield / 1_uT, outer)) {
CORSIKA_LOG_DEBUG(
- "********************\n TEST section PID={}, Bfield={} "
- "uT, outer (?)={}",
- PID, Bfield / 1_uT, outer);
+ "********************\n TEST algo={} section PID={}, "
+ "Bfield={} "
+ "uT, start_outside={}",
+ boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT, outer);
const int chargeNumber = get_charge_number(PID);
LengthType radius = 10_m;
int deflect = 0;
if (chargeNumber != 0 and Bfield != 0_T) {
deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection
- LengthType const gyroradius =
- P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(P0) *
+ constants::c / (abs(get_charge(PID)) * abs(Bfield)));
+ CORSIKA_LOG_DEBUG("Rg={} deflect={}", gyroradius, deflect);
radius = gyroradius;
}
@@ -91,6 +113,12 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
TestType tracking;
Point const center(cs, {0_m, 0_m, 0_m});
auto target = setup::Environment::createNode<Sphere>(center, radius);
+ // every particle should hit target_2
+ auto target_2 = setup::Environment::createNode<Sphere>(
+ Point(cs, {-radius * 3 / 4, 0_m, 0_m}), radius * 0.2);
+ // only neutral particles hit_target_neutral
+ auto target_neutral = setup::Environment::createNode<Sphere>(
+ Point(cs, {radius / 2, 0_m, 0_m}), radius * 0.1);
using MyHomogeneousModel = MediumPropertyModel<
UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
@@ -99,8 +127,18 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
target->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ target_neutral->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ target_2->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
auto* targetPtr = target.get();
+ auto* targetPtr_2 = target_2.get();
+ auto* targetPtr_neutral = target_neutral.get();
worldPtr->addChild(std::move(target));
+ targetPtr->addChild(std::move(target_2));
+ targetPtr->addChild(std::move(target_neutral));
auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs);
{ [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
@@ -123,18 +161,29 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
particle.setMomentum(traj.getDirection(1) * particle.getMomentum().getNorm());
if (outer) {
// now we know we are in target volume, depending on "outer"
- CHECK(traj.getLength(1) == 0_m);
+ CHECK(traj.getLength(1) / 1_m == Approx(0).margin(1e-3));
CHECK(nextVol == targetPtr);
}
// move forward, until we leave target volume
- while (nextVol == targetPtr) {
+ bool hit_neutral = false;
+ bool hit_2nd = false;
+ while (nextVol != worldPtr) {
+ if (nextVol == targetPtr_neutral) hit_neutral = true;
+ if (nextVol == targetPtr_2) hit_2nd = true;
const auto [traj2, nextVol2] = tracking.getTrack(particle);
nextVol = nextVol2;
particle.setNode(nextVol);
particle.setPosition(traj2.getPosition(1));
particle.setMomentum(traj2.getDirection(1) * particle.getMomentum().getNorm());
+ CORSIKA_LOG_TRACE("pos={}, p={}, |p|={} |v|={}, delta-l={}, delta-t={}",
+ particle.getPosition(), particle.getMomentum(),
+ particle.getMomentum().getNorm(),
+ particle.getVelocity().getNorm(), traj2.getLength(1),
+ traj2.getLength(1) / particle.getVelocity().getNorm());
}
CHECK(nextVol == worldPtr);
+ CHECK(hit_2nd == true);
+ CHECK(hit_neutral == (deflect == 0 ? true : false));
Point pointCheck(cs, (deflect == 0 ? radius : 0_m), (deflect * radius), 0_m);
@@ -147,3 +196,166 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
Approx(0).margin(1e-3));
}
}
+
+/** specifc test for curved leap-frog algorithm **/
+
+TEST_CASE("TrackingLeapFrogCurved") {
+
+ logging::set_level(logging::level::info);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ corsika::Code PID = Code::MuPlus;
+
+ using MyHomogeneousModel = MediumPropertyModel<
+ UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+ SECTION("infinite sphere / universe") {
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, 100_uT);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ tracking_leapfrog_curved::Tracking tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, worldPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in X-direction
+
+ particle.setPosition(Point(cs, 0_m, 0_m, 0_m));
+
+ Sphere sphere(center, 1_km * std::numeric_limits<double>::infinity());
+
+ auto intersections = tracking.intersect(particle, sphere);
+ // this must be a "linear trajectory" with no curvature
+
+ CHECK(intersections.hasIntersections() == false);
+ }
+
+ SECTION("momentum along field") {
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, 100_uT);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ tracking_leapfrog_curved::Tracking tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+ auto target = setup::Environment::createNode<Sphere>(center, 10_km);
+
+ MagneticFieldVector magneticfield(cs, 100_T, 0_T, 0_uT);
+ target->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ auto* targetPtr = target.get();
+ worldPtr->addChild(std::move(target));
+
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; } // prevent warning
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in X-direction
+
+ particle.setNode(targetPtr); // set particle outside "target" volume
+ particle.setPosition(Point(cs, 0_m, 0_m, 0_m));
+
+ auto [traj, nextVol] = tracking.getTrack(particle);
+ { [[maybe_unused]] auto const& dummy = nextVol; } // prevent warning
+
+ // this must be a "linear trajectory" with no curvature
+ CHECK(traj.getDirection(0).getComponents() == traj.getDirection(1).getComponents());
+ }
+}
+
+TEMPLATE_TEST_CASE("TrackingFail", "doesntwork", tracking_leapfrog_curved::Tracking,
+ tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
+
+ logging::set_level(logging::level::info);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, 1_mT);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ TestType tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+
+ auto [stack, viewPtr] =
+ setup::testing::setup_stack(Code::Proton, 0, 0, P0, worldPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ NonExistingDummyObject const dummy;
+ CHECK_THROWS(tracking.intersect(particle, dummy));
+}
+
+TEMPLATE_TEST_CASE("TrackingPlane", "plane", tracking_leapfrog_curved::Tracking,
+ tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
+
+ logging::set_level(logging::level::trace);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Gamma);
+ // for algorithms that know magnetic deflections choose: +-50uT, 0uT
+ // otherwise just 0uT
+ auto Bfield = GENERATE(filter(
+ []([[maybe_unused]] MagneticFluxType v) {
+ if constexpr (std::is_same_v<TestType, tracking_line::Tracking>)
+ return v == 0_uT;
+ else
+ return true;
+ },
+ values<MagneticFluxType>({50_uT, 0_uT, -50_uT})));
+
+ SECTION(fmt::format("Tracking PID={}, Bfield={} uT", PID, Bfield / 1_uT)) {
+
+ CORSIKA_LOG_DEBUG(
+ "********************\n TEST algo={} section PID={}, "
+ "Bfield={}uT ",
+ boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT);
+
+ const int chargeNumber = get_charge_number(PID);
+ LengthType radius = 10_m;
+ int deflect = 0;
+ if (chargeNumber != 0 and Bfield != 0_T) {
+ deflect = 1;
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(P0) *
+ constants::c / (abs(get_charge(PID)) * abs(Bfield)));
+ CORSIKA_LOG_DEBUG("Rg={} deflect={}", gyroradius, deflect);
+ radius = gyroradius;
+ }
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, Bfield);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ TestType tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, worldPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in Z-direction
+
+ particle.setPosition(Point(cs, -radius, 0_m, 0_m));
+
+ // put plane where we expect deflection by radius/2
+ Plane const plane(Point(cs, radius * (1 - sqrt(3. / 4.)), 0_m, 0_m),
+ DirectionVector(cs, {-1., 0., 0.}));
+ Intersections const hit = tracking.intersect(particle, plane);
+
+ CORSIKA_LOG_DEBUG("entry={} exit={}", hit.getEntry(), hit.getExit());
+
+ CHECK(hit.hasIntersections());
+ CHECK(hit.getEntry() / 1_s == Approx(0.00275 * deflect).margin(0.0003));
+ CHECK(hit.getExit() == 1_s * std::numeric_limits<double>::infinity());
+ }
+}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 0f9f4a86243ec0f7f5b0c16e989b9459a5e0a293..f55aa032948c22026b170776d8946dc074c6e688 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -59,7 +59,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("UrQMD") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("conversion") {
CHECK_THROWS(corsika::urqmd::convertFromUrQMD(106, 0));
@@ -99,6 +98,7 @@ TEST_CASE("UrQMD") {
{ [[maybe_unused]] auto const& env_dummy = env; }
auto [stack, view] =
setup::testing::setup_stack(Code::Proton, 0, 0, 100_GeV, nodePtr, cs);
+ [[maybe_unused]] setup::StackView& viewRef = *(view.get());
CHECK(urqmd.getInteractionLength(stack->getNextParticle()) / 1_g * square(1_cm) ==
Approx(105).margin(5));
}
@@ -109,6 +109,7 @@ TEST_CASE("UrQMD") {
{ [[maybe_unused]] auto const& env_dummy = env; }
auto [stack, view] =
setup::testing::setup_stack(Code::Neutron, 0, 0, 100_GeV, nodePtr, cs);
+ [[maybe_unused]] setup::StackView& viewRef = *(view.get());
CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle()));
}
@@ -121,6 +122,7 @@ TEST_CASE("UrQMD") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Nucleus, A, Z, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
+ [[maybe_unused]] setup::StackView& viewRef = *(secViewPtr.get());
CHECK(stackPtr->getEntries() == 1);
CHECK(secViewPtr->getEntries() == 0);
@@ -172,6 +174,7 @@ TEST_CASE("UrQMD") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
+ [[maybe_unused]] auto particle = stackPtr->first();
CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); // Code::Proton not a valid target
}
@@ -227,22 +230,4 @@ TEST_CASE("UrQMD") {
CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
Approx(0).margin(1e-2));
}
-
- SECTION("K0Short projectile") {
- auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
- [[maybe_unused]] auto const& env_dummy = env; // against warnings
- [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
-
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::K0Short, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
- CHECK(stackPtr->getEntries() == 1);
- CHECK(secViewPtr->getEntries() == 0);
-
- // must be assigned to variable, cannot be used as rvalue?!
- auto projectile = secViewPtr->getProjectile();
- auto const projectileMomentum = projectile.getMomentum();
-
- CHECK_THROWS(urqmd.doInteraction(*secViewPtr));
- }
}