diff --git a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
index 9ff574321434aa6c7c495ea3fc8db432e4a47879..af31f8209eedad47bd9de9cbf1344c17c9ef758e 100644
--- a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
+++ b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
@@ -12,6 +12,7 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/QuarticSolver.hpp>
namespace corsika {
@@ -35,27 +36,49 @@ namespace corsika {
}
inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
- return initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+ return getDirection(u) * initialVelocity_.getNorm();
}
inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
- return getVelocity(u).normalized();
+ return (initialDirection_ +
+ initialDirection_.cross(magneticfield_) * timeStep_ * u * k_)
+ .normalized();
}
inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
- return u * timeStep_ *
- (double(getVelocity(u).getNorm() / initialVelocity_.getNorm()) + 1.0) / 2;
+ TimeType const step = timeStep_ * u;
+ double const correction = 1;
+ // the eventual (delta-L to L) correction factor is:
+ // (initialDirection_ + initialDirection_.cross(magneticfield_) * step *
+ // k_).getNorm();
+ return step / 2 * (correction + 1);
}
inline LengthType LeapFrogTrajectory::getLength(double const u) const {
- return timeStep_ * initialVelocity_.getNorm() * u;
+ return getDuration(u) * initialVelocity_.getNorm();
}
inline void LeapFrogTrajectory::setLength(LengthType const limit) {
- if (initialVelocity_.getNorm() == 0_m / 1_s) setDuration(0_s);
+ if (initialVelocity_.getNorm() == SpeedType::zero()) setDuration(0_s);
setDuration(limit / initialVelocity_.getNorm());
}
- inline void LeapFrogTrajectory::setDuration(TimeType const limit) { timeStep_ = limit; }
+ inline void LeapFrogTrajectory::setDuration(TimeType const limit) {
+ /*
+ initial attempt to calculate delta-L from assumed full-leap-frog-length L:
+
+ Note: often return 0. Not good enough yet.
+
+ LengthType const L = initialVelocity_.getNorm() * limit; // distance
+ double const a = (initialVelocity_.cross(magneticfield_) * k_).getSquaredNorm() / 4 /
+ square(1_m) * static_pow<4>(1_s);
+ double const d = L * initialVelocity_.getNorm() / square(1_m) * 1_s;
+ double const e = -square(L) / square(1_m);
+ std::vector<double> solutions = solve_quartic_real(a, 0, 0, d, e);
+ CORSIKA_LOG_DEBUG("setDuration limit={} L={} solution={}", limit, L,
+ fmt::join(solutions, ", "));
+ */
+ timeStep_ = limit;
+ }
} // namespace corsika
diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
index c15594e02a02bef2683e4d3ae4d4ad02895ef747..6baa0a39d3515453158c42ce09d580763d845882 100644
--- a/corsika/detail/framework/utility/CubicSolver.inl
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -18,6 +18,8 @@ namespace corsika {
//---------------------------------------------------------------------------
// solve cubic equation A x^3 + B*x^2 + C*x + D = 0
+ // x^3 + a*x^2 + b*x + c = 0
+ // mainly along WolframAlpha formulas
inline std::vector<double> solve_cubic_real_analytic(long double A, long double B,
long double C, long double D,
double const epsilon) {
@@ -29,8 +31,8 @@ namespace corsika {
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 q = (3 * b - a2) / 9;
+ long double r = (a * (9 * b - 2 * a2) - 27 * c) / 54;
long double q3 = q * q * q;
// disc = q**3 + r**2
@@ -38,8 +40,8 @@ namespace corsika {
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);
+ 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;
@@ -51,30 +53,34 @@ namespace corsika {
au = f - au;
ab = u - ab;
// sum all terms into final result
- long double const disc = -(((e + uf) + au) + ab) + v;
+ long double const disc = (((e + uf) + au) + ab) + v;
- if (disc >= 0) {
- long double t = r / std::sqrt(q3);
+ CORSIKA_LOG_TRACE("disc={} {}", disc, q3 + r * r);
+
+ if (std::abs(disc) < epsilon) {
+
+ a /= 3;
+ long double const cbrtR = std::cbrt(r);
+ return {double(2 * cbrtR - a), double(-cbrtR - a)}; // 2nd solution is doublet
+
+ } else if (disc > 0) {
+
+ long double const S = std::cbrt(r + std::sqrt(disc));
+ long double const T = std::cbrt(r - std::sqrt(disc));
+ a /= 3;
+ return {double((S + T) - a)}; // plus two imaginary solution
+
+ } else { // 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);
+ 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)};
+ double(q * std::cos((t + 4 * M_PI) / 3) - a)};
}
}
} // namespace andre
@@ -228,67 +234,95 @@ namespace corsika {
}
/**
- * Iterative approach.
+ * Iterative approach. https://en.wikipedia.org/wiki/Halley%27s_method
+ * Halley's method
*/
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),
+ CORSIKA_LOG_TRACE("cubic_iterative: 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);
+ auto pre_opt = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
+ long double x1 = 0; // start value
- 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);
+ if (pre_opt.size()) {
+ x1 = pre_opt[0]; //*std::max_element(pre_opt.begin(), pre_opt.end());
+#ifdef DEBUG
+ for (long double test_v : pre_opt) {
+ CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
+ cubic_function(test_v, a, b, c, d));
}
-
- 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);
+#endif
+ } else {
+ // this is only if the former solve_cubic_real_analytic would not result
+ // in any solution. We have no test case for this. This is excluded from tests:
+ // LCOV_EXCL_START
+ long double const dist = std::fma(b / a, b / a, -3 * c / a);
+ long double const xinfl = -b / (a * 3);
+
+ x1 = xinfl;
+ long double f_test = cubic_function(xinfl, a, b, c, d);
+
+ if (std::abs(f_test) > epsilon) {
+ if (std::abs(dist) < epsilon) {
+ x1 = xinfl - std::cbrt(f_test);
+ } else if (dist > 0) {
+ if (f_test > 0)
+ x1 = xinfl - 2 / 3 * std::sqrt(dist);
+ else
+ x1 = xinfl + 2 / 3 * std::sqrt(dist);
}
- 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));
+ }
+ // LCOV_EXCL_STOP
+ }
- CORSIKA_LOG_TRACE("niter={}", niter);
+ long double f_x1 = cubic_function(x1, a, b, c, d);
+ long double dx1 = 0;
+
+ 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) {
+ dx1 = std::cbrt(f_x1);
+ } else {
+ dx1 = f_x1 * f_prime_x1 * 2 / (f_prime_x1 * f_prime_x1 * 2 - f_x1 * f_prime2_x1);
+ }
+ x1 -= dx1;
+ f_x1 = cubic_function(x1, a, b, c, d);
+ CORSIKA_LOG_TRACE(
+ "niter={} x1={:.20f} f_x1={:.20f} f_prime={:.20f} f_prime2={:.20f} dx1={}",
+ niter, x1, f_x1, f_prime_x1, f_prime2_x1,
+ f_x1 * f_prime_x1 / (f_prime_x1 * f_prime_x1 - f_x1 * f_prime2_x1 * 0.5));
+ } while ((++niter < maxiter) && (std::abs(f_x1) > epsilon * 1000) &&
+ (std::abs(dx1) > epsilon));
+
+ CORSIKA_LOG_TRACE("niter={}", niter);
+ if (niter >= maxiter) {
+ CORSIKA_LOG_DEBUG("niter reached max iterations {}", niter);
+ return andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
}
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);
+ std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, 1e-3);
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
} // namespace corsika
diff --git a/corsika/detail/framework/utility/QuadraticSolver.inl b/corsika/detail/framework/utility/QuadraticSolver.inl
index 9c9759e975ae6d322dee7a478f935c095bdf0ff9..a335b00d682f540b48fa549f2f496af3efdbf2e0 100644
--- a/corsika/detail/framework/utility/QuadraticSolver.inl
+++ b/corsika/detail/framework/utility/QuadraticSolver.inl
@@ -62,7 +62,7 @@ namespace corsika {
long double f = std::fma(b, b, -w);
long double radicant = f + e;
- CORSIKA_LOG_TRACE(" radicant={} {} ", radicant, b * b - a * c * 4);
+ CORSIKA_LOG_TRACE("radicant={} {} ", radicant, b * b - a * c * 4);
if (std::abs(radicant) < epsilon) { // just one real solution
return {double(-b / (2 * a))};
diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
index 4562b39c0aa0fea4330fc686deab57d48683d98f..81bda646e5b21aba7534df5265cf329997e35ed0 100644
--- a/corsika/detail/framework/utility/QuarticSolver.inl
+++ b/corsika/detail/framework/utility/QuarticSolver.inl
@@ -35,6 +35,9 @@ namespace corsika {
// y^3 − c*y^2 + (bd−4e)*y − b^2*e−d^2+4*c*e = 0
std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
+ if (!x3.size()) {
+ return {}; // no solution, numeric problem (LCOV_EXCL_LINE)
+ }
long double y = x3[0]; // there is always at least one solution
// The essence - choosing Y with maximal absolute value.
if (x3.size() == 3) {
@@ -74,8 +77,8 @@ namespace corsika {
// 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);
+ std::vector<double> quad1 = solve_quadratic_real(1, p1, q1, 1e-5);
+ std::vector<double> quad2 = solve_quadratic_real(1, p2, q2, 1e-5);
if (quad2.size() > 0) {
for (auto val : quad2) quad1.push_back(val);
}
@@ -109,17 +112,37 @@ namespace corsika {
}
CORSIKA_LOG_TRACE("check m={}", m);
if (m == 0) { return {0}; }
-
- CORSIKA_LOG_TRACE("check m={}", m);
-
+ if (m < 0) {
+ // this is a rare numerical instability
+ // first: try analytical solution, second: discard (curved->straight tracking)
+ std::vector<double> const resolve_cubic_analytic =
+ andre::solve_cubic_real_analytic(1, p, p2 / 4 - r, -q2 / 8, epsilon);
+
+ CORSIKA_LOG_TRACE("andre::resolve_cubic_analytic: N={}, m=[{}]",
+ resolve_cubic_analytic.size(),
+ fmt::join(resolve_cubic_analytic, ", "));
+
+ if (!resolve_cubic_analytic.size()) return {};
+
+ for (auto const& v : resolve_cubic_analytic) {
+ 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}; }
+ if (m < 0) {
+ return {}; // now we are out of options, cannot solve: curved->straight tracking
+ }
+ }
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);
std::vector<double> z_quad1 =
- solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, epsilon);
+ solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, 1e-5);
std::vector<double> z_quad2 =
- solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, epsilon);
+ solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, 1e-5);
for (auto const& z : z_quad2) z_quad1.push_back(z);
return z_quad1;
}
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 967b9314198af14e7eb76b9e99d654d60c00c573..b3d4b916c5ea6af18ea40d23166d296ae173095b 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -27,16 +27,22 @@ namespace corsika {
The current step did not yet reach the ObservationPlane, do nothing now and wait:
*/
if (!stepLimit) {
-#ifdef DEBUG
+ // @todo this is actually needed to fix small instabilities of the leap-frog
+ // tracking: Note, this is NOT a general solution and should be clearly revised with
+ // a more robust tracking. #ifdef DEBUG
if (deleteOnHit_) {
+ // since this is basically a bug, it cannot be tested LCOV_EXCL_START
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;
+ CORSIKA_LOG_WARN("PARTICLE AVOIDED OBSERVATIONPLANE {}", check);
+ CORSIKA_LOG_WARN("Temporary fix: write and remove particle.");
+ } else
+ return ProcessReturn::Ok;
+ // LCOV_EXCL_STOP
+ } else
+ // #endif
+ return ProcessReturn::Ok;
}
HEPEnergyType const energy = particle.getEnergy();
@@ -63,8 +69,9 @@ namespace corsika {
inline LengthType ObservationPlane<TTracking, TOutput>::getMaxStepLength(
TParticle const& particle, TTrajectory const& trajectory) {
- CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}", particle.asString(),
- particle.getPosition(), particle.getDirection(), plane_.asString());
+ CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, plane={}",
+ particle.asString(), particle.getPosition(),
+ particle.getDirection(), plane_.asString());
auto const intersection = TTracking::intersect(particle, plane_);
@@ -77,10 +84,10 @@ namespace corsika {
return std::numeric_limits<double>::infinity() * 1_m;
}
double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
- auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
- auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm();
- CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m", dist / 1_m);
- return dist;
+ CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength dist={} m, pos={}",
+ trajectory.getLength(fractionOfIntersection) / 1_m,
+ trajectory.getPosition(fractionOfIntersection));
+ return trajectory.getLength(fractionOfIntersection);
}
template <typename TTracking, typename TOutput>
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index a2f643cc58d5f72f77e64579157e9a351aa411ed..7772f26609fdd65a3a1216235f4788ba5eda31a8 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -73,7 +73,7 @@ namespace corsika {
// Sternheimer parameterization, density corrections towards high energies
// NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
// MISSING
- CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m{}=",
+ CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m={}",
p.getMomentum().getNorm() / m);
double const x = log10(p.getMomentum().getNorm() / m);
double delta = 0;
diff --git a/corsika/detail/modules/tracking/Intersect.inl b/corsika/detail/modules/tracking/Intersect.inl
index c653cc5aa698107cd25eb3c009e634f5f9bd8f8a..e627213515fcae2c97635931cd7e3a338afaa6ad 100644
--- a/corsika/detail/modules/tracking/Intersect.inl
+++ b/corsika/detail/modules/tracking/Intersect.inl
@@ -44,9 +44,10 @@ namespace corsika {
fmt::ptr(&volumeNode), fmt::ptr(volumeNode.getParent()),
time_intersections_curr.hasIntersections());
if (time_intersections_curr.hasIntersections()) {
- CORSIKA_LOG_DEBUG("intersection times with currentLogicalVolumeNode: {} s and {} s",
- time_intersections_curr.getEntry() / 1_s,
- time_intersections_curr.getExit() / 1_s);
+ CORSIKA_LOG_DEBUG(
+ "intersection times with currentLogicalVolumeNode: entry={} s and exit={} s",
+ time_intersections_curr.getEntry() / 1_s,
+ time_intersections_curr.getExit() / 1_s);
if (time_intersections_curr.getExit() <= minTime) {
minTime =
time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 1ede8cc9825bb926c6423511dcc701005cc26a40..16b8369674bd7f11fc31a70d04df1891e444821d 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -29,10 +29,10 @@ namespace corsika {
namespace tracking_leapfrog_curved {
template <typename TParticle>
- inline auto make_LeapFrogStep(TParticle const& particle, LengthType steplength) {
+ inline auto Tracking::makeStep(TParticle const& particle, LengthType steplength) {
if (particle.getMomentum().getNorm() == 0_GeV) {
return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
- double(0));
+ double(0) * 1_m);
} // charge of the particle
ElectricChargeType const charge = particle.getCharge();
auto const* currentLogicalVolumeNode = particle.getNode();
@@ -40,9 +40,14 @@ namespace corsika {
currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
VelocityVector velocity = particle.getVelocity();
- decltype(meter / (second * volt)) k =
- charge * constants::cSquared * 1_eV /
- (velocity.getNorm() * 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 * meter)) const k =
+ charge / p_norm; // * initialVelocity.getNorm();
+
DirectionVector direction = velocity.normalized();
auto position = particle.getPosition(); // First Movement
// assuming magnetic field does not change during movement
@@ -94,6 +99,8 @@ namespace corsika {
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
+ CORSIKA_LOG_TRACE("p_perp={} eV", p_perp / 1_eV);
+
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.
@@ -101,10 +108,20 @@ namespace corsika {
return getLinearTrajectory(particle);
}
- LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
- constants::c / (abs(charge) * magnitudeB));
+ LengthType const gyroradius = convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
+ constants::c / (abs(charge) * magnitudeB);
- double const maxRadians = 0.01; // maximal allowed deflection
+ if (gyroradius > 1e9_m) {
+ // this cannot be really unit-tested. It is hidden. LCOV_EXCL_START
+ CORSIKA_LOG_WARN(
+ "CurvedLeapFrog is not very stable for extremely high gyroradius steps. "
+ "Rg={} -> straight tracking.",
+ gyroradius);
+ return getLinearTrajectory(particle);
+ // LCOV_EXCL_STOP
+ }
+
+ double const maxRadians = 0.01; // maximally allowed deflection
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
TimeType const steplimit_time = steplimit / initialVelocity.getNorm();
CORSIKA_LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
@@ -121,19 +138,20 @@ namespace corsika {
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)
+ // 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
- minNode); // next volume node
+ minTime), // --> trajectory
+ minNode); // --> next volume node
}
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Sphere const& sphere) {
- if (sphere.getRadius() == 1_km * std::numeric_limits<double>::infinity()) {
+ LengthType const radius = sphere.getRadius();
+ if (radius == 1_km * std::numeric_limits<double>::infinity()) {
return Intersections();
}
@@ -150,36 +168,68 @@ namespace corsika {
auto const projectedDirectionSqrNorm = projectedDirection.getSquaredNorm();
bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
- if ((charge == 0 * constants::e) || magneticfield.getNorm() == 0_T || isParallel) {
+ CORSIKA_LOG_TRACE("projectedDirectionSqrNorm={} T^2",
+ projectedDirectionSqrNorm / square(1_T));
+ if (isParallel) {
+ // particle moves parallel to field -> no deflection
return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
}
bool const numericallyInside = sphere.contains(particle.getPosition());
CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
+ Vector<length_d> const deltaPos = position - sphere.getCenter();
+
+ { // check extreme cases we don't want to solve analytically explicit
+ HEPMomentumType const p_perp =
+ (particle.getMomentum() -
+ particle.getMomentum().getParallelProjectionOnto(magneticfield))
+ .getNorm();
+
+ LengthType const gyroradius =
+ (convert_HEP_to_SI<MassType::dimension_type>(p_perp) * constants::c /
+ (abs(charge) * magneticfield.getNorm()));
+
+ LengthType const trackDist = abs(deltaPos.getNorm() - radius);
+ if (trackDist > gyroradius) {
+ // there is never a solution
+ return Intersections();
+ }
+
+ if (gyroradius > 1000 * trackDist) {
+ // the bending is negligible, use straight intersections instead
+ return tracking_line::Tracking::intersect(particle, sphere);
+ }
+ }
+
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 = charge / p_norm;
+ decltype(1 / (tesla * meter)) const k = charge / p_norm;
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_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 /
+ double const d = (deltaPosLength + radius) * (deltaPosLength - radius) * denom /
(1_m * 1_m * 1_m * 1_m);
CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
- std::vector<double> solutions = andre::solve_quartic_real(1, 0, b, c, d);
+ // solutions of deltaL are obtained from quartic equation. Note, deltaL/2 is the
+ // length of each half step, however, the second half step is slightly longer
+ // because of the non-conservation of norm/velocity.
+ // The leap-frog length L is deltaL/2 * (1+|u_{n+1}|)
+ std::vector<double> solutions = solve_quartic_real(1, 0, b, c, d);
+ if (!solutions.size()) { return Intersections(); }
LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0;
for (auto solution : solutions) {
LengthType const dist = solution * 1_m;
- CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
+ CORSIKA_LOG_TRACE(
+ "Solution (real) for current Volume: deltaL/2*2={} (deltaL/2*2/v={}) ", dist,
+ dist / absVelocity);
if (numericallyInside) {
// there must be an entry (negative) and exit (positive) solution
if (dist < 0.0001_m) { // security margin to assure
@@ -237,6 +287,8 @@ namespace corsika {
first_entry, first_exit);
return Intersections();
}
+ // return in units of time
+
return Intersections(d_enter / absVelocity, d_exit / absVelocity);
}
@@ -248,7 +300,7 @@ namespace corsika {
ElectricChargeType const charge = particle.getCharge();
- if (charge != 0 * constants::e) {
+ if (charge != ElectricChargeType::zero()) {
auto const* currentLogicalVolumeNode = particle.getNode();
VelocityVector const velocity = particle.getVelocity();
@@ -279,6 +331,10 @@ namespace corsika {
plane.getNormal().dot(position - plane.getCenter())) // unit: kg*m/s *m
/ (1_m * 1_m * 1_kg) * 1_s;
+ // deltaL from quadratic solution return half-step length deltaL/2 for leap-frog
+ // algorithmus. Note, the leap-frog length L is longer by (1+|u_{n_1}|)/2 because
+ // the direction norm of the second half step is >1.
+
std::vector<double> const deltaLs = solve_quadratic_real(denom, p, q);
CORSIKA_LOG_TRACE("deltaLs=[{}]", fmt::join(deltaLs, ", "));
@@ -305,15 +361,11 @@ namespace corsika {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
}
- CORSIKA_LOG_TRACE("maxStepLength={}", maxStepLength);
+ CORSIKA_LOG_TRACE("maxStepLength={} s", maxStepLength / 1_s);
- // with final length correction
- auto const corr =
- (direction + direction_x_B * maxStepLength * charge / (p_norm * 2))
- .getNorm() /
- absVelocity; // unit: s/m
+ // with final length correction, |direction| becomes >1 during step
- return Intersections(maxStepLength * corr); // unit: s
+ return Intersections(maxStepLength / absVelocity); // unit: s
} // no charge
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
index 4e7efed775468a6525a236fe7b5a72b45896e518..675636a173490fa57c064ea3c7d893c2d4207f42 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
@@ -63,7 +63,7 @@ namespace corsika {
// check, where the first halve-step direction has geometric intersections
auto const [initialTrack, initialTrackNextVolume] =
tracking_line::Tracking::getTrack(particle);
- { [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
+ //{ [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
auto const initialTrackLength = initialTrack.getLength(1);
CORSIKA_LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
@@ -82,8 +82,8 @@ namespace corsika {
.getNorm();
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.
+ // particle travels along, parallel to magnetic field. Rg is
+ // "0", return straight track here.
CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
@@ -105,8 +105,9 @@ namespace corsika {
DirectionVector const direction = initialVelocity.normalized();
// avoid any intersections within first halve steplength
+ // aim for intersection in second halve step
LengthType const firstHalveSteplength =
- std::min(steplimit, initialTrackLength * firstFraction_);
+ std::min(steplimit / 2, initialTrackLength * firstFraction_);
CORSIKA_LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
firstHalveSteplength, steplimit, initialTrackLength);
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index cad8ac39d99ebf67af35c19eea7974897806b3f9..3bdbb73df381967627ac9a596a1188aac8252490 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -24,6 +24,15 @@
namespace corsika::tracking_line {
+ template <typename TParticle>
+ inline auto Tracking::makeStep(TParticle const& particle, LengthType steplength) {
+ if (particle.getMomentum().getNorm() == 0_GeV) {
+ return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV);
+ } // charge of the particle
+ DirectionVector const dir = particle.getDirection();
+ return std::make_tuple(particle.getPosition() + dir * steplength, dir.normalized());
+ }
+
template <typename TParticle>
inline auto Tracking::getTrack(TParticle const& particle) {
VelocityVector const initialVelocity =
@@ -92,9 +101,10 @@ namespace corsika::tracking_line {
CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
particle.getMomentum());
- return Intersections(n_dot_v.magnitude() == 0
- ? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
- : n.dot(delta) / n_dot_v);
+ if (n_dot_v.magnitude() == 0)
+ return Intersections();
+ else
+ return Intersections(n.dot(delta) / n_dot_v);
}
} // namespace corsika::tracking_line
diff --git a/corsika/framework/geometry/Intersections.hpp b/corsika/framework/geometry/Intersections.hpp
index c9951213ef162fb273626f77dd8e3707f17403c2..ce0131ee39eee304190f9dfed7413b0502c8baf4 100644
--- a/corsika/framework/geometry/Intersections.hpp
+++ b/corsika/framework/geometry/Intersections.hpp
@@ -43,10 +43,20 @@ namespace corsika {
bool hasIntersections() const { return has_intersections_; }
///! where did the trajectory currently enter the volume
- TimeType getEntry() const { return intersections_.first; }
+ TimeType getEntry() const {
+ if (has_intersections_)
+ return intersections_.first;
+ else
+ return std::numeric_limits<TimeType::value_type>::infinity() * second;
+ }
///! where did the trajectory currently exit the volume
- TimeType getExit() const { return intersections_.second; }
+ TimeType getExit() const {
+ if (has_intersections_)
+ return intersections_.second;
+ else
+ return std::numeric_limits<TimeType::value_type>::infinity() * second;
+ }
private:
bool has_intersections_;
diff --git a/corsika/framework/geometry/LeapFrogTrajectory.hpp b/corsika/framework/geometry/LeapFrogTrajectory.hpp
index ec281e5f71d56c77602b2dd64b22d5e3e47ea8c9..ddccc83c80165713130e216b7118e984b5875281 100644
--- a/corsika/framework/geometry/LeapFrogTrajectory.hpp
+++ b/corsika/framework/geometry/LeapFrogTrajectory.hpp
@@ -19,7 +19,7 @@ namespace corsika {
/**
* The LeapFrogTrajectory stores information on one leap-frog step.
*
- * The leap-frog algorithm uses a half-step and is used in magnetic
+ * The leap-frog algorithm uses two half-steps and is used in magnetic
* field tracking. The LeapFrogTrajectory will solve the leap-frog
* algorithm equation for a given constant $k$ that has to be
* specified during construction (essentially fixing the magnetic
@@ -28,6 +28,17 @@ namespace corsika {
* direction as calculated by the algorithm for any steplength, or
* intermediate position.
*
+ * One complete leap-frog step is
+ * $
+ * \vec{x}(t_{i+0.5}) = \vec{x}(t_{i}) + \vec{v(t_{i})} * \Delta t / 2 \\
+ * \vec{v}(t_{i+1}) = \vec{v}(t_{i}) + \vec{v}(t_{i})\cross\vec{B}(x_{i}, t_{i}) *
+ * \Delta t \\
+ * \vec{x}(t_{i+1}) = \vec{x}(t_{i+0.5}) + \vec{v}(t_{i+1}) * \Delta t /2 \\
+ * $
+ *
+ * The volocity update has the characteristics $|\vec{v}(t_{i+1})|>1$, thus final
+ * velocities are renormalised. The full leap-frog steplength is thus
+ * $L = |\vec{v}(t_{i+1})| * \Delta t / 2 + |\vec{v}(t_{i+1})| *\Delta t / 2$
**/
class LeapFrogTrajectory : public BaseTrajectory {
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
index 513ffbaf509d2c05a7020f038d0ab320592c4712..0dc06ca410ee8d08972f2aa6fedbce51f7f6884a 100644
--- a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -29,24 +29,21 @@ namespace corsika {
namespace tracking_leapfrog_curved {
/**
- * \file TrackingLeapFrogCurved.hpp
- *
- * Performs one leap-frog step consistent of two halve-steps with steplength/2
- * The step is caluculated analytically precisely to reach to the next volume
- *boundary.
- **/
- template <typename TParticle>
- auto make_LeapFrogStep(TParticle const& particle, LengthType steplength);
+ * \file TrackingLeapFrogCurved.hpp The leap-frog tracking.
+ */
/**
- *
* The class tracking_leapfrog_curved::Tracking is based on the
* Bachelor thesis of Andre Schmidt (KIT). It implements a
* two-step leap-frog algorithm, but with analytically exact geometric
* intersections between leap-frog steps and geometric volumes
* (spheres, planes).
*
- **/
+ * Note that leap-frog times and length always reflect the actual properties of
+ * the final step. The internal steplength is slightly shorter, because the second
+ * halve steps of the algorithm is slightly longer than the first one (in principle
+ * violating |v|=const).
+ */
class Tracking : public Intersect<Tracking> {
@@ -59,16 +56,49 @@ namespace corsika {
template <typename TParticle>
auto getTrack(TParticle const& particle);
- //! find intersection of Sphere with Track
+ /**
+ * Performs one leap-frog step consistent of two halve-steps with steplength/2
+ * Due to the nature of the algorithm the second halve step is slightly longer than
+ * the first halve step.
+ */
+ template <typename TParticle>
+ static auto makeStep(TParticle const& particle, LengthType const steplength);
+
+ /**
+ * find intersection of Sphere with Track
+ *
+ * Returns intersection of particle assuming a curved leap-frog step, with a sphere.
+ * Entry and exit times are calculated, where the velocity is constant and the
+ * steplength is the geometric steplength of the leap-frog.
+ *
+ * @param particle Current particle state
+ * @param sphere Sphere object
+ */
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Sphere const& sphere);
- //! find intersection of Volume node with Track of particle
+ /**
+ * find intersection of any Volume node with particle
+ *
+ * The intersection time(s) of a particle, assuming a curved leap-frog
+ * step, are calculated for any volume type.
+ */
template <typename TParticle, typename TBaseNodeType>
static Intersections intersect(TParticle const& particle,
TBaseNodeType const& node);
- //! find intersection of Plane with Track
+ /**
+ * find intersection of Plane with Track
+ *
+ * Intersection times of particle are caculated with a plane, assuming a curved leap
+ * frog trajectory. The intersection time is assuming constant velocity (no change)
+ * along the geometric leap-frog step.
+ *
+ * @tparam TParticle Type of particle object on stack.
+ * @param particle Particle initial state.
+ * @param plane Plane.
+ * @return Intersections in time units.
+ */
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Plane const& plane);
@@ -80,7 +110,6 @@ namespace corsika {
* Use internally stored class tracking_line::Tracking to
* perform a straight line tracking, if no magnetic bendig was
* detected.
- *
*/
template <typename TParticle>
auto getLinearTrajectory(TParticle& particle);
diff --git a/corsika/modules/tracking/TrackingStraight.hpp b/corsika/modules/tracking/TrackingStraight.hpp
index a5800cd7ed64fa4fee8afe89a32fe927366ad1de..3ebc18b37429eb78733fc3886ec8b07e11f648e5 100644
--- a/corsika/modules/tracking/TrackingStraight.hpp
+++ b/corsika/modules/tracking/TrackingStraight.hpp
@@ -36,6 +36,12 @@ namespace corsika::tracking_line {
using Intersect<Tracking>::nextIntersect;
public:
+ /**
+ * Performs one straight step of length steplength.
+ */
+ template <typename TParticle>
+ static auto makeStep(TParticle const& particle, LengthType const steplength);
+
//! Determine track of particle
template <typename TParticle>
auto getTrack(TParticle const& particle);
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 7aa0ad850d41aad0cf6e19a84422f0b74845cf81..e3d056e22202e7af786c3f1eec00fcc36e077a61 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -59,3 +59,5 @@ CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
add_executable (corsika corsika.cpp)
target_link_libraries (corsika CORSIKA8::CORSIKA8)
+
+
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 4f43447844e848359e24a73bb3c75e7a0aca445f..e065c4ee8ddc787e2bd56d2f8502f37cfb7d0bd1 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -67,7 +67,7 @@
using namespace corsika;
using namespace std;
-void registerRandomStreams(int seed) {
+void registerRandomStreams(uint64_t seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("qgsjet");
RNGManager<>::getInstance().registerRandomStream("sibyll");
@@ -82,6 +82,34 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().setSeed(seed);
}
+class TrackCheck : public ContinuousProcess<TrackCheck> {
+
+public:
+ TrackCheck(Plane const& plane)
+ : plane_(plane) {}
+
+ template <typename TParticle, typename TTrack>
+ ProcessReturn doContinuous(TParticle const& particle, TTrack const&, bool const) {
+ auto const delta = particle.getPosition() - plane_.getCenter();
+ auto const n = plane_.getNormal();
+ auto const proj = n.dot(delta);
+ if (proj < -1_m) {
+ CORSIKA_LOG_INFO("particle {} failes: proj={}, delta={}, p={}", particle.asString(),
+ proj, delta, particle.getPosition());
+ throw std::runtime_error("particle below obs level");
+ }
+ return ProcessReturn::Ok;
+ }
+
+ template <typename TParticle, typename TTrack>
+ LengthType getMaxStepLength(TParticle const&, TTrack const&) const {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+
+private:
+ Plane plane_;
+};
+
template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
@@ -98,8 +126,8 @@ int main(int argc, char** argv) {
}
feenableexcept(FE_INVALID);
- int seed = 0;
- if (argc > 4) seed = std::stoi(std::string(argv[4]));
+ uint64_t seed = 0;
+ if (argc > 4) seed = std::stol(std::string(argv[4]));
// initialize random number sequence(s)
registerRandomStreams(seed);
@@ -209,7 +237,7 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
- ParticleCut cut{60_GeV, false, true};
+ ParticleCut cut{3_GeV, false, true};
BetheBlochPDG eLoss{showerAxis};
CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0,
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index 1335974cdd70104ced3ff639b8afb3610fffcec3..cf9dc523f4a3fe571e3d3323b957534ac790b285 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -369,7 +369,6 @@ private:
TEST_CASE("ProcessSequence General", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
SECTION("BaseProcess") {
@@ -889,7 +888,6 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
TEST_CASE("ProcessSequence Indexing", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
SECTION("Indexing") {
diff --git a/tests/framework/testSolver.cpp b/tests/framework/testSolver.cpp
index 089335ed5ee6d42d25ec8ddee4a47e459c503325..197061b813ce2b7fdee427e9a450d6291f1a7785 100644
--- a/tests/framework/testSolver.cpp
+++ b/tests/framework/testSolver.cpp
@@ -98,35 +98,64 @@ TEST_CASE("Solver") {
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, ", "));
+ {
+ 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))));
+ }
+ }
+ }
- 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))));
+ // cubic with x^3 * 0 should result in the same
+ {
+ std::vector<double> s1 = solve_cubic_real(0, a, b, c);
+ remove_duplicates(s1, epsilon_check);
+
+ CORSIKA_LOG_INFO("quadratic/cubic: 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))));
+
+ // quartic with x^4 * 0 + x^3 * 0 must be the same
+ {
+ std::vector<double> s1 = solve_quartic_real(0, 0, a, b, c);
+ remove_duplicates(s1, epsilon_check);
+
+ CORSIKA_LOG_INFO("quadratic/quartic: 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))));
+ }
}
}
- }*/
}
}
}
@@ -167,23 +196,48 @@ TEST_CASE("Solver") {
"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, ", "));
+ {
+ vector<double> s1 = solve_cubic_real(a, b, c, d);
+ remove_duplicates(s1, epsilon_check * 10);
+
+ 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))));
+ }
+ }
+ }
- 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))));
+ // quartic with x^4 *0 must be the same
+ {
+ vector<double> s1 = solve_quartic_real(0, a, b, c, d);
+ remove_duplicates(s1, epsilon_check * 10);
+
+ CORSIKA_LOG_INFO("(quartic) 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))));
+ }
}
}
}
@@ -193,7 +247,7 @@ TEST_CASE("Solver") {
SECTION("quartic") {
- epsilon_check = 1e-4; // for catch2 asserts
+ epsilon_check = 1e-2; // for catch2 asserts
// **clang-format-off**
// tests of type:
@@ -251,10 +305,12 @@ TEST_CASE("Solver") {
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);
+ remove_duplicates(s1, epsilon_check * 10);
+ vector<double> s2 = solve_quartic_real(a, b, c, d, e);
+ remove_duplicates(s2, epsilon_check * 5);
CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", "));
+ CORSIKA_LOG_INFO("N={}, s2=[{}]", s2.size(), fmt::join(s2, ", "));
CHECK(s1.size() == idegree + 1);
for (double value : s1) {
@@ -272,8 +328,27 @@ TEST_CASE("Solver") {
(value == Approx(z4).epsilon(epsilon_check))));
}
}
+
+ // this is a bit less precise
+ CHECK(s2.size() == idegree + 1);
+ for (double value : s2) {
+ 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 * 5)) ||
+ (value == Approx(z2).margin(epsilon_check * 5)) ||
+ (value == Approx(z3).margin(epsilon_check * 5)) ||
+ (value == Approx(z4).margin(epsilon_check * 5))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check * 5)) ||
+ (value == Approx(z2).epsilon(epsilon_check * 5)) ||
+ (value == Approx(z3).epsilon(epsilon_check * 5)) ||
+ (value == Approx(z4).epsilon(epsilon_check * 5))));
+ }
+ }
}
}
}
+
} // quartic
}
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index e6cd55e8f51fe9d96c1f92c508362e68de0b9969..f6809d4a61b4e6015e625a5c302ab8e7fd8ad713 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -118,10 +118,12 @@ TEST_CASE("Pythia8Interface", "modules") {
corsika::pythia8::Decay decay(particleList);
+ CORSIKA_LOG_INFO("stack: {} {}", stack.asString(), particle.asString());
[[maybe_unused]] const TimeType time = decay.getLifetime(particle);
double const gamma = particle.getEnergy() / get_mass(Code::PiPlus);
CHECK(time == get_lifetime(Code::PiPlus) * gamma);
decay.doDecay(*secViewPtr);
+ CORSIKA_LOG_INFO("piplus->{}", stack.asString());
CHECK(stack.getEntries() == 3); // piplus, muplu, numu
auto const pSum = sumMomentum(view, cs);
CHECK((pSum - plab).getNorm() / 1_GeV == Approx(0).margin(1e-4));
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index fb398e795668ecb87c230ed03a047fad17af9049..bbd94cb347e35adfec4b2eada8199258912b6d43 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -70,34 +70,37 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
const HEPEnergyType P0 = 10_GeV;
auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Photon);
+ // test also special case: movement parallel to field (along x)
+ auto isParallel = GENERATE(as<bool>{}, true, false);
// for algorithms that know magnetic deflections choose: +-50uT, 0uT
// otherwise just 0uT
- auto Bfield = GENERATE(filter(
- []([[maybe_unused]] MagneticFluxType v) {
+ auto Bfield = GENERATE_COPY(filter(
+ [isParallel]([[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})));
+ values<MagneticFluxType>({50_uT, 0_uT, (isParallel ? 0_uT : -50_uT)})));
// particle --> (world) --> | --> (target)
// true: start inside "world" volume
// false: start inside "target" volume
auto outer = GENERATE(as<bool>{}, true, false);
- SECTION(fmt::format("Tracking PID={}, Bfield={} uT, from outside={}", PID,
- Bfield / 1_uT, outer)) {
+ SECTION(fmt::format("Tracking PID={}, Bfield={} uT, isParallel={}, from outside={}",
+ PID, Bfield / 1_uT, isParallel, outer)) {
CORSIKA_LOG_DEBUG(
"********************\n TEST algo={} section PID={}, "
"Bfield={} "
- "uT, start_outside={}",
- boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT, outer);
+ "uT, field is parallel={}, start_outside={}",
+ boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT,
+ isParallel, outer);
const int chargeNumber = get_charge_number(PID);
LengthType radius = 10_m;
int deflect = 0;
- if (chargeNumber != 0 and Bfield != 0_T) {
+ if (chargeNumber != 0 && Bfield != 0_T && !isParallel) {
deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection
LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(P0) *
constants::c / (abs(get_charge(PID)) * abs(Bfield)));
@@ -113,13 +116,25 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_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
+ // it is very close to injection and not so small
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
+ // this is far from injection and really small
auto target_neutral = setup::Environment::createNode<Sphere>(
Point(cs, {radius / 2, 0_m, 0_m}), radius * 0.1);
+ // target to be overlapped entirely by target_2
+ auto target_2_behind = setup::Environment::createNode<Sphere>(
+ Point(cs, {-radius * 3 / 4, 0_m, 0_m}), radius * 0.1);
+
+ // target to be overlapped partly by target_2
+ auto target_2_partly_behind = setup::Environment::createNode<Sphere>(
+ Point(cs, {-radius * 3 / 4 + radius * 0.1, 0_m, 0_m}), radius * 0.2);
+
using MyHomogeneousModel = MediumPropertyModel<
UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
@@ -133,12 +148,24 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
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.}));
+ target_2_behind->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_partly_behind->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();
+ auto* targetPtr_2_behind = target_2_behind.get();
+ auto* targetPtr_2_partly_behind = target_2_partly_behind.get();
+ target_2_behind->excludeOverlapWith(target_2);
+ target_2_partly_behind->excludeOverlapWith(target_2);
worldPtr->addChild(std::move(target));
targetPtr->addChild(std::move(target_2));
targetPtr->addChild(std::move(target_neutral));
+ targetPtr->addChild(std::move(target_2_behind));
+ targetPtr->addChild(std::move(target_2_partly_behind));
auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs);
{ [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
@@ -167,9 +194,13 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
// move forward, until we leave target volume
bool hit_neutral = false;
bool hit_2nd = false;
+ bool hit_2nd_behind = false;
+ [[maybe_unused]] bool hit_2nd_partly_behind = false;
while (nextVol != worldPtr) {
if (nextVol == targetPtr_neutral) hit_neutral = true;
if (nextVol == targetPtr_2) hit_2nd = true;
+ if (nextVol == targetPtr_2_behind) hit_2nd_behind = true;
+ if (nextVol == targetPtr_2_partly_behind) hit_2nd_partly_behind = true;
const auto [traj2, nextVol2] = tracking.getTrack(particle);
nextVol = nextVol2;
particle.setNode(nextVol);
@@ -181,6 +212,10 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
particle.getVelocity().getNorm(), traj2.getLength(1),
traj2.getLength(1) / particle.getVelocity().getNorm());
}
+ CHECK_FALSE(hit_2nd_behind); // this can never happen
+ // the next line is maybe an actual BUG: this should be investigated and eventually
+ // fixed:
+ // CHECK(hit_2nd == hit_2nd_partly_behind); // if one is hit, the other also must be
CHECK(nextVol == worldPtr);
CHECK(hit_2nd == true);
CHECK(hit_neutral == (deflect == 0 ? true : false));