From cf2d0c5f2ec4d71069b10b20d62ff3ba528cf779 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 2 Feb 2021 11:23:29 +0100
Subject: [PATCH] fixed missing curved-plane intersection code and a few typos
---
corsika/detail/framework/core/Cascade.inl | 4 -
.../tracking/TrackingLeapFrogCurved.inl | 77 ++++++++++++++++++-
.../tracking/TrackingLeapFrogStraight.inl | 2 -
.../tracking/TrackingLeapFrogCurved.hpp | 12 ++-
corsika/setup/SetupTrajectory.hpp | 6 +-
tests/common/SetupTestTrajectory.hpp | 36 ++++-----
tests/framework/testCascade.cpp | 13 ++--
tests/media/testMedium.cpp | 2 +-
8 files changed, 112 insertions(+), 40 deletions(-)
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 0e9161259..75c74d2b7 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -18,9 +18,6 @@
#include <corsika/framework/stack/SecondaryView.hpp>
#include <corsika/media/Environment.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
#include <cassert>
#include <cmath>
#include <iostream>
@@ -139,7 +136,6 @@ namespace corsika {
distance_interact / 1_m, continuous_max_dist / 1_m);
// here the particle is actually moved along the trajectory to new position:
- // std::visit(setup::ParticleUpdate<particle_type>{vParticle}, step);
step.setLength(min_distance);
vParticle.setPosition(step.getPosition(1));
// assumption: tracking does not change absolute momentum:
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 8545d2bbd..1691b0486 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -30,7 +30,7 @@ namespace corsika {
template <typename TParticle>
inline auto make_LeapFrogStep(TParticle const& particle, LengthType steplength) {
- if (particle.getMomentum().norm() == 0_GeV) {
+ if (particle.getMomentum().getNorm() == 0_GeV) {
return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
double(0));
} // charge of the particle
@@ -225,6 +225,81 @@ namespace corsika {
return Intersections(d_enter / absVelocity, d_exit / absVelocity);
}
+
+ template <typename TParticle, typename TMedium>
+ Intersections Tracking::intersect(TParticle const& particle, Plane const& plane,
+ TMedium const& medium) {
+
+ int chargeNumber;
+ if (is_nucleus(particle.getPID())) {
+ chargeNumber = particle.getNuclearZ();
+ } else {
+ chargeNumber = get_charge_number(particle.getPID());
+ }
+ auto const* currentLogicalVolumeNode = particle.getNode();
+ auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
+ particle.getPosition());
+ VelocityVector const velocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+ auto const absVelocity = velocity.getNorm();
+ DirectionVector const direction =
+ velocity.normalized(); // determine steplength to next volume
+ Point const position = particle.getPosition();
+
+ if (chargeNumber != 0 &&
+ abs(plane.getNormal().dot(velocity.cross(magneticfield))) *
+ 1_s / 1_m / 1_T > 1e-6) {
+
+ auto const* currentLogicalVolumeNode = particle.getNode();
+ auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
+ particle.getPosition());
+ auto k = chargeNumber * constants::c * 1_eV /
+ (particle.getMomentum().getNorm() * 1_V);
+
+ if (direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
+ (plane.getNormal().dot(position - plane.getCenter()) *
+ plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
+ 0) {
+ return Intersections(std::numeric_limits<double>::infinity() * 1_s);
+ }
+
+ LengthType const MaxStepLength1 =
+ (sqrt(direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
+ (plane.getNormal().dot(position - plane.getCenter()) *
+ plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane.getNormal()) / direction.getNorm()) /
+ (plane.getNormal().dot(direction.cross(magneticfield)) * k);
+
+ LengthType const MaxStepLength2 =
+ (-sqrt(
+ direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
+ (plane.getNormal().dot(position - plane.getCenter()) *
+ plane.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane.getNormal()) / direction.getNorm()) /
+ (plane.getNormal().dot(direction.cross(magneticfield)) * k);
+
+ if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
+ return Intersections(std::numeric_limits<double>::infinity() * 1_s);
+ } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
+ CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength2);
+ return Intersections(MaxStepLength2 *
+ (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
+ .getNorm() / absVelocity);
+ } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
+ CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength1);
+ return Intersections(MaxStepLength1 *
+ (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
+ .getNorm() / absVelocity);
+ }
+
+ CORSIKA_LOG_WARN("Particle wasn't tracked with curved trajectory -> straight");
+
+ } // end if curved-tracking
+
+ return tracking_line::Tracking::intersect(particle, plane, medium);
+ }
+
+
template <typename TParticle, typename TBaseNodeType>
inline Intersections Tracking::intersect(const TParticle& particle,
const TBaseNodeType& volumeNode) {
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
index 132a87e5d..7f913bd17 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
@@ -9,8 +9,6 @@
#pragma once
#include <corsika/framework/geometry/Line.hpp>
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/StraightTrajectory.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
index dff4bc18e..6fe3af378 100644
--- a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -61,12 +61,16 @@ namespace corsika {
auto getTrack(TParticle const& particle);
template <typename TParticle, typename TMedium>
- static Intersections intersect(const TParticle& particle, const Sphere& sphere,
- const TMedium& medium);
+ static Intersections intersect(TParticle const& particle, Sphere const& sphere,
+ TMedium const& medium);
template <typename TParticle, typename TBaseNodeType>
- static Intersections intersect(const TParticle& particle,
- const TBaseNodeType& volumeNode);
+ static Intersections intersect(TParticle const& particle,
+ TBaseNodeType const& volumeNode);
+
+ template <typename TParticle, typename TMedium>
+ static Intersections intersect(TParticle const& particle, Plane const& plane,
+ TMedium const& medium);
protected:
/**
diff --git a/corsika/setup/SetupTrajectory.hpp b/corsika/setup/SetupTrajectory.hpp
index f531414f6..95a91115b 100644
--- a/corsika/setup/SetupTrajectory.hpp
+++ b/corsika/setup/SetupTrajectory.hpp
@@ -38,8 +38,8 @@ namespace corsika::setup {
The default tracking algorithm.
*/
- typedef corsika::process::tracking_leapfrog_curved::Tracking Tracking;
- // typedef corsika::process::tracking_leapfrog_straight::Tracking Tracking;
+ typedef corsika::tracking_leapfrog_curved::Tracking Tracking;
+ // typedef corsika::tracking_leapfrog_straight::Tracking Tracking;
// typedef corsika::tracking_line::Tracking Tracking;
/**
@@ -47,6 +47,6 @@ namespace corsika::setup {
*/
/// definition of Trajectory base class, to be used in tracking and cascades
// typedef StraightTrajectory Trajectory;
- typedef corsika::geometry::LeapFrogTrajectory Trajectory;
+ typedef corsika::LeapFrogTrajectory Trajectory;
} // namespace corsika::setup
diff --git a/tests/common/SetupTestTrajectory.hpp b/tests/common/SetupTestTrajectory.hpp
index 0c06d7f39..a33be655e 100644
--- a/tests/common/SetupTestTrajectory.hpp
+++ b/tests/common/SetupTestTrajectory.hpp
@@ -13,33 +13,33 @@
#include <corsika/framework/geometry/Helix.hpp>
#include <corsika/framework/geometry/StraightTrajectory.hpp>
-// #include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
-// #include <corsika/modules/TrackingLine.hpp>
-// #include <corsika/modules/TrackingCurved.hpp> // simple leap-frog implementation
-// #include <corsika/modules/TrackingLeapFrog.hpp> // more complete leap-frog
-// implementation
+#include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
+//#include <corsika/modules/TrackingLine.hpp>
+//#include <corsika/modules/TrackingCurved.hpp> // simple leap-frog implementation
+//#include <corsika/modules/TrackingLeapFrog.hpp> // more complete leap-frog
+ // implementation
namespace corsika::setup::testing {
template <typename TTrack>
- TTrack make_track(Line const& line, TimeType const tEnd);
+ TTrack make_track(Line const line,
+ TimeType const tEnd = std::numeric_limits<TimeType::value_type>::infinity() * 1_s);
template <>
- inline StraightTrajectory make_track<StraightTrajectory>(Line const& line,
+ inline StraightTrajectory make_track<StraightTrajectory>(Line const line,
TimeType const tEnd) {
return StraightTrajectory(line, tEnd);
}
- /*
- template <>
- inline LeapFrogTrajectory make_track<LeapFrogTrajectory>(Line const& line,
- TimeType const tEnd) {
-
- auto const k = square(0_m) / (square(1_s) * 1_V);
- return LeapFrogTrajectory(
- line.getStartPoint(), line.getVelocity(),
- MagneticFieldVector{line.getStartPoint().getCoordinateSystem(), 0_T, 0_T, 0_T},
- k, tEnd);
+ template <>
+ inline LeapFrogTrajectory make_track<LeapFrogTrajectory>(Line const line,
+ TimeType const tEnd) {
+
+ auto const k = square(0_m) / (square(1_s) * 1_V);
+ return LeapFrogTrajectory(
+ line.getStartPoint(), line.getVelocity(),
+ MagneticFieldVector{line.getStartPoint().getCoordinateSystem(), 0_T, 0_T, 0_T},
+ k, tEnd);
}
- */
+
} // namespace corsika::setup::testing
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index 8a7978b84..9070f0d89 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -24,6 +24,8 @@
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
+#include <SetupTestTrajectory.hpp>
+
#include <catch2/catch.hpp>
using namespace corsika;
@@ -68,14 +70,11 @@ public:
auto getTrack(TParticle const& particle) {
VelocityVector const initialVelocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
+ Line const theLine = Line(particle.getPosition(), initialVelocity);
+ TimeType const tEnd = std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
return std::make_tuple(
- StraightTrajectory(
- Line(particle.getPosition(), initialVelocity),
- std::numeric_limits<TimeType::value_type>::infinity() * 1_s), // trajectory,
- // just
- // go
- // ahead
- // forever
+ corsika::setup::testing::make_track<setup::Trajectory>(theLine, tEnd),
+ // trajectory: just go ahead forever
particle.getNode()); // next volume node
}
};
diff --git a/tests/media/testMedium.cpp b/tests/media/testMedium.cpp
index 670c25d92..d303c820f 100644
--- a/tests/media/testMedium.cpp
+++ b/tests/media/testMedium.cpp
@@ -101,7 +101,7 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
auto const tEnd = 1_s;
// and the associated trajectory
- setup::Trajectory const trajectory(line, tEnd);
+ setup::Trajectory const trajectory = corsika::setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
--
GitLab