diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 90df44a82e23d15fac2f22cdb233f943ac552bc2..9ba4c4bc87abb84f2a025872a8b095c0f4719fde 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -227,12 +227,11 @@ namespace corsika::cascade { std::cout << "distance_max=" << distance_max << std::endl; // take minimum of geometry, interaction, decay for next step - std::cout << "Interaction: " << distance_interact << std::endl; - std::cout << "Decay: " << distance_decay << std::endl; - std::cout << "ObsPlane: " << distance_max << std::endl; - std::cout << "Transition: " << geomMaxLength << std::endl; - auto const min_distance = std::min( + auto min_distance = std::min( {distance_interact, distance_decay, distance_max, geomMaxLength}); + if (min_distance == geomMaxLength) { + min_distance = 1.001 * geomMaxLength; + } C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m); @@ -260,11 +259,11 @@ namespace corsika::cascade { // Second Movement position = position + directionAfter * min_distance / 2; auto distance = position - vParticle.GetPosition(); - //distance.norm() != min_distance if q != 0 - //small error can be neglected + // distance.norm() != min_distance if q != 0 + // small error can be neglected if (distance.norm() != 0_m) { velocity = distance.normalized() * velocity.norm(); - } //no velocity update for very small steps + } // no velocity update for very small steps // here the particle is actually moved along the trajectory to new position: // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step); diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc index 298b53b8815f84873093c433c6d803f8e8ce52dc..d12a55f5ad3c5d808ae0c30f6a812276f2e7e777 100644 --- a/Processes/ObservationPlane/ObservationPlane.cc +++ b/Processes/ObservationPlane/ObservationPlane.cc @@ -88,6 +88,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) - velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) / (plane_.GetNormal().dot(velocity.cross(magneticfield)) * k); + if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) { return std::numeric_limits<double>::infinity() * 1_m; } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) { diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h index fd89415abac9318d5979e596f117e8915c550d28..eff54e7e2e126700b386a9195ce10b5e31c3ba8c 100644 --- a/Processes/TrackingLine/TrackingLine.h +++ b/Processes/TrackingLine/TrackingLine.h @@ -135,7 +135,7 @@ namespace corsika::process { (position - currentPosition).GetNorm(); velocity1 = direction * velocity.GetNorm(); } // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed - // line has some errors + // using line has some errors for huge steps geometry::Line line(currentPosition, velocity1); if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) { @@ -232,9 +232,10 @@ namespace corsika::process { min * k; // Second Movement position = position + directionAfter * velocity.norm() * min / 2; - geometry::Vector<dimensionless_d> const direction = (position - currentPosition) / - (position - currentPosition).GetNorm(); - velocity = direction * velocity.GetNorm(); + if ((position - currentPosition).GetNorm() != 0_m) { + geometry::Vector<dimensionless_d> const direction = (position - currentPosition).normalized(); + velocity = direction * velocity.norm(); + } // no velocity update for very small steps geometry::Line lineWithB(currentPosition, velocity); return std::make_tuple(geometry::Trajectory<geometry::Line>(lineWithoutB, min),