diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 021b8c50aa1ae690df165b92ce68a0021df59a96..32127908c714d79dfdd7e22e74615aed449b1b7c 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -269,8 +269,8 @@ namespace corsika {
if (actual_decay_time * 0.99 > initial_inv_decay_time) {
CORSIKA_LOG_WARN(
"Decay time decreased during step! This leads to un-physical step length. "
- "initial_decay_time={}, actual_decay_time={}",
- 1 / initial_inv_decay_time, 1 / actual_decay_time);
+ "delta_inverse_decay_time={}",
+ 1 / initial_inv_decay_time - 1 / actual_decay_time);
}
#endif
@@ -303,9 +303,8 @@ namespace corsika {
if (actual_inv_length * 0.99 > initial_inv_int_length) {
CORSIKA_LOG_WARN(
"Interaction length decreased during step! This leads to un-physical step "
- "length. "
- "initial_inv_int_length={}, actual_inv_length={}",
- 1 / initial_inv_int_length, 1 / actual_inv_length);
+ "length. delta_innverse_interaction_length={}",
+ 1 / initial_inv_int_length - 1 / actual_inv_length);
}
#endif
diff --git a/corsika/detail/framework/utility/LinearSolver.inl b/corsika/detail/framework/utility/LinearSolver.inl
index e9583acf2f08514e881ddc4a30b2db0b8852fb25..b2fc96ee7ae9a2a6947f1858b7c097e001b920bc 100644
--- a/corsika/detail/framework/utility/LinearSolver.inl
+++ b/corsika/detail/framework/utility/LinearSolver.inl
@@ -14,7 +14,7 @@
namespace corsika {
- inline std::vector<double> solve_linear_real(double a, double b, double const epsilon) {
+ inline std::vector<double> solve_linear_real(double a, double b) {
if (a == 0) {
return {}; // no (b!=0), or infinite number (b==0) of solutions....
@@ -23,8 +23,7 @@ namespace corsika {
return {-b / a};
}
- inline std::vector<std::complex<double>> solve_linear(double a, double b,
- double const epsilon) {
+ inline std::vector<std::complex<double>> solve_linear(double a, double b) {
if (std::abs(a) == 0) {
return {}; // no (b!=0), or infinite number (b==0) of solutions....
diff --git a/corsika/detail/framework/utility/QuadraticSolver.inl b/corsika/detail/framework/utility/QuadraticSolver.inl
index d3b3f86714af097745e9879c3543244ddf4fbcb4..9e67ea031c11d6160602a520746c4e3aead33a48 100644
--- a/corsika/detail/framework/utility/QuadraticSolver.inl
+++ b/corsika/detail/framework/utility/QuadraticSolver.inl
@@ -11,10 +11,10 @@ namespace corsika {
inline std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
long double c,
double const epsilon) {
- if (std::abs(a) < epsilon) { return solve_linear(b, c, epsilon); }
+ if (std::abs(a) < epsilon) { return solve_linear(b, c); }
if (std::abs(c) < epsilon) {
- std::vector<std::complex<double>> lin_result = solve_linear(a, b, epsilon);
+ std::vector<std::complex<double>> lin_result = solve_linear(a, b);
lin_result.push_back({0.});
return lin_result;
}
@@ -45,9 +45,9 @@ namespace corsika {
CORSIKA_LOG_TRACE("quadratic: a={} b={} c={}", a, b, c);
- if (std::abs(a) < epsilon) { return solve_linear_real(b, c, epsilon); }
+ if (std::abs(a) < epsilon) { return solve_linear_real(b, c); }
if (std::abs(c) < epsilon) {
- std::vector<double> lin_result = solve_linear_real(a, b, epsilon);
+ std::vector<double> lin_result = solve_linear_real(a, b);
lin_result.push_back(0.);
return lin_result;
}
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index cad673519f8de048c7df4bdaeb3a62e55be126b1..35ebad050be46604a48a364ecefe8b9f7f35e449 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -61,14 +61,14 @@ namespace corsika {
auto const position = particle.getPosition();
CORSIKA_LOG_DEBUG(
- "Tracking pid: {}"
- " , E = {} GeV",
- particle.getPID(), particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking pos: {}", position.getCoordinates());
- CORSIKA_LOG_DEBUG("Tracking E: {} GeV", particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking p: {} GeV",
- particle.getMomentum().getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking v: {} ", initialVelocity.getComponents());
+ "TrackingLeapfrog_Curved pid: {}"
+ " , E = {} GeV \n"
+ "\tTracking pos: {} \n"
+ "\tTracking p: {} GeV \n"
+ "\tTracking v: {}",
+ particle.getPID(), particle.getEnergy() / 1_GeV, position.getCoordinates(),
+ particle.getMomentum().getComponents() / 1_GeV,
+ initialVelocity.getComponents());
typedef
typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
@@ -84,7 +84,10 @@ namespace corsika {
ElectricChargeType const charge = particle.getCharge();
bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T;
- if (no_deflection) { return getLinearTrajectory(particle); }
+ if (no_deflection) {
+ CORSIKA_LOG_TRACE("no_deflection");
+ return getLinearTrajectory(particle);
+ }
HEPMomentumType const p_perp =
(particle.getMomentum() -
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
index 77cef175de6d136b8681456752fdbaec725c03d4..4e7efed775468a6525a236fe7b5a72b45896e518 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
@@ -29,15 +29,17 @@ namespace corsika {
particle.getMomentum() / particle.getEnergy() * constants::c;
Point const initialPosition = particle.getPosition();
+
CORSIKA_LOG_DEBUG(
- "TrackingB pid: {}"
- " , E = {} GeV",
- particle.getPID(), particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB pos: {}", initialPosition.getCoordinates());
- CORSIKA_LOG_DEBUG("TrackingB E: {} GeV", particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB p: {} GeV",
- particle.getMomentum().getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB v: {} ", initialVelocity.getComponents());
+ "TrackingLeapFrogStraight pid: {}"
+ " , E = {} GeV \n"
+ "\tTracking pos: {} \n"
+ "\tTracking p: {} GeV \n"
+ "\tTracking v: {} ",
+ particle.getPID(), particle.getEnergy() / 1_GeV,
+ initialPosition.getCoordinates(),
+ particle.getMomentum().getComponents() / 1_GeV,
+ initialVelocity.getComponents());
typedef decltype(particle.getNode()) node_type;
node_type const volumeNode = particle.getNode();
@@ -93,7 +95,7 @@ namespace corsika {
// need to follow strongly curved trajectories segment-wise,
// at least if we don't employ concepts as "Helix
// Trajectories" or similar
- double const maxRadians = 0.01;
+ double const maxRadians = 0.001;
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 36e0975605bbdc48b5335a36516de059a5a41d34..7f41291885e628bf62360c5e870f95ddd78f2873 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -31,14 +31,13 @@ namespace corsika::tracking_line {
auto const initialPosition = particle.getPosition();
CORSIKA_LOG_DEBUG(
- "Tracking pid: {}"
- " , E = {} GeV",
- particle.getPID(), particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking pos: {}", initialPosition.getCoordinates());
- CORSIKA_LOG_DEBUG("Tracking E: {} GeV", particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking p: {} GeV",
- particle.getMomentum().getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking v: {} ", initialVelocity.getComponents());
+ "TrackingStraight pid: {}"
+ " , E = {} GeV \n"
+ "\tTracking pos: {} \n"
+ "\tTracking p: {} GeV \n"
+ "\tTracking v: {}",
+ particle.getPID(), particle.getEnergy() / 1_GeV, initialPosition.getCoordinates(),
+ particle.getMomentum().getComponents() / 1_GeV, initialVelocity.getComponents());
// traverse the environment volume tree and find next
// intersection
@@ -52,7 +51,8 @@ namespace corsika::tracking_line {
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Sphere const& sphere) {
- auto const delta = particle.getPosition() - sphere.getCenter();
+ auto const position = particle.getPosition();
+ auto const delta = position - sphere.getCenter();
auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
auto const vSqNorm = velocity.getSquaredNorm();
auto const R = sphere.getRadius();
@@ -64,6 +64,9 @@ namespace corsika::tracking_line {
if (discriminant.magnitude() > 0) {
auto const sqDisc = sqrt(discriminant);
auto const invDenom = 1 / vSqNorm;
+
+ bool const numericallyInside = sphere.contains(position);
+ CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
return Intersections((-vDotDelta - sqDisc) * invDenom,
(-vDotDelta + sqDisc) * invDenom);
}
diff --git a/corsika/framework/core/Logging.hpp b/corsika/framework/core/Logging.hpp
index 67f059c2fe6d1dc5d4b228e551f4ebbbcd1acf9e..d127040b036dab150ce06a40644050f20a3fab2a 100644
--- a/corsika/framework/core/Logging.hpp
+++ b/corsika/framework/core/Logging.hpp
@@ -49,7 +49,8 @@ namespace corsika {
/*
* The default pattern for CORSIKA8 loggers.
*/
- const std::string default_pattern{"[%n:%^%-8l%$] %v"};
+ const std::string minimal_pattern{"[%n:%^%-8l%$] %v"};
+ const std::string default_pattern{"[%n:%^%-8l%$(%s:%#)] %v"};
const std::string source_pattern{"[%n:%^%-8l%$(%s:%!:%#)] %v"};
/**
@@ -91,7 +92,6 @@ namespace corsika {
*/
static inline std::shared_ptr<spdlog::logger> corsika_logger =
get_logger("corsika", true);
- // corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// many of these free functions are special to the logging
// infrastructure so we hide them in the corsika::logging namespace.
diff --git a/corsika/framework/utility/LinearSolver.hpp b/corsika/framework/utility/LinearSolver.hpp
index 95322958d09c744ccc4971930fec2616a7d1dde2..754039571287a1a816c77ec741b7d850b87c46da 100644
--- a/corsika/framework/utility/LinearSolver.hpp
+++ b/corsika/framework/utility/LinearSolver.hpp
@@ -14,10 +14,9 @@
namespace corsika {
- std::vector<double> solve_linear_real(double a, double b, double const epsilon = 1e-12);
+ std::vector<double> solve_linear_real(double a, double b);
- std::vector<std::complex<double>> solve_linear(double a, double b,
- double const epsilon = 1e-12);
+ std::vector<std::complex<double>> solve_linear(double a, double b);
} // namespace corsika
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index 3f1d6ecaa9978bdf0a9f927cb6e1c760e45b7e08..f09297e7699483f6fb25b9a9d6514eefcb76d6d9 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -20,15 +20,22 @@
#include <catch2/catch.hpp>
+#include <boost/type_index.hpp>
+
using namespace corsika;
+struct NonExistingDummyObject : public IVolume {
+ NonExistingDummyObject const& getVolume() const { return *this; }
+ bool contains(Point const&) const { return false; }
+};
+
template <typename T>
int sgn(T val) {
return (T(0) < val) - (val < T(0));
}
/**
- \file testTracking.cpp
+ @file testTracking.cpp
This is the unified and commond unit test for all Tracking algorithms:
@@ -36,14 +43,29 @@ int sgn(T val) {
- tracking_leapfrog_straight::Tracking
- tracking_line::Tracking
+
+ The main part of tests are to inject particles at 10GeV momentum at
+ (-Rg,0,0) in +x direction into a sphere of radius Rg, where Rg is
+ the gyroradius (or 10m for neutral particles). Then it is checked
+ where the particles leaves the sphere for different charges
+ (-1,0,+1) and field strength (-50uT, 0T, +50uT).
+
+ Each test is perfromed once, with the particles starting logically
+ outside of the Rg sphere (thus it first has to enter insides) and a
+ second time with the particle already logically inside the sphere.
+
+ There is a second smaller, internal sphere at +z displacement. Only
+ neutral particles are allowed and expected to hit this.
+
+ All those tests are parameterized, thus, they can be easily extended
+ or applied to new algorithms.
+
*/
-TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
- tracking_leapfrog_curved::Tracking,
+TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
logging::set_level(logging::level::info);
- logging::set_level(logging::level::trace);
const HEPEnergyType P0 = 10_GeV;
@@ -67,9 +89,10 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
Bfield / 1_uT, outer)) {
CORSIKA_LOG_DEBUG(
- "********************\n TEST section PID={}, Bfield={} "
- "uT, outer (?)={}",
- PID, Bfield / 1_uT, outer);
+ "********************\n TEST algo={} section PID={}, "
+ "Bfield={} "
+ "uT, start_outside={}",
+ boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT, outer);
const int chargeNumber = get_charge_number(PID);
LengthType radius = 10_m;
@@ -90,6 +113,12 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
TestType tracking;
Point const center(cs, {0_m, 0_m, 0_m});
auto target = setup::Environment::createNode<Sphere>(center, radius);
+ // every particle should hit target_2
+ auto target_2 = setup::Environment::createNode<Sphere>(
+ Point(cs, {-radius * 3 / 4, 0_m, 0_m}), radius * 0.2);
+ // only neutral particles hit_target_neutral
+ auto target_neutral = setup::Environment::createNode<Sphere>(
+ Point(cs, {radius / 2, 0_m, 0_m}), radius * 0.1);
using MyHomogeneousModel = MediumPropertyModel<
UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
@@ -98,8 +127,18 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
target->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ target_neutral->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ target_2->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
auto* targetPtr = target.get();
+ auto* targetPtr_2 = target_2.get();
+ auto* targetPtr_neutral = target_neutral.get();
worldPtr->addChild(std::move(target));
+ targetPtr->addChild(std::move(target_2));
+ targetPtr->addChild(std::move(target_neutral));
auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs);
{ [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
@@ -126,7 +165,11 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
CHECK(nextVol == targetPtr);
}
// move forward, until we leave target volume
- while (nextVol == targetPtr) {
+ bool hit_neutral = false;
+ bool hit_2nd = false;
+ while (nextVol != worldPtr) {
+ if (nextVol == targetPtr_neutral) hit_neutral = true;
+ if (nextVol == targetPtr_2) hit_2nd = true;
const auto [traj2, nextVol2] = tracking.getTrack(particle);
nextVol = nextVol2;
particle.setNode(nextVol);
@@ -139,6 +182,8 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
traj2.getLength(1) / particle.getVelocity().getNorm());
}
CHECK(nextVol == worldPtr);
+ CHECK(hit_2nd == true);
+ CHECK(hit_neutral == (deflect == 0 ? true : false));
Point pointCheck(cs, (deflect == 0 ? radius : 0_m), (deflect * radius), 0_m);
@@ -151,3 +196,165 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
Approx(0).margin(1e-3));
}
}
+
+/** specifc test for curved leap-frog algorithm **/
+
+TEST_CASE("TrackingLeapFrogCurved") {
+
+ logging::set_level(logging::level::info);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ corsika::Code PID = Code::MuPlus;
+
+ using MyHomogeneousModel = MediumPropertyModel<
+ UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+ SECTION("infinite sphere / universe") {
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, 100_uT);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ tracking_leapfrog_curved::Tracking tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, worldPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in X-direction
+
+ particle.setPosition(Point(cs, 0_m, 0_m, 0_m));
+
+ Sphere sphere(center, 1_km * std::numeric_limits<double>::infinity());
+
+ auto intersections = tracking.intersect(particle, sphere);
+ // this must be a "linear trajectory" with no curvature
+
+ CHECK(intersections.hasIntersections() == false);
+ }
+
+ SECTION("momentum along field") {
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, 100_uT);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ tracking_leapfrog_curved::Tracking tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+ auto target = setup::Environment::createNode<Sphere>(center, 10_km);
+
+ MagneticFieldVector magneticfield(cs, 100_T, 0_T, 0_uT);
+ target->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ auto* targetPtr = target.get();
+ worldPtr->addChild(std::move(target));
+
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in X-direction
+
+ particle.setNode(targetPtr); // set particle outside "target" volume
+ particle.setPosition(Point(cs, 0_m, 0_m, 0_m));
+
+ auto [traj, nextVol] = tracking.getTrack(particle);
+ // this must be a "linear trajectory" with no curvature
+
+ CHECK(traj.getDirection(0).getComponents() == traj.getDirection(1).getComponents());
+ }
+}
+
+TEMPLATE_TEST_CASE("TrackingFail", "doesntwork", tracking_leapfrog_curved::Tracking,
+ tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
+
+ logging::set_level(logging::level::info);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, 1_mT);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ TestType tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+
+ auto [stack, viewPtr] =
+ setup::testing::setup_stack(Code::Proton, 0, 0, P0, worldPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ NonExistingDummyObject const dummy;
+ CHECK_THROWS(tracking.intersect(particle, dummy));
+}
+
+TEMPLATE_TEST_CASE("TrackingPlane", "plane", tracking_leapfrog_curved::Tracking,
+ tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
+
+ logging::set_level(logging::level::trace);
+
+ const HEPEnergyType P0 = 10_GeV;
+
+ auto PID = GENERATE(as<Code>{}, Code::MuPlus, Code::MuPlus, Code::Gamma);
+ // for algorithms that know magnetic deflections choose: +-50uT, 0uT
+ // otherwise just 0uT
+ auto Bfield = GENERATE(filter(
+ []([[maybe_unused]] MagneticFluxType v) {
+ if constexpr (std::is_same_v<TestType, tracking_line::Tracking>)
+ return v == 0_uT;
+ else
+ return true;
+ },
+ values<MagneticFluxType>({50_uT, 0_uT, -50_uT})));
+
+ SECTION(fmt::format("Tracking PID={}, Bfield={} uT", PID, Bfield / 1_uT)) {
+
+ CORSIKA_LOG_DEBUG(
+ "********************\n TEST algo={} section PID={}, "
+ "Bfield={}uT ",
+ boost::typeindex::type_id<TestType>().pretty_name(), PID, Bfield / 1_uT);
+
+ const int chargeNumber = get_charge_number(PID);
+ LengthType radius = 10_m;
+ int deflect = 0;
+ if (chargeNumber != 0 and Bfield != 0_T) {
+ deflect = 1;
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(P0) *
+ constants::c / (abs(get_charge(PID)) * abs(Bfield)));
+ CORSIKA_LOG_DEBUG("Rg={} deflect={}", gyroradius, deflect);
+ radius = gyroradius;
+ }
+
+ auto [env, csPtr, worldPtr] =
+ corsika::setup::testing::setup_environment(Code::Oxygen, Bfield);
+ { [[maybe_unused]] const auto& env_dummy = env; }
+ auto const& cs = *csPtr;
+
+ TestType tracking;
+ Point const center(cs, {0_m, 0_m, 0_m});
+
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, worldPtr, cs);
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ auto particle = stack->first();
+ // Note: momentum in X-direction
+ // magnetic field in Z-direction
+
+ particle.setPosition(Point(cs, -radius, 0_m, 0_m));
+
+ // put plane where we expect deflection by radius/2
+ Plane const plane(Point(cs, radius * (1 - sqrt(3. / 4.)), 0_m, 0_m),
+ DirectionVector(cs, {-1., 0., 0.}));
+ Intersections const hit = tracking.intersect(particle, plane);
+
+ CORSIKA_LOG_DEBUG("entry={} exit={}", hit.getEntry(), hit.getExit());
+
+ CHECK(hit.hasIntersections());
+ CHECK(hit.getEntry() / 1_s == Approx(0.00275 * deflect).margin(0.0003));
+ CHECK(hit.getExit() == 1_s * std::numeric_limits<double>::infinity());
+ }
+}