diff --git a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
index 9ff574321434aa6c7c495ea3fc8db432e4a47879..fd0e29d44532fa19ad5e71569b9bbceee41acfbc 100644
--- a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
+++ b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
@@ -35,7 +35,9 @@ namespace corsika {
}
inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
- return initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+ return (initialDirection_ +
+ initialDirection_.cross(magneticfield_) * timeStep_ * u * k_) *
+ initialVelocity_.getNorm();
}
inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
@@ -44,18 +46,23 @@ namespace corsika {
inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
return u * timeStep_ *
- (double(getVelocity(u).getNorm() / initialVelocity_.getNorm()) + 1.0) / 2;
+ (1. + fabs(0.5 * initialDirection_.cross(magneticfield_).getNorm() * u *
+ timeStep_ * k_));
}
inline LengthType LeapFrogTrajectory::getLength(double const u) const {
- return timeStep_ * initialVelocity_.getNorm() * u;
+ return getDuration(u) * initialVelocity_.getNorm();
}
inline void LeapFrogTrajectory::setLength(LengthType const limit) {
- if (initialVelocity_.getNorm() == 0_m / 1_s) setDuration(0_s);
+ if (initialVelocity_.getNorm() == SpeedType::zero()) setDuration(0_s);
setDuration(limit / initialVelocity_.getNorm());
}
- inline void LeapFrogTrajectory::setDuration(TimeType const limit) { timeStep_ = limit; }
+ inline void LeapFrogTrajectory::setDuration(TimeType const limit) {
+ double const correction =
+ (1. + fabs(0.5 * initialDirection_.cross(magneticfield_).getNorm() * limit * k_));
+ timeStep_ = limit / correction;
+ }
} // namespace corsika
diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
index edfce000082c34a6645c459325f8b45f0a6fcb4e..a48a0da0c64d7c8836ee4761d7788b5410b2b71d 100644
--- a/corsika/detail/framework/utility/CubicSolver.inl
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -288,8 +288,8 @@ namespace corsika {
CORSIKA_LOG_TRACE("niter={}", niter);
if (niter >= maxiter) {
- CORSIKA_LOG_TRACE("failure, no solution");
- // return std::vector<double>{};
+ // CORSIKA_LOG_TRACE("failure, no solution");
+ // return std::vector<double>{};
}
}
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 967b9314198af14e7eb76b9e99d654d60c00c573..2cb39341246d5f6b8bf3fab96515d27fa75516f1 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -63,8 +63,9 @@ namespace corsika {
inline LengthType ObservationPlane<TTracking, TOutput>::getMaxStepLength(
TParticle const& particle, TTrajectory const& trajectory) {
- CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}", particle.asString(),
- particle.getPosition(), particle.getDirection(), plane_.asString());
+ CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, plane={}",
+ particle.asString(), particle.getPosition(),
+ particle.getDirection(), plane_.asString());
auto const intersection = TTracking::intersect(particle, plane_);
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index a2f643cc58d5f72f77e64579157e9a351aa411ed..7772f26609fdd65a3a1216235f4788ba5eda31a8 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -73,7 +73,7 @@ namespace corsika {
// Sternheimer parameterization, density corrections towards high energies
// NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
// MISSING
- CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m{}=",
+ CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m={}",
p.getMomentum().getNorm() / m);
double const x = log10(p.getMomentum().getNorm() / m);
double delta = 0;
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 749693bc7c7a327732ea30d77872006f63eff5da..009be7c63610710d1ab7c1d0113518ff38dc7b68 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -94,6 +94,8 @@ namespace corsika {
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
+ CORSIKA_LOG_TRACE("p_perp={} eV", p_perp / 1_eV);
+
if (p_perp < 1_eV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
@@ -121,12 +123,17 @@ namespace corsika {
decltype(1 / (tesla * second)) const k =
charge / p_norm *
initialVelocity.getNorm(); // since we use steps in time and not length
- // units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1/(T s)
+ // units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1 / (T*s)
+
+ // correction factor to move from leap-frog length L to half-step deltaL/2
+ double const correction =
+ 1. + fabs(0.5 * particle.getDirection().cross(magneticfield).getNorm() * k *
+ minTime);
return std::make_tuple(
LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
- minTime), // trajectory
- minNode); // next volume node
+ minTime / correction), // --> trajectory
+ minNode); // --> next volume node
}
template <typename TParticle>
@@ -164,7 +171,7 @@ namespace corsika {
constants::c * convert_HEP_to_SI<MassType::dimension_type>(
particle.getMomentum().getNorm()); // km * m /s
// this is: k = q/|p|
- auto const k = charge / p_norm;
+ decltype(1 / (tesla * meter)) const k = charge / p_norm;
MagneticFieldVector const direction_x_B = directionBefore.cross(magneticfield);
auto const denom = 4. / (direction_x_B.getSquaredNorm() * k * k);
@@ -176,6 +183,10 @@ namespace corsika {
(deltaPosLength - sphere.getRadius()) * denom /
(1_m * 1_m * 1_m * 1_m);
CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
+ // solutions of deltaL are obtained from quartic equation. Note, deltaL/2 is the
+ // length of each half step, however, the second half step is slightly longer
+ // because of the non-conservation of norm/velocity.
+ // The leap-frog length L is deltaL/2 * (1+|u_{n+1}|)
std::vector<double> solutions = solve_quartic_real(1, 0, b, c, d);
if (!solutions.size()) {
return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
@@ -184,7 +195,11 @@ namespace corsika {
int first = 0, first_entry = 0, first_exit = 0;
for (auto solution : solutions) {
LengthType const dist = solution * 1_m;
- CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
+ CORSIKA_LOG_TRACE(
+ "Solution (real) for current Volume: deltaL/2*2={} (deltaL/2*2/v={}, "
+ "corrected={}) ",
+ dist, dist / absVelocity,
+ dist / absVelocity * (1. + fabs(0.5 * direction_x_B.getNorm() * k * dist)));
if (numericallyInside) {
// there must be an entry (negative) and exit (positive) solution
if (dist < 0.0001_m) { // security margin to assure
@@ -242,7 +257,12 @@ namespace corsika {
first_entry, first_exit);
return Intersections();
}
- return Intersections(d_enter / absVelocity, d_exit / absVelocity);
+ // return in units of time, correct to leap-frog length L
+ double const corr_enter = 1. + fabs(0.5 * direction_x_B.getNorm() * k * d_enter);
+ double const corr_exit = 1. + fabs(0.5 * direction_x_B.getNorm() * k * d_exit);
+
+ return Intersections(d_enter * corr_enter / absVelocity,
+ d_exit * corr_exit / absVelocity);
}
template <typename TParticle>
@@ -253,7 +273,7 @@ namespace corsika {
ElectricChargeType const charge = particle.getCharge();
- if (charge != 0 * constants::e) {
+ if (charge != ElectricChargeType::zero()) {
auto const* currentLogicalVolumeNode = particle.getNode();
VelocityVector const velocity = particle.getVelocity();
@@ -284,6 +304,10 @@ namespace corsika {
plane.getNormal().dot(position - plane.getCenter())) // unit: kg*m/s *m
/ (1_m * 1_m * 1_kg) * 1_s;
+ // deltaL from quadratic solution return half-step length deltaL/2 for leap-frog
+ // algorithmus. Note, the leap-frog length L is longer by (1+|u_{n_1}|)/2 because
+ // the direction norm of the second half step is >1.
+
std::vector<double> const deltaLs = solve_quadratic_real(denom, p, q);
CORSIKA_LOG_TRACE("deltaLs=[{}]", fmt::join(deltaLs, ", "));
@@ -310,17 +334,14 @@ namespace corsika {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
}
- CORSIKA_LOG_TRACE("maxStepLength={}", maxStepLength);
-
- // with final length correction
- auto const corr =
- (direction + direction_x_B * maxStepLength * charge / (p_norm * 2))
- .getNorm() /
- absVelocity; // unit: s/m
+ CORSIKA_LOG_TRACE("maxStepLength={} s", maxStepLength / 1_s);
- return Intersections(maxStepLength * corr); // unit: s
+ // with final length correction, |direction| becomes >1 during step
+ double const corr =
+ 1. + fabs(0.5 * direction_x_B.getNorm() * maxStepLength * charge / p_norm);
- } // no charge
+ return Intersections(maxStepLength * corr / absVelocity); // unit: s
+ } // no charge
CORSIKA_LOG_TRACE("(plane) straight tracking with charge={}, B={}", charge,
particle.getNode()->getModelProperties().getMagneticField(
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
index 513ffbaf509d2c05a7020f038d0ab320592c4712..c4ed5aafd3883b1ff4f5ba859c8a418018fc50b3 100644
--- a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -33,8 +33,8 @@ namespace corsika {
*
* Performs one leap-frog step consistent of two halve-steps with steplength/2
* The step is caluculated analytically precisely to reach to the next volume
- *boundary.
- **/
+ * boundary.
+ */
template <typename TParticle>
auto make_LeapFrogStep(TParticle const& particle, LengthType steplength);
@@ -59,16 +59,42 @@ namespace corsika {
template <typename TParticle>
auto getTrack(TParticle const& particle);
- //! find intersection of Sphere with Track
+ /**
+ * find intersection of Sphere with Track
+ *
+ * Returns intersection of particle assuming a curved leap-frog step, with a sphere.
+ * Entry and exit times are calculated, where the velocity is constant and the
+ * steplength is the geometric steplength of the leap-frog.
+ *
+ * @param particle Current particle state
+ * @param sphere Sphere object
+ */
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Sphere const& sphere);
- //! find intersection of Volume node with Track of particle
+ /**
+ * find intersection of any Volume node with particle
+ *
+ * The intersection time(s) of a particle, assuming a curved leap-frog
+ * step, are calculated for any volume type.
+ *
+ */
template <typename TParticle, typename TBaseNodeType>
static Intersections intersect(TParticle const& particle,
TBaseNodeType const& node);
- //! find intersection of Plane with Track
+ /**
+ * find intersection of Plane with Track
+ *
+ * Intersection times of particle are caculated with a plane, assuming a curved leap
+ * frog trajectory. The intersection time is assuming constant velocity (no change)
+ * along the geometric leap-frog step.
+ *
+ * @tparam TParticle Type of particle object on stack.
+ * @param particle Particle initial state.
+ * @param plane Plane.
+ * @return Intersections in time units.
+ */
template <typename TParticle>
static Intersections intersect(TParticle const& particle, Plane const& plane);
@@ -80,7 +106,6 @@ namespace corsika {
* Use internally stored class tracking_line::Tracking to
* perform a straight line tracking, if no magnetic bendig was
* detected.
- *
*/
template <typename TParticle>
auto getLinearTrajectory(TParticle& particle);
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index abafdf0c3ab305d008c0ef6d13cb3885e6d08810..f283679183598d43151fed5841ec2cc69581f516 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -202,8 +202,8 @@ int main(int argc, char** argv) {
<< std::endl;
OutputManager output("hybrid_MC_outputs");
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env, true,
- 1000};
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
+ false, 1000};
// setup processes, decays and interactions