diff --git a/corsika/detail/modules/tracking/Intersect.inl b/corsika/detail/modules/tracking/Intersect.inl
index 79cfa6984f2bdc33ee5570c24be7e0efd7c8c371..c653cc5aa698107cd25eb3c009e634f5f9bd8f8a 100644
--- a/corsika/detail/modules/tracking/Intersect.inl
+++ b/corsika/detail/modules/tracking/Intersect.inl
@@ -59,9 +59,10 @@ namespace corsika {
for (auto const& node : volumeNode.getChildNodes()) {
Intersections const time_intersections = TDerived::intersect(particle, *node);
+ CORSIKA_LOG_TRACE("intersection times with child volume {}", fmt::ptr(node));
if (!time_intersections.hasIntersections()) { continue; }
- CORSIKA_LOG_DEBUG("intersection times with child volume {} : enter {} s, exit {} s",
- fmt::ptr(node), time_intersections.getEntry() / 1_s,
+ CORSIKA_LOG_TRACE(" : enter {} s, exit {} s",
+ time_intersections.getEntry() / 1_s,
time_intersections.getExit() / 1_s);
auto const t_entry = time_intersections.getEntry();
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 8dc9f1f3a2ce5ae9330135956629c2178868ecd9..415c1ada386c5840dd2c9ded27a28fdf7c3f2236 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -149,21 +149,25 @@ namespace corsika {
}
bool const numericallyInside = sphere.contains(particle.getPosition());
+ CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
auto const absVelocity = velocity.getNorm();
- auto energy = particle.getEnergy();
- auto k = chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V);
+ auto const energy = particle.getEnergy();
+ // this is: k = q/|p|
+ auto const k =
+ chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V);
auto const direction_B_perp = directionBefore.cross(magneticfield);
- auto const denom = direction_B_perp.getSquaredNorm() * k * k;
+ auto const denom = 4. / (direction_B_perp.getSquaredNorm() * k * k);
Vector<length_d> const deltaPos = position - sphere.getCenter();
- double const a = (direction_B_perp.dot(deltaPos) * k + 1) * 4 / (1_m * 1_m * denom);
- double const b = directionBefore.dot(deltaPos) * 8 / (denom * 1_m * 1_m * 1_m);
- double const c =
- (deltaPos.getSquaredNorm() - (sphere.getRadius() * sphere.getRadius())) * 4 /
- (denom * 1_m * 1_m * 1_m * 1_m);
- CORSIKA_LOG_TRACE("denom={}, a={}, b={}, c={}", denom, a, b, c);
- std::complex<double>* solutions = quartic_solver::solve_quartic(0, a, b, c);
+ double const b = (direction_B_perp.dot(deltaPos) * k + 1) * denom / (1_m * 1_m);
+ double const c = directionBefore.dot(deltaPos) * 2 * denom / (1_m * 1_m * 1_m);
+ LengthType const deltaPosLength = deltaPos.getNorm();
+ double const d = (deltaPosLength + sphere.getRadius()) *
+ (deltaPosLength - sphere.getRadius()) * denom /
+ (1_m * 1_m * 1_m * 1_m);
+ CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
+ std::complex<double> const* solutions = quartic_solver::solve_quartic(0, b, c, d);
LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0;
for (int i = 0; i < 4; i++) {
@@ -172,8 +176,8 @@ namespace corsika {
CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
if (numericallyInside) {
// there must be an entry (negative) and exit (positive) solution
- if (dist < -0.0001_m) { // security margin to assure transfer to next
- // logical volume
+ if (dist < 0.0001_m) { // security margin to assure transfer to next
+ // logical volume
if (first_entry == 0) {
d_enter = dist;
} else {
@@ -181,7 +185,7 @@ namespace corsika {
}
first_entry++;
- } else { // thus, dist >= -0.0001_m
+ } else { // thus, dist > -0.0001_m
if (first_exit == 0) {
d_exit = dist;
@@ -196,7 +200,7 @@ namespace corsika {
// both physical solutions (entry, exit) must be positive, and as small as
// possible
- if (dist < -0.0001_m) { // need small numerical margin, to assure transport
+ if (dist < 0.0001_m) { // need small numerical margin, to assure transport
// into next logical volume
continue;
}
@@ -268,6 +272,9 @@ namespace corsika {
LengthType const MaxStepLength1 = (sqrSqrtArg - norm_projected) / denom;
LengthType const MaxStepLength2 = (-sqrSqrtArg - norm_projected) / denom;
+ CORSIKA_LOG_TRACE("MaxStepLength1={}, MaxStepLength2={}", MaxStepLength1,
+ MaxStepLength2);
+
// check: both intersections in past
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
@@ -295,6 +302,9 @@ namespace corsika {
} // end if curved-tracking
+ CORSIKA_LOG_TRACE("straight tracking with chargeNumber={}, B={}", chargeNumber,
+ magneticfield);
+
return tracking_line::Tracking::intersect(particle, plane);
}
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 411ab54ff1b049e0b646599edb0259e4bedcc700..9a4659250286414cc88f837e96c3ac5898952006 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -84,13 +84,14 @@ namespace corsika::tracking_line {
auto const delta = plane.getCenter() - particle.getPosition();
auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
auto const n = plane.getNormal();
- auto const c = n.dot(velocity);
+ auto const n_dot_v = n.dot(velocity);
- CORSIKA_LOG_TRACE("c={}, delta={}, momentum={}", c, delta, particle.getMomentum());
+ CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
+ particle.getMomentum());
- return Intersections(c.magnitude() == 0
+ return Intersections(n_dot_v.magnitude() == 0
? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
- : n.dot(delta) / c);
+ : n.dot(delta) / n_dot_v);
}
} // namespace corsika::tracking_line
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index c04ae2c4c6933d16259851a24afe49dbdc299d51..641fb95f85ce3e61e1779b21037ffc0ae6b3003c 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -6,6 +6,8 @@
* the license.
*/
+#define DEBUG 1
+
/* clang-format off */
// InteractionCounter used boost/histogram, which
// fails if boost/type_traits have been included before. Thus, we have
@@ -129,9 +131,25 @@ int main(int argc, char** argv) {
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
- builder.addLinearLayer(1e9_cm, 112.8_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km+constants::EarthRadius::Mean);
builder.assemble(env);
+ CORSIKA_LOG_DEBUG(
+ "environment setup: universe={}, layer1={}, layer2={}, layer3={}, layer4={}, "
+ "layer5={}",
+ fmt::ptr(env.getUniverse()->getContainingNode(
+ Point(rootCS, {constants::EarthRadius::Mean + 130_km, 0_m, 0_m}))),
+ fmt::ptr(env.getUniverse()->getContainingNode(
+ Point(rootCS, {constants::EarthRadius::Mean + 110_km, 0_m, 0_m}))),
+ fmt::ptr(env.getUniverse()->getContainingNode(
+ Point(rootCS, {constants::EarthRadius::Mean + 50_km, 0_m, 0_m}))),
+ fmt::ptr(env.getUniverse()->getContainingNode(
+ Point(rootCS, {constants::EarthRadius::Mean + 20_km, 0_m, 0_m}))),
+ fmt::ptr(env.getUniverse()->getContainingNode(
+ Point(rootCS, {constants::EarthRadius::Mean + 5_km, 0_m, 0_m}))),
+ fmt::ptr(env.getUniverse()->getContainingNode(
+ Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m}))));
+
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
@@ -159,7 +177,7 @@ int main(int argc, char** argv) {
<< ", norm = " << plab.getNorm() << endl;
auto const observationHeight = 0_km + builder.getEarthRadius();
- auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+ auto const injectionHeight = 111.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));