From dbf6e99568a9fced83758456c61ad4bff7defcc0 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 23 Nov 2020 09:09:30 +0100
Subject: [PATCH] made leap-frog more stabel and robust, also in extrem
magnetic field cases
---
Documentation/Examples/boundary_example.cc | 3 +-
Documentation/Examples/cascade_example.cc | 3 +-
.../Examples/cascade_proton_example.cc | 3 +-
Documentation/Examples/em_shower.cc | 3 +-
Documentation/Examples/hybrid_MC.cc | 3 +-
Documentation/Examples/vertical_EAS.cc | 13 ++--
Environment/BaseExponential.h | 6 +-
Environment/DensityFunction.h | 1 -
Environment/FlatExponential.h | 7 +-
Environment/HomogeneousMedium.h | 6 +-
Environment/IMediumModel.h | 7 +-
Environment/InhomogeneousMedium.h | 4 +-
Environment/LinearApproximationIntegrator.h | 8 +-
Environment/ShowerAxis.cc | 26 ++++---
Environment/ShowerAxis.h | 5 +-
Environment/SlidingPlanarExponential.h | 4 +-
Environment/testEnvironment.cc | 21 +++---
Environment/testShowerAxis.cc | 4 +-
Framework/Cascade/Cascade.h | 58 +++++++++------
Framework/Geometry/Trajectory.h | 31 ++++++--
.../ObservationPlane/ObservationPlane.cc | 1 +
.../ObservationPlane/testObservationPlane.cc | 6 +-
Processes/ParticleCut/testParticleCut.cc | 4 +-
Processes/Tracking/Intersect.hpp | 32 +++++---
Processes/TrackingLeapFrogCurved/Tracking.h | 61 ++++++++++++----
Processes/TrackingLeapFrogStraight/Tracking.h | 73 +++++++++++++------
Processes/TrackingLine/Tracking.h | 3 +-
Setup/SetupTrajectory.h | 42 ++++++++++-
28 files changed, 291 insertions(+), 147 deletions(-)
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 21cdbe353..c40823c14 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 2707e3dde..755c7c7c1 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 11b3dac77..909cad004 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 d834afe77..31298aa28 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 a9f7432ae..75ed47a19 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 235c5b3f5..834989466 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 94b454bf6..bf8c16d96 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 5f61e26d7..3ac33c3f2 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 31f330994..8aeafbce4 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 d630d8a6b..f39cc7061 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 388f79c43..7a3b95171 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 87002d6b1..76378d5d5 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 230d756fc..3a6c37ccb 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 3cc5779e0..76fdc6672 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 bfc6f1212..052aa0bff 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 b37fff2d1..aae3ccef1 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 37ba3049a..eecae6292 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 232b734cb..e2862176b 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 e1366cbcc..bc29bb620 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 013c92c56..1580a0db3 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 6153c6356..dd8b79743 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 ac5573e0b..ddaaaad96 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 95e217bc7..c2343f763 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 277c5f633..717ad6ca9 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 2f3b770c1..9993f7b70 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 d9e7c57de..b88a6fc46 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 0d099ca97..d9690ec0d 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 a441bfcab..a8b86e030 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
--
GitLab