From 1ed8143063e381849421fbc4aef470607dfc6635 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Fri, 20 Nov 2020 09:19:10 +0100
Subject: [PATCH] fixed some remaining tracking issues, mostly with infinite
tracks
---
Documentation/Examples/boundary_example.cc | 2 +-
.../Examples/cascade_proton_example.cc | 2 +-
Framework/Cascade/Cascade.h | 9 +--
Framework/Cascade/testCascade.cc | 11 +++-
Framework/Geometry/Trajectory.h | 56 ++++++++++++++++---
5 files changed, 62 insertions(+), 18 deletions(-)
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 010cb8db2..21cdbe353 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -85,7 +85,7 @@ private:
//
int main() {
- logging::SetLevel(logging::level::info);
+ logging::SetLevel(logging::level::trace);
C8LOG_INFO("boundary_example");
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index 6946f2490..11b3dac77 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -60,7 +60,7 @@ using namespace corsika::units::si;
//
int main() {
- logging::SetLevel(logging::level::info);
+ logging::SetLevel(logging::level::trace);
std::cout << "cascade_proton_example" << std::endl;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 6dc602de2..e1366cbcc 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -213,9 +213,10 @@ namespace corsika::cascade {
"transport particle by : {} m "
"Medium transition after: {} m "
"Decay after: {} m "
- "Interaction after: {} m",
+ "Interaction after: {} m"
+ "Continuous limit: {} m",
min_distance / 1_m, geomMaxLength / 1_m, distance_decay / 1_m,
- distance_interact / 1_m);
+ distance_interact / 1_m, distance_max / 1_m);
// here the particle is actually moved along the trajectory to new position:
step.SetLength(min_distance);
@@ -236,9 +237,9 @@ namespace corsika::cascade {
}
C8LOG_DEBUG("sth. happening before geometric limit ? {}",
- ((min_distance <= geomMaxLength) ? "yes" : "no"));
+ ((min_distance < geomMaxLength) ? "yes" : "no"));
- if (min_distance <= geomMaxLength) { // interaction to happen within geometric limit
+ if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
// check whether decay or interaction limits this step the
// outcome of decay or interaction MAY be a) new particles in
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 2ca4dd7f8..a4863d65e 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -75,9 +75,14 @@ public:
geometry::Vector<SpeedType::dimension_type> const initialVelocity =
particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
return std::make_tuple(
- geometry::LineTrajectory(geometry::Line(particle.GetPosition(), initialVelocity),
- 0_s), // trajectory
- particle.GetNode()); // next volume node
+ geometry::LineTrajectory(
+ geometry::Line(particle.GetPosition(), initialVelocity),
+ std::numeric_limits<TimeType::value_type>::infinity() * 1_s), // trajectory,
+ // just
+ // go
+ // ahead
+ // forever
+ particle.GetNode()); // next volume node
}
};
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 4f1a84242..013c92c56 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -44,12 +44,26 @@ namespace corsika::geometry {
LineTrajectory(const LineTrajectory&) = default;
LineTrajectory(LineTrajectory&&) = default;
LineTrajectory& operator=(const LineTrajectory&) = delete;
+
+ /**
+ * \param theLine The geometric \sa Line object that represents a straight-line
+ * connection \param timeLength The time duration to traverse the straight trajectory
+ * in units of \sa TimeType
+ */
LineTrajectory(Line const& theLine, corsika::units::si::TimeType timeLength)
: line_(theLine)
, timeLength_(timeLength)
, timeStep_(timeLength)
- , initialVelocity_(theLine.GetVelocity(timeLength * 0))
- , finalVelocity_(theLine.GetVelocity(timeLength * 1)) {}
+ , initialVelocity_(theLine.GetVelocity(corsika::units::si::TimeType::zero()))
+ , finalVelocity_(theLine.GetVelocity(timeLength)) {}
+
+ /**
+ * \param theLine The geometric \sa Line object that represents a straight-line
+ * connection \param timeLength The time duration to traverse the straight trajectory
+ * in units of \sa TimeType \param timeStep Time duration to folow eventually curved
+ * trajectory in units of \sa TimesType \param initialV Initial velocity vector at
+ * start of trajectory \param finalV Final velocity vector at start of trajectory
+ */
LineTrajectory(
Line const& theLine,
corsika::units::si::TimeType timeLength, // length of theLine (straight)
@@ -77,6 +91,11 @@ namespace corsika::geometry {
corsika::units::si::LengthType GetLength(double u = 1) const {
using namespace corsika::units::si;
if (timeLength_ == 0_s) return 0_m;
+ if (timeStep_ ==
+ std::numeric_limits<corsika::units::si::TimeType::value_type>::infinity() * 1_s)
+ return std::numeric_limits<
+ corsika::units::si::LengthType::value_type>::infinity() *
+ 1_m;
return GetDistance(u) * timeStep_ / timeLength_;
}
@@ -90,14 +109,32 @@ namespace corsika::geometry {
void SetDuration(corsika::units::si::TimeType limit) {
using namespace corsika::units::si;
if (timeStep_ == 0_s) {
- timeLength_ *= 0;
+ timeLength_ = 0_s;
SetFinalVelocity(GetVelocity(0));
timeStep_ = limit;
} else {
- const double scale = limit / timeStep_;
- timeLength_ *= scale;
- SetFinalVelocity(GetVelocity(scale));
- timeStep_ = limit;
+ // for infinite steps there can't be a difference between
+ // curved and straight trajectory, this is fundamentally
+ // undefined: assume they are the same (which, i.e. is always correct for a
+ // straight line trajectory).
+ //
+ // Final note: only straight-line trajectories should have
+ // infinite steps! Everything else is ill-defined.
+ if (timeStep_ == std::numeric_limits<
+ corsika::units::si::TimeType::value_type>::infinity() *
+ 1_s ||
+ timeLength_ == std::numeric_limits<
+ corsika::units::si::TimeType::value_type>::infinity() *
+ 1_s) {
+ timeLength_ = limit;
+ timeStep_ = limit;
+ // ...and don't touch velocity
+ } else {
+ const double scale = limit / timeStep_;
+ timeLength_ *= scale;
+ SetFinalVelocity(GetVelocity(scale));
+ timeStep_ = limit;
+ }
}
}
@@ -113,8 +150,9 @@ namespace corsika::geometry {
private:
Line line_;
- corsika::units::si::TimeType timeLength_;
- corsika::units::si::TimeType timeStep_;
+ corsika::units::si::TimeType
+ timeLength_; ///! length of straight step (shortest connecting line)
+ corsika::units::si::TimeType timeStep_; ///! length of bend step (curved)
VelocityVec initialVelocity_;
VelocityVec finalVelocity_;
};
--
GitLab