From f79ed358ce81deee0e361e0563d1177bd5d06d33 Mon Sep 17 00:00:00 2001
From: Fan <fan_hu@pku.edu.cn>
Date: Sun, 28 Nov 2021 15:39:10 +0800
Subject: [PATCH] fix logical bugs
---
corsika/detail/modules/ObservationCubic.inl | 6 +++---
corsika/detail/modules/tracking/TrackingStraight.inl | 12 +++++++++---
corsika/framework/geometry/Cubic.hpp | 4 ++--
3 files changed, 14 insertions(+), 8 deletions(-)
diff --git a/corsika/detail/modules/ObservationCubic.inl b/corsika/detail/modules/ObservationCubic.inl
index 6e7428db2..1967e3dd6 100644
--- a/corsika/detail/modules/ObservationCubic.inl
+++ b/corsika/detail/modules/ObservationCubic.inl
@@ -18,7 +18,7 @@ ObservationCubic<TTracking, TOutput>::ObservationCubic(
template <typename TTracking, typename TOutput>
template <typename TParticle, typename TTrajectory>
inline ProcessReturn ObservationCubic<TTracking, TOutput>::doContinuous(
- TParticle &particle, TTrajectory &, bool const stepLimit) {
+ TParticle &particle, TTrajectory & step, bool const stepLimit) {
/*
The current step did not yet reach the ObservationCubic, do nothing now and
wait:
@@ -42,7 +42,7 @@ inline ProcessReturn ObservationCubic<TTracking, TOutput>::doContinuous(
}
HEPEnergyType const energy = particle.getEnergy();
- Point const pointOfIntersection = particle.getPosition();
+ Point const pointOfIntersection = step.getPosition(1);
DirectionVector const dirction = particle.getDirection();
// add our particles to the output file stream
@@ -69,7 +69,7 @@ inline LengthType ObservationCubic<TTracking, TOutput>::getMaxStepLength(
CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, cubic={}",
particle.asString(), particle.getPosition(),
- particle.getDirection(), cubic_.asString());
+ particle.getDirection(), asString());
auto const intersection =
TTracking::intersect(particle, static_cast<Cubic const>(*this));
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index a88aed654..545ee6c8c 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -100,10 +100,14 @@ inline Intersections Tracking::intersect(TParticle const &particle,
SpeedType vx = velocity.getX(cs);
SpeedType vy = velocity.getY(cs);
SpeedType vz = velocity.getZ(cs);
+ CORSIKA_LOG_TRACE("particle in cubic coordinate: position: ({:.3f}, {:.3f}, "
+ "{:.3f}) m, veolocity: ({:.3f}, {:.3f}, {:.3f}) m/ns",
+ x0 / 1_m, y0 / 1_m, z0 / 1_m, vx / (1_m / 1_ns),
+ vy / (1_m / 1_ns), vz / (1_m / 1_ns));
auto get_intersect_min_max = [](LengthType x0, SpeedType v0, LengthType dx) {
- auto t1 = (x0 - dx) / v0;
- auto t2 = (x0 + dx) / v0;
+ auto t1 = (dx - x0) / v0;
+ auto t2 = (-dx - x0) / v0;
if (t1 > t2)
return std::make_pair(t1, t2);
else
@@ -115,8 +119,10 @@ inline Intersections Tracking::intersect(TParticle const &particle,
auto [tz_max, tz_min] = get_intersect_min_max(z0, vz, cubic.getZ());
TimeType t_exit = std::min(std::min(tx_max, ty_max), tz_max);
- TimeType t_enter = std::min(std::min(tx_min, ty_min), tz_min);
+ TimeType t_enter = std::max(std::max(tx_min, ty_min), tz_min);
+ CORSIKA_LOG_DEBUG("t_enter: {} ns, t_exit: {} ns", t_enter / 1_ns,
+ t_exit / 1_ns);
if ((t_exit > t_enter)) {
if (t_enter < 0_s && t_exit > 0_s)
CORSIKA_LOG_DEBUG("numericallyInside={}", cubic.contains(position));
diff --git a/corsika/framework/geometry/Cubic.hpp b/corsika/framework/geometry/Cubic.hpp
index 114965432..77b49a00e 100644
--- a/corsika/framework/geometry/Cubic.hpp
+++ b/corsika/framework/geometry/Cubic.hpp
@@ -12,12 +12,12 @@ namespace corsika {
// a CoordinateSystemPtr to specify the orintation of coordinate
Cubic(Point const& center, CoordinateSystemPtr cs, LengthType const x, LengthType const y, LengthType const z)
: center_(center)
- , cs_(make_translation(cs, center.getCoordinates(center.getCoordinateSystem())))
+ , cs_(make_translation(cs, center.getCoordinates(cs)))
, x_(x), y_(y), z_(z) {}
Cubic(Point const& center, CoordinateSystemPtr cs, LengthType const side)
: center_(center)
- , cs_(make_translation(cs, center.getCoordinates(center.getCoordinateSystem())))
+ , cs_(make_translation(cs, center.getCoordinates(cs)))
, x_(side/2), y_(side/2), z_(side/2) {}
--
GitLab