check for particles below obs level in hybridMC

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);

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) {

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) {

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;