From ade578778bdaa8f9bfe30af4341fc9aeda88f5f0 Mon Sep 17 00:00:00 2001
From: Andre Schmidt <schmidt-a@iklx288.ikp.kit.edu>
Date: Sun, 16 Aug 2020 15:01:41 +0200
Subject: [PATCH] mal 1.0001
---
Framework/Cascade/Cascade.h | 15 +++++++--------
Processes/ObservationPlane/ObservationPlane.cc | 1 +
Processes/TrackingLine/TrackingLine.h | 9 +++++----
3 files changed, 13 insertions(+), 12 deletions(-)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 90df44a82..9ba4c4bc8 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 298b53b88..d12a55f5a 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 fd89415ab..eff54e7e2 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),
--
GitLab