From 8463b28ad807edb61748a919c613f263538e924e Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Fri, 23 Apr 2021 09:22:03 +0200
Subject: [PATCH] rebase to master
---
.../detail/framework/utility/CubicSolver.inl | 3 ---
.../framework/utility/QuarticSolver.inl | 24 +++++++++----------
.../modules/tracking/TrackingStraight.inl | 3 +--
examples/vertical_EAS.cpp | 8 ++-----
4 files changed, 14 insertions(+), 24 deletions(-)
diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
index efeb48c34..0b73a606f 100644
--- a/corsika/detail/framework/utility/CubicSolver.inl
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -33,7 +33,6 @@ namespace corsika {
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 r2 = r * r;
long double q3 = q * q * q;
// disc = q**3 + r**2
@@ -56,8 +55,6 @@ namespace corsika {
// sum all terms into final result
long double const disc = -(((e + uf) + au) + ab) + v;
- CORSIKA_LOG_TRACE("solveP3 disc={} r2-r3={}", disc, r2 - q3);
-
if (disc >= 0) {
long double t = r / std::sqrt(q3);
if (t < -1) t = -1;
diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
index 2bed4396c..fc613fc86 100644
--- a/corsika/detail/framework/utility/QuarticSolver.inl
+++ b/corsika/detail/framework/utility/QuarticSolver.inl
@@ -61,9 +61,9 @@ namespace corsika {
}
} else {
if (Det < 0) return {};
- long double sqDet = sqrt(Det);
- q1 = (y + sqDet) * 0.5;
- q2 = (y - sqDet) * 0.5;
+ 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);
@@ -71,32 +71,30 @@ namespace corsika {
// solving quadratic eqs. x^2 + p1*x + q1 = 0
// x^2 + p2*x + q2 = 0
+#ifdef DEBUG
{
- long double Det1 = p1 * p1 - 4 * q1;
- long double Det2 = p2 * p2 - 4 * q2;
+ [[maybe_unused]] long double Det1 = p1 * p1 - 4 * q1;
+ [[maybe_unused]] long double Det2 = p2 * p2 - 4 * q2;
CORSIKA_LOG_TRACE("Det1={} Det2={}", Det1, Det2);
if (Det1 >= 0 && Det2 >= 0) {
- long double sqDet1 = sqrt(Det1);
- long double sqDet2 = sqrt(Det2);
- // return {double(-p1 + sqDet1) * 0.5, double(-p1 - sqDet1) * 0.5, double(-p2 +
- // sqDet2) * 0.5, double(-p2 - sqDet2) * 0.5};
+ [[maybe_unused]] long double sqDet1 = sqrt(Det1);
+ [[maybe_unused]] long double sqDet2 = sqrt(Det2);
CORSIKA_LOG_TRACE("check1 {} {} {} {}", double(-p1 + sqDet1) * 0.5,
double(-p1 - sqDet1) * 0.5, double(-p2 + sqDet2) * 0.5,
double(-p2 - sqDet2) * 0.5);
}
if (Det1 >= 0) {
- long double sqDet1 = sqrt(Det1);
- // return {double(-p1 + sqDet1) * 0.5, double(-p1 - sqDet1) * 0.5};
+ [[maybe_unused]] long double sqDet1 = sqrt(Det1);
CORSIKA_LOG_TRACE("check2 {} {} ", double(-p1 + sqDet1) * 0.5,
double(-p1 - sqDet1) * 0.5);
}
if (Det2 >= 0) {
- long double sqDet2 = sqrt(Det2);
- // return {double(-p2 + sqDet2) * 0.5, double(-p2 - sqDet2) * 0.5};
+ [[maybe_unused]] long double sqDet2 = sqrt(Det2);
CORSIKA_LOG_TRACE("check3 {} {}", double(-p1 + sqDet2) * 0.5,
double(-p1 - sqDet2) * 0.5);
}
}
+#endif
std::vector<double> quad1 = solve_quadratic_real(1, p1, q1);
std::vector<double> quad2 = solve_quadratic_real(1, p2, q2);
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 7f4129188..2799cde00 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -65,8 +65,7 @@ namespace corsika::tracking_line {
auto const sqDisc = sqrt(discriminant);
auto const invDenom = 1 / vSqNorm;
- bool const numericallyInside = sphere.contains(position);
- CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
+ CORSIKA_LOG_TRACE("numericallyInside={}", sphere.contains(position));
return Intersections((-vDotDelta - sqDisc) * invDenom,
(-vDotDelta + sqDisc) * invDenom);
}
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index ae123f4cf..1459b709a 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -103,7 +103,7 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("vertical_EAS");
if (argc < 5) {
- std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> <Nevt> [seed] [output-dir] \n"
+ std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> <Nevt> [seed] \n"
" if A=0, Z is interpreted as PDG code \n"
" if no seed is given, a random seed is chosen \n"
<< std::endl;
@@ -116,8 +116,6 @@ int main(int argc, char** argv) {
if (argc > 5) { seed = std::stoi(std::string(argv[5])); }
- CORSIKA_LOG_INFO("output_dir={} seed={}", output_dir, seed);
-
// initialize random number sequence(s)
registerRandomStreams(seed);
@@ -400,7 +398,7 @@ int main(int argc, char** argv) {
BetheBlochPDG emContinuous(showerAxis);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
- TrackWriter trackWriter(tracks_dir);
+ TrackWriter trackWriter(tracks_file);
LongitudinalProfile longprof{showerAxis};
@@ -433,8 +431,6 @@ int main(int argc, char** argv) {
cut.reset();
// emContinuous.reset();
- longprof.save(longprof_dat);
-
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
--
GitLab