diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 21cdbe3533df655fd047cbc74acc1de49114bcb5..c40823c142d1640484e916467ea76d9d81c567cf 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -8,7 +8,6 @@
 
 #include <corsika/cascade/Cascade.h>
 #include <corsika/process/ProcessSequence.h>
-#include <corsika/process/tracking_line/Tracking.h>
 
 #include <corsika/setup/SetupEnvironment.h>
 #include <corsika/setup/SetupStack.h>
@@ -123,7 +122,7 @@ int main() {
   universe.AddChild(std::move(world));
 
   // setup processes, decays and interactions
-  tracking_line::Tracking tracking;
+  setup::Tracking tracking;
 
   random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
   process::sibyll::Interaction sibyll;
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 2707e3dde5a55180016ba5fe7edc402582970449..755c7c7c14d0c18d920753a1ed1c024ef49cb0b3 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -10,7 +10,6 @@
 #include <corsika/process/ProcessSequence.h>
 #include <corsika/process/energy_loss/EnergyLoss.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/Tracking.h>
 
 #include <corsika/setup/SetupEnvironment.h>
 #include <corsika/setup/SetupStack.h>
@@ -135,7 +134,7 @@ int main() {
   }
 
   // setup processes, decays and interactions
-  tracking_line::Tracking tracking;
+  setup::Tracking tracking;
   stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
 
   random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index 11b3dac77a54b21990abe3eb0cdc196b34d59d83..909cad004543a93c090481289a9ff134b5751af0 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -10,7 +10,6 @@
 #include <corsika/process/ProcessSequence.h>
 #include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/Tracking.h>
 
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
@@ -121,7 +120,7 @@ int main() {
   }
 
   // setup processes, decays and interactions
-  tracking_line::Tracking tracking;
+  setup::Tracking tracking;
   stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
 
   random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index d834afe77891afba77c490f0b21fab6652db5c57..31298aa288e442e49eb87f0f007a2fc88acfa087 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -21,7 +21,6 @@
 #include <corsika/process/proposal/ContinuousProcess.h>
 #include <corsika/process/proposal/Interaction.h>
 #include <corsika/process/track_writer/TrackWriter.h>
-#include <corsika/process/tracking_line/Tracking.h>
 #include <corsika/random/RNGManager.h>
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
@@ -157,7 +156,7 @@ int main(int argc, char** argv) {
   auto sequence = process::sequence(proposalCounted, em_continuous, longprof, cut,
                                     observationLevel, trackWriter);
   // define air shower object, run simulation
-  tracking_line::Tracking tracking;
+  setup::Tracking tracking;
   cascade::Cascade EAS(env, tracking, sequence, stack);
 
   // to fix the point of first interaction, uncomment the following two lines:
diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc
index a9f7432ae1db26c1e5f4c3e98717142de7e61589..75ed47a1914b7d2b7323f294ae2b071f3ba10ec7 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/Documentation/Examples/hybrid_MC.cc
@@ -34,7 +34,6 @@
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
 #include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/process/tracking_line/Tracking.h>
 #include <corsika/process/urqmd/UrQMD.h>
 #include <corsika/random/RNGManager.h>
 #include <corsika/setup/SetupStack.h>
@@ -242,7 +241,7 @@ int main(int argc, char** argv) {
                                     eLoss, cut, conex, longprof, observationLevel);
 
   // define air shower object, run simulation
-  tracking_line::Tracking tracking;
+  setup::Tracking tracking;
   cascade::Cascade EAS(env, tracking, sequence, stack);
 
   // to fix the point of first interaction, uncomment the following two lines:
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 235c5b3f5156224acef15a301b3cd84a1a306dca..834989466164aa4a4880585c27ed7d58e3059676 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -38,9 +38,9 @@
 #include <corsika/process/proposal/Interaction.h>
 #include <corsika/process/pythia/Decay.h>
 #include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/process/sibyll/Interaction.h>
 #include <corsika/process/sibyll/NuclearInteraction.h>
-#include <corsika/process/tracking_line/Tracking.h>
 #include <corsika/process/urqmd/UrQMD.h>
 #include <corsika/random/RNGManager.h>
 #include <corsika/setup/SetupStack.h>
@@ -85,7 +85,7 @@ using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagnetic
 
 int main(int argc, char** argv) {
 
-  logging::SetLevel(logging::level::info);
+  logging::SetLevel(logging::level::trace);
 
   C8LOG_INFO("vertical_EAS");
 
@@ -110,7 +110,7 @@ int main(int argc, char** argv) {
       setup::EnvironmentInterface,
       MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
                           environment::Medium::AirDry1Atm,
-                          geometry::Vector{rootCS, 0_T, 50_uT, 0_T});
+                          geometry::Vector{rootCS, 0_T, 5000_mT, 0_T});
   builder.setNuclearComposition(
       {{particles::Code::Nitrogen, particles::Code::Oxygen},
        {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
@@ -255,12 +255,13 @@ int main(int argc, char** argv) {
       process::select(urqmdCounted, process::sequence(sibyllNucCounted, sibyllCounted),
                       EnergySwitch(55_GeV));
   auto decaySequence = process::sequence(decayPythia, decaySibyll);
+  stack_inspector::StackInspector<setup::Stack> stackInspect(1000, false, E0);
   auto sequence =
-      process::sequence(hadronSequence, reset_particle_mass, decaySequence,
-                        proposalCounted, em_continuous, cut, observationLevel, longprof);
+    process::sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
+                        proposalCounted, em_continuous, cut, trackWriter, observationLevel, longprof);
 
   // define air shower object, run simulation
-  tracking_line::Tracking tracking;
+  setup::Tracking tracking;
   cascade::Cascade EAS(env, tracking, sequence, stack);
 
   // to fix the point of first interaction, uncomment the following two lines:
diff --git a/Environment/BaseExponential.h b/Environment/BaseExponential.h
index 94b454bf686f4b0657b2e5692df8fb0e3b287816..bf8c16d961c4cfdeb28827851dd8e20e8040c17e 100644
--- a/Environment/BaseExponential.h
+++ b/Environment/BaseExponential.h
@@ -10,9 +10,9 @@
 
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupTrajectory.h>
 
 #include <limits>
 
@@ -47,7 +47,7 @@ namespace corsika::environment {
      */
     // clang-format on
     units::si::GrammageType IntegratedGrammage(
-        geometry::LineTrajectory const& vLine, units::si::LengthType vL,
+        setup::Trajectory const& vLine, units::si::LengthType vL,
         geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
       if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); }
 
@@ -80,7 +80,7 @@ namespace corsika::environment {
      */
     // clang-format on
     units::si::LengthType ArclengthFromGrammage(
-        geometry::LineTrajectory const& vLine, units::si::GrammageType vGrammage,
+        setup::Trajectory const& vLine, units::si::GrammageType vGrammage,
         geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
       auto const uDotA = vLine.GetDirection(0).dot(vAxis).magnitude();
       auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetLine().GetR0());
diff --git a/Environment/DensityFunction.h b/Environment/DensityFunction.h
index 5f61e26d7daf2502e47d8f6bc596d9cd484e5295..3ac33c3f21ebc0a8fbe792d9641b3a37c4dec7c4 100644
--- a/Environment/DensityFunction.h
+++ b/Environment/DensityFunction.h
@@ -11,7 +11,6 @@
 #include <corsika/environment/LinearApproximationIntegrator.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
 
 namespace corsika::environment {
 
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
index 31f3309943f3cc2c927b80dc4f29a32d8c61c82d..8aeafbce4a776c1ae1dd4b335c3578962a21e215 100644
--- a/Environment/FlatExponential.h
+++ b/Environment/FlatExponential.h
@@ -12,10 +12,11 @@
 #include <corsika/environment/NuclearComposition.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/setup/SetupTrajectory.h>
+
 namespace corsika::environment {
 
   //clang-format off
@@ -51,13 +52,13 @@ namespace corsika::environment {
 
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
 
-    units::si::GrammageType IntegratedGrammage(geometry::LineTrajectory const& vLine,
+    units::si::GrammageType IntegratedGrammage(setup::Trajectory const& vLine,
                                                units::si::LengthType vTo) const override {
       return Base::IntegratedGrammage(vLine, vTo, fAxis);
     }
 
     units::si::LengthType ArclengthFromGrammage(
-        geometry::LineTrajectory const& vLine,
+        setup::Trajectory const& vLine,
         units::si::GrammageType vGrammage) const override {
       return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
     }
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index d630d8a6b403959289d74184f48839595d5cbfce..f39cc70610c7e9fdca38172e51f1d9c8179f494b 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -16,6 +16,8 @@
 #include <corsika/random/RNGManager.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/setup/SetupTrajectory.h>
+
 #include <cassert>
 
 /**
@@ -42,14 +44,14 @@ namespace corsika::environment {
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
 
     corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::LineTrajectory const&,
+        corsika::setup::Trajectory const&,
         corsika::units::si::LengthType pTo) const override {
       using namespace corsika::units::si;
       return pTo * fDensity;
     }
 
     corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::LineTrajectory const&,
+        corsika::setup::Trajectory const&,
         corsika::units::si::GrammageType pGrammage) const override {
       return pGrammage / fDensity;
     }
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index 388f79c4366ef59f923b1ad89d3d111fc2fcd62b..7a3b951711b8fe19f928eed3c5cf3d21a17e73d6 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -11,9 +11,10 @@
 #include <corsika/environment/NuclearComposition.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/setup/SetupTrajectory.h>
+
 namespace corsika::environment {
 
   class IMediumModel {
@@ -26,11 +27,11 @@ namespace corsika::environment {
     // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
     // approach; for now, only lines are supported
     virtual corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::LineTrajectory const&,
+        corsika::setup::Trajectory const&,
         corsika::units::si::LengthType) const = 0;
 
     virtual corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::LineTrajectory const&,
+        corsika::setup::Trajectory const&,
         corsika::units::si::GrammageType) const = 0;
 
     virtual NuclearComposition const& GetNuclearComposition() const = 0;
diff --git a/Environment/InhomogeneousMedium.h b/Environment/InhomogeneousMedium.h
index 87002d6b109ef83cbd427e2680af8e9ddeff2cf7..76378d5d5202527b54d32811ec688f00807b0382 100644
--- a/Environment/InhomogeneousMedium.h
+++ b/Environment/InhomogeneousMedium.h
@@ -41,13 +41,13 @@ namespace corsika::environment {
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
 
     corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::LineTrajectory const& pLine,
+        corsika::setup::Trajectory const& pLine,
         corsika::units::si::LengthType pTo) const override {
       return fDensityFunction.IntegrateGrammage(pLine, pTo);
     }
 
     corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::LineTrajectory const& pLine,
+        corsika::setup::Trajectory const& pLine,
         corsika::units::si::GrammageType pGrammage) const override {
       return fDensityFunction.ArclengthFromGrammage(pLine, pGrammage);
     }
diff --git a/Environment/LinearApproximationIntegrator.h b/Environment/LinearApproximationIntegrator.h
index 230d756fc9fcb7893f51b72605625945c8da8b6e..3a6c37ccb6c7b92227093322a3024a864907fe0e 100644
--- a/Environment/LinearApproximationIntegrator.h
+++ b/Environment/LinearApproximationIntegrator.h
@@ -12,7 +12,7 @@
 
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
+#include <corsika/setup/SetupTrajectory.h>
 
 namespace corsika::environment {
   template <class TDerived>
@@ -20,7 +20,7 @@ namespace corsika::environment {
     auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
 
   public:
-    auto IntegrateGrammage(corsika::geometry::LineTrajectory const& line,
+    auto IntegrateGrammage(corsika::setup::Trajectory const& line,
                            corsika::units::si::LengthType length) const {
       auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0));
       auto const c1 = GetImplementation().fRho.FirstDerivative(line.GetPosition(0),
@@ -28,7 +28,7 @@ namespace corsika::environment {
       return (c0 + 0.5 * c1 * length) * length;
     }
 
-    auto ArclengthFromGrammage(corsika::geometry::LineTrajectory const& line,
+    auto ArclengthFromGrammage(corsika::setup::Trajectory const& line,
                                corsika::units::si::GrammageType grammage) const {
       auto const c0 = GetImplementation().fRho(line.GetPosition(0));
       auto const c1 = GetImplementation().fRho.FirstDerivative(line.GetPosition(0),
@@ -37,7 +37,7 @@ namespace corsika::environment {
       return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
     }
 
-    auto MaximumLength(corsika::geometry::LineTrajectory const& line,
+    auto MaximumLength(corsika::setup::Trajectory const& line,
                        [[maybe_unused]] double relError) const {
       using namespace corsika::units::si;
       [[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
diff --git a/Environment/ShowerAxis.cc b/Environment/ShowerAxis.cc
index 3cc5779e08a6c411a39a7d5aeec0831f5c1e4f2b..76fdc66724ba57850dd404732f0b7940984176fe 100644
--- a/Environment/ShowerAxis.cc
+++ b/Environment/ShowerAxis.cc
@@ -16,14 +16,17 @@ using namespace corsika::units::si;
 using namespace corsika;
 
 GrammageType ShowerAxis::X(LengthType l) const {
-  auto const fractionalBin = l / steplength_;
+  double const fractionalBin = l / steplength_;
   int const lower = fractionalBin; // indices of nearest X support points
-  auto const lambda = fractionalBin - lower;
-  decltype(X_.size()) const upper = lower + 1;
+  double const frac = fractionalBin - lower;
+  unsigned int const upper = lower + 1;
 
-  if (lower < 0) {
+  if (fractionalBin < 0) {
     C8LOG_ERROR("cannot extrapolate to points behind point of injection l={} m", l / 1_m);
-    throw std::runtime_error("cannot extrapolate to points behind point of injection");
+    if (throw_) {
+      throw std::runtime_error("cannot extrapolate to points behind point of injection");
+    }
+    return minimumX();
   }
 
   if (upper >= X_.size()) {
@@ -31,16 +34,19 @@ GrammageType ShowerAxis::X(LengthType l) const {
         fmt::format("shower axis too short, cannot extrapolate (l / max_length_ = {} )",
                     l / max_length_);
     C8LOG_ERROR(err);
-    throw std::runtime_error(err.c_str());
+    if (throw_) { throw std::runtime_error(err.c_str()); }
+    return maximumX();
   }
 
-  assert(0 <= lambda && lambda <= 1.);
+  C8LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}", frac,
+              fractionalBin, lower, upper);
+  assert(0 <= frac && frac <= 1.);
 
-  C8LOG_TRACE("ShowerAxis::X l={} m, lower={}, lambda={}, upper={}", l / 1_m, lower,
-              lambda, upper);
+  C8LOG_TRACE("ShowerAxis::X l={} m, lower={}, frac={}, upper={}", l / 1_m, lower, frac,
+              upper);
 
   // linear interpolation between X[lower] and X[upper]
-  return X_[upper] * lambda + X_[lower] * (1 - lambda);
+  return X_[upper] * frac + X_[lower] * (1 - frac);
 }
 
 LengthType ShowerAxis::steplength() const { return steplength_; }
diff --git a/Environment/ShowerAxis.h b/Environment/ShowerAxis.h
index bfc6f1212022535a3e727c6ce8afef1a6ee771e6..052aa0bff6ac286039b4fbe580b4093b55d3e396 100644
--- a/Environment/ShowerAxis.h
+++ b/Environment/ShowerAxis.h
@@ -45,9 +45,11 @@ namespace corsika::environment {
     template <typename TEnvModel>
     ShowerAxis(geometry::Point const& pStart,
                geometry::Vector<units::si::length_d> length,
-               environment::Environment<TEnvModel> const& env, int steps = 10'000)
+               environment::Environment<TEnvModel> const& env, bool doThrow = false,
+               int steps = 10'000)
         : pointStart_(pStart)
         , length_(length)
+        , throw_(doThrow)
         , max_length_(length_.norm())
         , steplength_(max_length_ / steps)
         , axis_normalized_(length / max_length_)
@@ -104,6 +106,7 @@ namespace corsika::environment {
   private:
     geometry::Point const pointStart_;
     geometry::Vector<units::si::length_d> const length_;
+    bool throw_ = false;
     units::si::LengthType const max_length_, steplength_;
     geometry::Vector<units::si::dimensionless_d> const axis_normalized_;
     std::vector<units::si::GrammageType> X_;
diff --git a/Environment/SlidingPlanarExponential.h b/Environment/SlidingPlanarExponential.h
index b37fff2d10cc7d520794ec89f5bfadd41986cf03..aae3ccef18ebd1abfafe48afe60963a3ed5da6f1 100644
--- a/Environment/SlidingPlanarExponential.h
+++ b/Environment/SlidingPlanarExponential.h
@@ -55,14 +55,14 @@ namespace corsika::environment {
     NuclearComposition const& GetNuclearComposition() const override { return nuclComp_; }
 
     units::si::GrammageType IntegratedGrammage(
-        geometry::LineTrajectory const& line,
+        setup::Trajectory const& line,
         units::si::LengthType l) const override {
       auto const axis = (line.GetLine().GetR0() - Base::fP0).normalized();
       return Base::IntegratedGrammage(line, l, axis);
     }
 
     units::si::LengthType ArclengthFromGrammage(
-        geometry::LineTrajectory const& line,
+        setup::Trajectory const& line,
         units::si::GrammageType grammage) const override {
       auto const axis = (line.GetLine().GetR0() - Base::fP0).normalized();
       return Base::ArclengthFromGrammage(line, grammage, axis);
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 37ba3049a47b0076bd8b6e39f283b342f82cc714..eecae62928a54b4876a164389c5d84626c93cef0 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -28,6 +28,8 @@
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/setup/SetupTrajectory.h>
+
 #include <catch2/catch.hpp>
 
 using namespace corsika::geometry;
@@ -61,8 +63,7 @@ TEST_CASE("FlatExponential") {
   SECTION("horizontal") {
     Line const line(gOrigin, Vector<SpeedType::dimension_type>(
                                  gCS, {20_cm / second, 0_m / second, 0_m / second}));
-    LineTrajectory const trajectory(line, tEnd);
-
+    setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
     CHECK((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
     CHECK((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
   }
@@ -70,7 +71,7 @@ TEST_CASE("FlatExponential") {
   SECTION("vertical") {
     Line const line(gOrigin, Vector<SpeedType::dimension_type>(
                                  gCS, {0_m / second, 0_m / second, 5_m / second}));
-    LineTrajectory const trajectory(line, tEnd);
+    setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
     LengthType const length = 2 * lambda;
     GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
 
@@ -81,7 +82,7 @@ TEST_CASE("FlatExponential") {
   SECTION("escape grammage") {
     Line const line(gOrigin, Vector<SpeedType::dimension_type>(
                                  gCS, {0_m / second, 0_m / second, -5_m / second}));
-    LineTrajectory const trajectory(line, tEnd);
+    setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
     GrammageType const escapeGrammage = rho0 * lambda;
 
@@ -93,7 +94,7 @@ TEST_CASE("FlatExponential") {
   SECTION("inclined") {
     Line const line(gOrigin, Vector<SpeedType::dimension_type>(
                                  gCS, {0_m / second, 5_m / second, 5_m / second}));
-    LineTrajectory const trajectory(line, tEnd);
+    setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
     double const cosTheta = M_SQRT1_2;
     LengthType const length = 2 * lambda;
     GrammageType const exact =
@@ -127,7 +128,7 @@ TEST_CASE("SlidingPlanarExponential") {
     Line const line({gCS, {0_m, 0_m, 1_m}},
                     Vector<SpeedType::dimension_type>(
                         gCS, {0_m / second, 0_m / second, 5_m / second}));
-    LineTrajectory const trajectory(line, tEnd);
+    setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
     CHECK(medium.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
           flat.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
@@ -166,7 +167,7 @@ TEST_CASE("InhomogeneousMedium") {
                          gCS, {20_m / second, 0_m / second, 0_m / second}));
 
   auto const tEnd = 5_s;
-  LineTrajectory const trajectory(line, tEnd);
+  setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
   Exponential const e;
   DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
@@ -292,7 +293,7 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
   auto const tEnd = 1_s;
 
   // and the associated trajectory
-  LineTrajectory const trajectory(line, tEnd);
+  setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
   // and check the integrated grammage
   CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
@@ -395,7 +396,7 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
   auto const tEnd = 1_s;
 
   // and the associated trajectory
-  LineTrajectory const trajectory(line, tEnd);
+  setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
   // and check the integrated grammage
   CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
@@ -469,7 +470,7 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
   auto const tEnd = 1_s;
 
   // and the associated trajectory
-  LineTrajectory const trajectory(line, tEnd);
+  setup::Trajectory const trajectory = setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
   // and check the integrated grammage
   CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
diff --git a/Environment/testShowerAxis.cc b/Environment/testShowerAxis.cc
index 232b734cb9f4153c73826c258d109d424b83e982..e2862176bc0864bac62238d6b06dae103bccf3b2 100644
--- a/Environment/testShowerAxis.cc
+++ b/Environment/testShowerAxis.cc
@@ -64,7 +64,9 @@ TEST_CASE("Homogeneous Density") {
   Point const injectionPos = showerCore + Vector<dimensionless_d>{cs, {0, 0, 1}} * t;
 
   environment::ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos),
-                                           *env, 20};
+                                           *env,
+                                           false, // -> do not throw exceptions
+                                           20};   // -> number of bins
 
   CHECK(showerAxis.steplength() == 500_m);
 
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index e1366cbcc8800790521599faa08d8cab1d5e26c7..bc29bb620fea26ae1225f4832395bb6834426fbd 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -111,7 +111,7 @@ namespace corsika::cascade {
           count_++;
           auto pNext = stack_.GetNextParticle();
           C8LOG_DEBUG(
-              "============== next particle : count={}, pid={}, "
+              "============== next particle : count={}, pid={} "
               ", stack entries={}"
               ", stack deleted={}",
               count_, pNext.GetPID(), stack_.getEntries(), stack_.getDeleted());
@@ -202,21 +202,20 @@ namespace corsika::cascade {
                                                                          next_interact);
 
       // determine the maximum geometric step length
-      LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step);
-      C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
+      LengthType const continuous_max_dist = process_sequence_.MaxStepLength(vParticle, step);
 
       // take minimum of geometry, interaction, decay for next step
       auto min_distance =
-          std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
+          std::min({distance_interact, distance_decay, continuous_max_dist, geomMaxLength});
 
       C8LOG_DEBUG(
           "transport particle by : {} m "
           "Medium transition after: {} m "
           "Decay after: {} m "
-          "Interaction after: {} m"
-          "Continuous limit: {} m",
+          "Interaction after: {} m "
+          "Continuous limit: {} m ",
           min_distance / 1_m, geomMaxLength / 1_m, distance_decay / 1_m,
-          distance_interact / 1_m, distance_max / 1_m);
+          distance_interact / 1_m, continuous_max_dist / 1_m);
 
       // here the particle is actually moved along the trajectory to new position:
       step.SetLength(min_distance);
@@ -248,7 +247,7 @@ namespace corsika::cascade {
 
         TStackView secondaries(vParticle);
 
-        if (min_distance != distance_max) {
+        if (min_distance < continuous_max_dist) {
           /*
             Create SecondaryView object on Stack. The data container
             remains untouched and identical, and 'projectil' is identical
@@ -261,10 +260,9 @@ namespace corsika::cascade {
 
           [[maybe_unused]] auto projectile = secondaries.GetProjectile();
 
-          if (min_distance == distance_interact) {
+          if (distance_interact < distance_decay) {
             interaction(secondaries);
           } else {
-            assert(min_distance == distance_decay);
             decay(secondaries);
             // make sure particle actually did decay if it should have done so
             if (secondaries.getSize() == 1 &&
@@ -289,20 +287,32 @@ namespace corsika::cascade {
                       fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
           return numericalNodeAfterStep == currentLogicalNode;
         };
-
-        assert(assertion()); // numerical and logical nodes don't match
-      } else {               // boundary crossing, step is limited by volume boundary
-        vParticle.SetNode(nextVol);
-        /*
-          DoBoundary may delete the particle (or not)
-
-          caveat: any changes to vParticle, or even the production
-          of new secondaries is currently not passed to ParticleCut,
-          thus, particles outside the desired phase space may be produced.
-
-          todo: this must be fixed.
-        */
-        process_sequence_.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+        assert(assertion()); // numerical and logical nodes should
+                             // match, we did not cross any volume
+                             // boundary
+
+      } else { // boundary crossing, step is limited by volume boundary
+
+	if (nextVol != currentLogicalNode) {
+	
+	  C8LOG_DEBUG("volume boundary crossing to {}", fmt::ptr(nextVol));
+
+	  if (nextVol == environment_.GetUniverse().get()) {
+	    C8LOG_DEBUG("particle left physics world, is now in unknown space -> delete");
+	    vParticle.Delete();
+	  }
+	  vParticle.SetNode(nextVol);
+	  /*
+	    DoBoundary may delete the particle (or not)
+	    
+	    caveat: any changes to vParticle, or even the production
+	    of new secondaries is currently not passed to ParticleCut,
+	    thus, particles outside the desired phase space may be produced.
+	    
+	    todo: this must be fixed.
+	  */
+	  process_sequence_.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+	}
       }
     }
 
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index 013c92c56a608aed154caa4000a7ded9ceb500a2..1580a0db36b6c8e6d973513ac27916bc27a0498e 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -47,7 +47,9 @@ namespace corsika::geometry {
 
     /**
      * \param theLine The geometric \sa Line object that represents a straight-line
-     * connection \param timeLength The time duration to traverse the straight trajectory
+     * connection 
+     *
+     * \param timeLength The time duration to traverse the straight trajectory
      * in units of \sa TimeType
      */
     LineTrajectory(Line const& theLine, corsika::units::si::TimeType timeLength)
@@ -59,9 +61,15 @@ namespace corsika::geometry {
 
     /**
      * \param theLine The geometric \sa Line object that represents a straight-line
-     * connection \param timeLength The time duration to traverse the straight trajectory
-     * in units of \sa TimeType \param timeStep Time duration to folow eventually curved
-     * trajectory in units of \sa TimesType \param initialV Initial velocity vector at
+     * connection 
+     * 
+     * \param timeLength The time duration to traverse the straight trajectory
+     * in units of \sa TimeType 
+     * 
+     * \param timeStep Time duration to folow eventually curved
+     * trajectory in units of \sa TimesType 
+     * 
+     * \param initialV Initial velocity vector at
      * start of trajectory \param finalV Final velocity vector at start of trajectory
      */
     LineTrajectory(
@@ -187,14 +195,23 @@ namespace corsika::geometry {
         , k_(k)
         , timeStep_(timeStep) {}
 
-    const Line GetLine(double u) const { return Line(GetPosition(u), GetVelocity(u)); }
+    const Line GetLine() const {
+      using namespace corsika::units::si;
+      auto D = GetPosition(1) - GetPosition(0);
+      auto d = D.norm();
+      auto v = initialVelocity_;
+      if (d>1_um) { // if trajectory is ultra-short, we do not
+		    // re-calculate velocity, just use initial
+		    // value. Otherwise, this is numerically unstable
+	v = D/d * GetVelocity(0).norm();
+      }
+      return Line(GetPosition(0), v);
+    }
     Point GetPosition(double u) const {
       Point position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
       VelocityVec velocity =
           initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
       return position + velocity * timeStep_ * u / 2;
-      //      auto steplength_true = steplength_true * (1.0 + double(direction.norm())) /
-      //      2;
     }
     VelocityVec GetVelocity(double u) const {
       return initialVelocity_ +
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index 6153c63567ff3efd56c207a04a1d57efd45c116f..dd8b797436a7bb1c892b75ca742ff6a9e0557ade 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -30,6 +30,7 @@ ObservationPlane::ObservationPlane(
 
 corsika::process::EProcessReturn ObservationPlane::DoContinuous(
     setup::Stack::ParticleType& particle, setup::Trajectory const& trajectory) {
+  
   TimeType const timeOfIntersection =
       (plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
       trajectory.GetLine().GetV0().dot(plane_.GetNormal());
diff --git a/Processes/ObservationPlane/testObservationPlane.cc b/Processes/ObservationPlane/testObservationPlane.cc
index ac5573e0b2df8552bd09ed246793e9452d75712b..ddaaaad96e4795df9cb2695fde09a3de24e0c069 100644
--- a/Processes/ObservationPlane/testObservationPlane.cc
+++ b/Processes/ObservationPlane/testObservationPlane.cc
@@ -19,9 +19,8 @@
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/units/PhysicalUnits.h>
 
-//#include <corsika/setup/SetupEnvironment.h>
 #include <corsika/setup/SetupStack.h>
-//#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/setup/SetupTrajectory.h>
 
 using namespace corsika::units::si;
 using namespace corsika::process::observation_plane;
@@ -51,7 +50,8 @@ TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
   Vector<units::si::SpeedType::dimension_type> vec(cs, 0_m / second, 0_m / second,
                                                    -units::constants::c);
   Line line(start, vec);
-  LineTrajectory track(line, 12_m / units::constants::c);
+  setup::Trajectory track =
+      setup::testing::make_track<setup::Trajectory>(line, 12_m / units::constants::c);
 
   particle.SetPosition(Point(cs, {1_m, 1_m, 10_m})); // moving already along -z
 
diff --git a/Processes/ParticleCut/testParticleCut.cc b/Processes/ParticleCut/testParticleCut.cc
index 95e217bc73017b0c39e84b415fa722bd82d287a7..c2343f76373acc2712f553e2803dd4ea10e02644 100644
--- a/Processes/ParticleCut/testParticleCut.cc
+++ b/Processes/ParticleCut/testParticleCut.cc
@@ -170,11 +170,11 @@ TEST_CASE("ParticleCut", "[processes]") {
     CHECK(cut.GetCutEnergy() == 0_GeV);
   }
 
-  corsika::setup::Trajectory const track{
+  corsika::setup::Trajectory const track = setup::testing::make_track<setup::Trajectory>(
       geometry::Line{point0,
                      geometry::Vector<units::si::SpeedType::dimension_type>{
                          rootCS, {0_m / second, 0_m / second, -units::constants::c}}},
-      12_m / units::constants::c};
+      12_m / units::constants::c);
 
   SECTION("cut on DoContinous, just invisibles") {
 
diff --git a/Processes/Tracking/Intersect.hpp b/Processes/Tracking/Intersect.hpp
index 277c5f6336478ce8dd2cbf0c84290cd72be8b6a0..717ad6ca99647ea1b68e5dce5ad645db0c67449d 100644
--- a/Processes/Tracking/Intersect.hpp
+++ b/Processes/Tracking/Intersect.hpp
@@ -67,15 +67,18 @@ namespace corsika::process::tracking {
       // thus, the last entry is always the exit point
       const Intersections time_intersections_curr =
           TDerived::Intersect(particle, volumeNode);
-      C8LOG_TRACE("curr node {}, parent node {} ", fmt::ptr(&volumeNode),
-                  fmt::ptr(volumeNode.GetParent()));
-      C8LOG_DEBUG("intersection times with currentLogicalVolumeNode: {} s and {} s",
-                  time_intersections_curr.getEntry() / 1_s,
-                  time_intersections_curr.getExit() / 1_s);
-      if (time_intersections_curr.getExit() <= minTime) {
-        minTime =
-            time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
-        minNode = volumeNode.GetParent();
+      C8LOG_TRACE("curr node {}, parent node {}, hasIntersections={} ",
+                  fmt::ptr(&volumeNode), fmt::ptr(volumeNode.GetParent()),
+                  time_intersections_curr.hasIntersections());
+      if (time_intersections_curr.hasIntersections()) {
+        C8LOG_DEBUG("intersection times with currentLogicalVolumeNode: {} s and {} s",
+                    time_intersections_curr.getEntry() / 1_s,
+                    time_intersections_curr.getExit() / 1_s);
+        if (time_intersections_curr.getExit() <= minTime) {
+          minTime =
+              time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
+          minNode = volumeNode.GetParent();
+        }
       }
 
       // where do we collide with any of the next-tree-level volumes
@@ -83,6 +86,7 @@ namespace corsika::process::tracking {
       for (const auto& node : volumeNode.GetChildNodes()) {
 
         const Intersections time_intersections = TDerived::Intersect(particle, *node);
+        if (!time_intersections.hasIntersections()) { continue; }
         C8LOG_DEBUG("intersection times with child volume {} : enter {} s, exit {} s",
                     fmt::ptr(node), time_intersections.getEntry() / 1_s,
                     time_intersections.getExit() / 1_s);
@@ -93,8 +97,13 @@ namespace corsika::process::tracking {
                     t_entry <= minTime);
         // note, theoretically t can even be smaller than 0 since we
         // KNOW we can't yet be in this volume yet, so we HAVE TO
-        // enter it IF exit point is not also in the "past"!
-        if (t_exit > 0_s && t_entry <= minTime) { // enter volumen child here
+        // enter it IF exit point is not also in the "past", AND if
+        // extry point is [[much much]] closer than exit point
+        // (because we might have already numerically "exited" it)!
+        if (t_exit > 0_s && t_entry <= minTime &&
+            -t_entry < t_exit) { // protection agains numerical problem, when we already
+                                 // _exited_ before
+                                 // enter chile volume here
           minTime = t_entry;
           minNode = node.get();
         }
@@ -105,6 +114,7 @@ namespace corsika::process::tracking {
       for (node_type* node : volumeNode.GetExcludedNodes()) {
 
         const Intersections time_intersections = TDerived::Intersect(particle, *node);
+        if (!time_intersections.hasIntersections()) { continue; }
         C8LOG_DEBUG("intersection times with exclusion volume {} : enter {} s, exit {} s",
                     fmt::ptr(node), time_intersections.getEntry() / 1_s,
                     time_intersections.getExit() / 1_s);
diff --git a/Processes/TrackingLeapFrogCurved/Tracking.h b/Processes/TrackingLeapFrogCurved/Tracking.h
index 2f3b770c1608ff66ecbe11770800e98ed56feb0f..9993f7b70084b167e21c539e86a37b0a8d570b54 100644
--- a/Processes/TrackingLeapFrogCurved/Tracking.h
+++ b/Processes/TrackingLeapFrogCurved/Tracking.h
@@ -156,6 +156,8 @@ namespace corsika::process {
           return tracking_line::Tracking::Intersect(particle, sphere, medium);
         }
 
+        bool const numericallyInside = sphere.Contains(particle.GetPosition());
+
         const geometry::Vector<SpeedType::dimension_type> velocity =
             particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
         const auto absVelocity = velocity.norm();
@@ -181,27 +183,58 @@ namespace corsika::process {
                           k * 1_m * 1_m * 1_m * 1_m);
         std::complex<double>* solutions = solve_quartic(0, a, b, c);
         LengthType d_enter, d_exit;
-        int first = 0;
+        int first = 0, first_entry = 0, first_exit = 0;
         for (int i = 0; i < 4; i++) {
           if (solutions[i].imag() == 0) {
-            LengthType time = solutions[i].real() * 1_m;
-            C8LOG_TRACE("Solutions for current Volume: {} ", time);
-            if (first == 0) {
-              d_enter = time;
-            } else {
-              if (time < d_enter) {
-                d_exit = d_enter;
-                d_enter = time;
+            LengthType const dist = solutions[i].real() * 1_m;
+            C8LOG_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 (first_entry == 0) {
+                  d_enter = dist;
+                } else {
+                  d_enter = std::max(d_enter, dist); // closest negative to zero (-1e-4) m
+                }
+                first_entry++;
+
+              } else { // thus, dist >= -0.0001_m
+
+                if (first_exit == 0) {
+                  d_exit = dist;
+                } else {
+                  d_exit = std::min(d_exit, dist); // closest positive to zero (-1e-4) m
+                }
+                first_exit++;
+              }
+              first = int(first_exit > 0) + int(first_entry > 0);
+
+            } else { // thus, numericallyInside == false
+
+              // 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
+                // into next logical volume
+                continue;
+              }
+              if (first == 0) {
+                d_enter = dist;
               } else {
-                d_exit = time;
+                if (dist < d_enter) {
+                  d_exit = d_enter;
+                  d_enter = dist;
+                } else {
+                  d_exit = dist;
+                }
               }
+              first++;
             }
-            first++;
-          }
+          } // loop over solutions
         }
         delete[] solutions;
 
-        if (first != 2) {
+        if (first != 2) { // entry and exit points found
           C8LOG_DEBUG("no intersection! count={}", first);
           return geometry::Intersections();
         }
@@ -219,7 +252,7 @@ namespace corsika::process {
         throw std::runtime_error(
             "The Volume type provided is not supported in Intersect(particle, node)");
       }
-    };
+    }; // namespace tracking_leapfrog_curved
 
   } // namespace tracking_leapfrog_curved
 
diff --git a/Processes/TrackingLeapFrogStraight/Tracking.h b/Processes/TrackingLeapFrogStraight/Tracking.h
index d9e7c57de709866bf410a8050016d70bdf62ebd3..b88a6fc46883115866171d9acf7e3b584700c889 100644
--- a/Processes/TrackingLeapFrogStraight/Tracking.h
+++ b/Processes/TrackingLeapFrogStraight/Tracking.h
@@ -18,7 +18,6 @@
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/process/tracking_line/Tracking.h>
-#include <corsika/utl/quartic.h>
 
 #include <type_traits>
 #include <utility>
@@ -89,12 +88,12 @@ namespace corsika::process {
         auto const momentumVerticalMag =
             particle.GetMomentum() -
             particle.GetMomentum().parallelProjectionOnto(magneticfield);
+        bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
         LengthType const gyroradius =
-            (chargeNumber == 0 || magnitudeB == 0_T
-                 ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
-                 : momentumVerticalMag.norm() * 1_V /
-                       (corsika::units::constants::c * abs(chargeNumber) * magnitudeB *
-                        1_eV));
+            (no_deflection ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
+                           : momentumVerticalMag.norm() * 1_V /
+                                 (corsika::units::constants::c * abs(chargeNumber) *
+                                  magnitudeB * 1_eV));
         const double maxRadians = 0.01;
         const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
         C8LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
@@ -104,7 +103,7 @@ namespace corsika::process {
         const auto absMomentum = initialMomentum.norm();
         const auto absVelocity = initialVelocity.norm();
         const geometry::Vector<dimensionless_d> direction = initialVelocity.normalized();
-        ;
+
         // check if particle is moving at all
         if (absVelocity * 1_s == 0_m) {
           return std::make_tuple(
@@ -123,8 +122,15 @@ namespace corsika::process {
                     initialTrack.GetPosition(0).GetCoordinates(),
                     initialTrack.GetPosition(1).GetCoordinates(), initialTrackLength);
 
+        // if particle is non-deflectable, we are done:
+        if (no_deflection) {
+          C8LOG_DEBUG("no deflection. tracking finished");
+          return std::make_tuple(initialTrack, initialTrackNextVolume);
+        }
+
         // avoid any intersections within first halve steplength
-        LengthType firstHalveSteplength = std::min(steplimit, initialTrackLength) / 2;
+        LengthType const firstHalveSteplength =
+            std::min(steplimit, initialTrackLength) / 2;
 
         C8LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
                     firstHalveSteplength, steplimit, initialTrackLength);
@@ -135,9 +141,12 @@ namespace corsika::process {
         const auto new_direction =
             direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
         const auto new_direction_norm = new_direction.norm(); // by design this is >1
-        C8LOG_DEBUG("position_mid={}, new_direction={}, new_direction_norm={}",
-                    position_mid.GetCoordinates(), new_direction.GetComponents(),
-                    new_direction_norm);
+        C8LOG_DEBUG(
+            "position_mid={}, new_direction={}, (new_direction_norm)={}, deflection={}",
+            position_mid.GetCoordinates(), new_direction.GetComponents(),
+            new_direction_norm,
+            acos(std::min(1.0, direction.dot(new_direction) / new_direction_norm)) * 180 /
+                M_PI);
 
         // check, where the second halve-step direction has geometric intersections
         particle.SetPosition(position_mid);
@@ -146,18 +155,38 @@ namespace corsika::process {
             tracking_line::Tracking::GetTrack(particle);
         particle.SetPosition(initialPosition); // this is not nice...
         particle.SetMomentum(initialMomentum); // this is not nice...
-        const auto finalTrackLength = finalTrack.GetLength(1);
 
-        C8LOG_DEBUG("finalTrack(0)={}, finalTrack(1)={}, finalTrackLength={}",
-                    finalTrack.GetPosition(0).GetCoordinates(),
-                    finalTrack.GetPosition(1).GetCoordinates(), finalTrackLength);
-
-        const LengthType secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
-        const LengthType secondHalveStepLength =
+        LengthType const finalTrackLength = finalTrack.GetLength(1);
+        LengthType const secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
+
+        // check if volume transition is obvious, OR
+        // for numerical reasons, particles slighly bend "away" from a
+        // volume boundary have a very hard time to cross the border,
+        // thus, if secondLeapFrogLength is just slighly shorter (1e-4m) than
+        // finalTrackLength we better just [extend the
+        // secondLeapFrogLength slightly and] force the volume
+        // crossing:
+        bool const switch_volume = finalTrackLength - 0.0001_m <= secondLeapFrogLength;
+        LengthType const secondHalveStepLength =
             std::min(secondLeapFrogLength, finalTrackLength);
 
+        C8LOG_DEBUG(
+            "finalTrack(0)={}, finalTrack(1)={}, finalTrackLength={}, "
+            "secondLeapFrogLength={}, secondHalveStepLength={}, "
+            "secondLeapFrogLength-finalTrackLength={}, "
+            "secondHalveStepLength-finalTrackLength={}, "
+            "nextVol={}, transition={}",
+            finalTrack.GetPosition(0).GetCoordinates(),
+            finalTrack.GetPosition(1).GetCoordinates(), finalTrackLength,
+            secondLeapFrogLength, secondHalveStepLength,
+            secondLeapFrogLength - finalTrackLength,
+            secondHalveStepLength - finalTrackLength, fmt::ptr(finalTrackNextVolume),
+            switch_volume);
+
         // perform the second halve-step
-        const Point finalPosition = position_mid + new_direction * secondHalveStepLength;
+        auto const new_direction_normalized = new_direction.normalized();
+        const Point finalPosition =
+            position_mid + new_direction_normalized * secondHalveStepLength;
 
         const LengthType totalStep = firstHalveSteplength + secondHalveStepLength;
         const auto delta_pos = finalPosition - initialPosition;
@@ -171,10 +200,8 @@ namespace corsika::process {
                 distance / absVelocity,  // straight distance
                 totalStep / absVelocity, // bend distance
                 initialVelocity,
-                new_direction.normalized() * absVelocity), // trajectory
-            (finalTrackLength > secondLeapFrogLength
-                 ? volumeNode
-                 : finalTrackNextVolume)); // next step volume
+                new_direction_normalized * absVelocity), // trajectory
+            (switch_volume ? finalTrackNextVolume : volumeNode));
       }
     };
 
diff --git a/Processes/TrackingLine/Tracking.h b/Processes/TrackingLine/Tracking.h
index 0d099ca97ea18acaf2a9a9341b0dc5acf066d6f9..d9690ec0d9c42e05a95fb8e41745d8242c2450b9 100644
--- a/Processes/TrackingLine/Tracking.h
+++ b/Processes/TrackingLine/Tracking.h
@@ -11,13 +11,12 @@
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Plane.h>
 #include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Trajectory.h>
 #include <corsika/geometry/Vector.h>
 #include <corsika/geometry/Intersections.hpp>
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/logging/Logging.h>
-#include <corsika/setup/SetupEnvironment.h>
 #include <corsika/process/tracking/Intersect.hpp>
+#include <corsika/geometry/Trajectory.h>
 
 #include <type_traits>
 #include <utility>
diff --git a/Setup/SetupTrajectory.h b/Setup/SetupTrajectory.h
index a441bfcab8002d2917e73a40d32f4a452c2f5dbd..a8b86e030fbb7eabc033f226b3693f13c1dea592 100644
--- a/Setup/SetupTrajectory.h
+++ b/Setup/SetupTrajectory.h
@@ -12,13 +12,49 @@
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Trajectory.h>
 
-#include <corsika/units/PhysicalUnits.h>
+#include <corsika/process/tracking_line/Tracking.h>
+#include <corsika/process/tracking_leapfrog_curved/Tracking.h>
+#include <corsika/process/tracking_leapfrog_straight/Tracking.h>
 
-// #include <variant>
+#include <corsika/units/PhysicalUnits.h>
 
 namespace corsika::setup {
 
+  // Note: Tracking and Trajectory must fit together !
+  typedef corsika::process::tracking_leapfrog_curved::Tracking Tracking;
+  // tracking_leapfrog_straight::Tracking tracking;
+
   /// definition of Trajectory base class, to be used in tracking and cascades
-  typedef corsika::geometry::LineTrajectory Trajectory;
+  // typedef corsika::geometry::LineTrajectory Trajectory;
+  typedef corsika::geometry::LeapFrogTrajectory Trajectory;
+
+  namespace testing {
+
+    template <typename TTrack>
+    TTrack make_track(const corsika::geometry::Line& line,
+                      const corsika::units::si::TimeType tEnd);
+
+    template <>
+    inline corsika::geometry::LineTrajectory
+    make_track<corsika::geometry::LineTrajectory>(
+        const corsika::geometry::Line& line, const corsika::units::si::TimeType tEnd) {
+      return corsika::geometry::LineTrajectory(line, tEnd);
+    }
+
+    template <>
+    inline corsika::geometry::LeapFrogTrajectory
+    make_track<corsika::geometry::LeapFrogTrajectory>(
+        const corsika::geometry::Line& line, const corsika::units::si::TimeType tEnd) {
+      using namespace corsika::units::si;
+      typedef corsika::geometry::Vector<magnetic_flux_density_d> MagneticFieldVector;
+
+      auto const k = square(0_m) / (square(1_s) * 1_V);
+      return corsika::geometry::LeapFrogTrajectory(
+          line.GetR0(), line.GetV0(),
+          MagneticFieldVector{line.GetR0().GetCoordinateSystem(), 0_T, 0_T, 0_T}, k,
+          tEnd);
+    }
+
+  } // namespace testing
 
 } // namespace corsika::setup