diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
index c15594e02a02bef2683e4d3ae4d4ad02895ef747..edfce000082c34a6645c459325f8b45f0a6fcb4e 100644
--- a/corsika/detail/framework/utility/CubicSolver.inl
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -234,10 +234,21 @@ namespace corsika {
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));
+#ifdef DEBUG
+ {
+ auto test = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
+
+ for (long double test_v : test) {
+ CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
+ cubic_function(test_v, a, b, c, d));
+ }
+ }
+#endif
+
if (std::abs(a) < epsilon) { // this is just a quadratic
return solve_quadratic_real(b, c, d, epsilon);
}
@@ -276,6 +287,10 @@ namespace corsika {
} while ((++niter < maxiter) && (std::abs(f_x1) > epsilon));
CORSIKA_LOG_TRACE("niter={}", niter);
+ if (niter >= maxiter) {
+ CORSIKA_LOG_TRACE("failure, no solution");
+ // return std::vector<double>{};
+ }
}
CORSIKA_LOG_TRACE("x1={} f_x1={}", x1, f_x1);
diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
index 4562b39c0aa0fea4330fc686deab57d48683d98f..764bc3a85e0e127e595801418cb2809ae5d617c1 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
+ }
long double y = x3[0]; // there is always at least one solution
// The essence - choosing Y with maximal absolute value.
if (x3.size() == 3) {
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 1ede8cc9825bb926c6423511dcc701005cc26a40..749693bc7c7a327732ea30d77872006f63eff5da 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -101,8 +101,8 @@ 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
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
@@ -150,6 +150,8 @@ namespace corsika {
auto const projectedDirectionSqrNorm = projectedDirection.getSquaredNorm();
bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
+ CORSIKA_LOG_TRACE("projectedDirectionSqrNorm={} T^2",
+ projectedDirectionSqrNorm / square(1_T));
if ((charge == 0 * constants::e) || magneticfield.getNorm() == 0_T || isParallel) {
return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
}
@@ -174,7 +176,10 @@ namespace corsika {
(deltaPosLength - sphere.getRadius()) * 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);
+ std::vector<double> solutions = solve_quartic_real(1, 0, b, c, d);
+ if (!solutions.size()) {
+ return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
+ }
LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0;
for (auto solution : solutions) {
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 159baed4f73255038b226236a674f0c916500570..abafdf0c3ab305d008c0ef6d13cb3885e6d08810 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -90,11 +90,12 @@ public:
template <typename TParticle, typename TTrack>
ProcessReturn doContinuous(TParticle const& particle, TTrack const&, bool const) {
- auto const delta = plane_.getCenter() - particle.getPosition();
+ auto const delta = particle.getPosition() - plane_.getCenter();
auto const n = plane_.getNormal();
auto const proj = n.dot(delta);
- if (proj < -0_mm) {
- CORSIKA_LOG_INFO("particle {} failes", particle.asString());
+ 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;
@@ -115,7 +116,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::trace);
CORSIKA_LOG_INFO("hybrid_MC");
@@ -160,7 +161,7 @@ int main(int argc, char** argv) {
unsigned short Z = std::stoi(std::string(argv[2]));
auto const mass = get_nucleus_mass(A, Z);
const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
- double theta = 0.;
+ double theta = 50.;
auto const thetaRad = theta / 180. * M_PI;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
@@ -268,8 +269,8 @@ int main(int argc, char** argv) {
make_sequence(sibyllNucCounted, sibyllCounted));
auto decaySequence = make_sequence(decayPythia, decaySibyll);
auto sequence =
- make_sequence(TrackCheck(obsPlane), hadronSequence, reset_particle_mass,
- decaySequence, eLoss, cut, conex_model, longprof, observationLevel);
+ make_sequence(hadronSequence, reset_particle_mass, decaySequence, eLoss, cut,
+ conex_model, longprof, observationLevel, TrackCheck(obsPlane));
// define air shower object, run simulation
setup::Tracking tracking;