diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5cd58a179c7a15e027af21b1bad2b652b2178918..2e9df74375dcef20916ccdf997378bd00b6089f7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -115,7 +115,7 @@ endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
#+++++++++++++++++++++++++++++
# Setup external dependencies.
#
-include (conan.cmake)
+include (conan) # self-provided in 'cmake' directory
#
# download and build all dependencies
conan_cmake_run (CONANFILE conanfile.txt
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 2d68f8c1019849314ac9c082bf8c34b4d08e8409..0a01b5202f6469f70e1bcbe3afc82a85929a1611 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -27,16 +27,6 @@
namespace corsika {
- template <typename TTracking, typename TProcessList, typename TStack,
- typename TStackView>
- void Cascade<TTracking, TProcessList, TStack, TStackView>::setNodes() {
- std::for_each(stack_.begin(), stack_.end(), [&](auto& p) {
- auto const* numericalNode =
- environment_.getUniverse()->getContainingNode(p.getPosition());
- p.setNode(numericalNode);
- });
- }
-
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
void Cascade<TTracking, TProcessList, TStack, TStackView>::run() {
@@ -68,6 +58,7 @@ namespace corsika {
typename TStackView>
void Cascade<TTracking, TProcessList, TStack, TStackView>::forceInteraction() {
CORSIKA_LOG_TRACE("forced interaction!");
+ setNodes();
auto vParticle = stack_.getNextParticle();
TStackView secondaries(vParticle);
interaction(secondaries);
@@ -79,16 +70,12 @@ namespace corsika {
typename TStackView>
void Cascade<TTracking, TProcessList, TStack, TStackView>::step(Particle& vParticle) {
- // determine geometric tracking
- auto [step, geomMaxLength, nextVol] = tracking_.getTrack(vParticle);
- [[maybe_unused]] auto const& dummy_nextVol = nextVol;
-
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
sequence_.getInverseInteractionLength(vParticle);
// sample random exponential step length in grammage
- corsika::ExponentialDistribution expDist(1 / total_inv_lambda);
+ ExponentialDistribution expDist(1 / total_inv_lambda);
GrammageType const next_interact = expDist(rng_);
CORSIKA_LOG_DEBUG(
@@ -101,22 +88,15 @@ namespace corsika {
// assert that particle stays outside void Universe if it has no
// model properties set
- assert(currentLogicalNode != &*environment_.getUniverse() ||
- environment_.getUniverse()->hasModelProperties());
-
- // convert next_step from grammage to length
- LengthType const distance_interact =
- currentLogicalNode->getModelProperties().getArclengthFromGrammage(step,
- next_interact);
-
- // determine the maximum geometric step length
- LengthType const distance_max = sequence_.getMaxStepLength(vParticle, step);
+ assert((currentLogicalNode != &*environment_.getUniverse() ||
+ environment_.getUniverse()->hasModelProperties()) &&
+ "FATAL: The environment model has no valid properties set!");
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime = sequence_.getInverseLifetime(vParticle);
// sample random exponential decay time
- corsika::ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
+ ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
TimeType const next_decay = expDistDecay(rng_);
CORSIKA_LOG_DEBUG(
@@ -128,19 +108,39 @@ namespace corsika {
LengthType const distance_decay = next_decay * vParticle.getMomentum().getNorm() /
vParticle.getEnergy() * constants::c;
+ // determine geometric tracking
+ auto [step, nextVol] = tracking_.getTrack(vParticle);
+ auto geomMaxLength = step.getLength(1);
+
+ // convert next_step from grammage to length
+ LengthType const distance_interact =
+ currentLogicalNode->getModelProperties().getArclengthFromGrammage(step,
+ next_interact);
+
+ // determine the maximum geometric step length
+ LengthType const continuous_max_dist =
+ process_sequence_.getMaxStepLength(vParticle, step);
+
// take minimum of geometry, interaction, decay for next step
auto const min_distance =
- std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
+ std::min({distance_interact, distance_decay, continuous_max_dist, geomMaxLength});
- CORSIKA_LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
+ CORSIKA_LOG_DEBUG(
+ "transport particle by : {} m "
+ "Medium transition after: {} m "
+ "Decay after: {} m "
+ "Interaction after: {} m "
+ "Continuous limit: {} m ",
+ min_distance / 1_m, geomMaxLength / 1_m, distance_decay / 1_m,
+ distance_interact / 1_m, continuous_max_dist / 1_m);
// here the particle is actually moved along the trajectory to new position:
// std::visit(setup::ParticleUpdate<particle_type>{vParticle}, step);
- vParticle.setPosition(step.getPositionFromArclength(min_distance));
- // .... also update time, momentum, direction, ...
- vParticle.setTime(vParticle.getTime() + min_distance / constants::c);
-
- step.getLimitEndTo(min_distance);
+ step.setLength(min_distance);
+ vParticle.setPosition(step.getPosition(1));
+ // assumption: tracking does not change absolute momentum:
+ vParticle.setMomentum(step.getDirection(1) * vParticle.getMomentum().norm());
+ vParticle.setTime(vParticle.getTime() + step.getDuration());
// apply all continuous processes on particle + track
if (sequence_.doContinuous(vParticle, step) == ProcessReturn::ParticleAbsorbed) {
@@ -162,7 +162,7 @@ namespace corsika {
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
@@ -175,10 +175,9 @@ namespace corsika {
[[maybe_unused]] auto projectile = secondaries.getProjectile();
- if (min_distance == distance_interact) {
+ if (distance_interat < 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 &&
@@ -192,6 +191,7 @@ namespace corsika {
} else { // step-length limitation within volume
CORSIKA_LOG_DEBUG("step-length limitation");
+ // no further physics happens here. just proceed to next step.
}
[[maybe_unused]] auto const assertion = [&] {
@@ -203,68 +203,91 @@ namespace corsika {
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)
+ assert(assertion()); // numerical and logical nodes should match, since
+ // we did not cross any volume boundary
- 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.
+ } else { // boundary crossing, step is limited by volume boundary
- \todo: this must be fixed.
- */
+ if (nextVol != currentLogicalNode) {
- sequence_.doBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+ CORSIKA_LOG_DEBUG("volume boundary crossing to {}", fmt::ptr(nextVol));
+
+ if (nextVol == environment_.getUniverse().get()) {
+ CORSIKA_LOG_DEBUG(
+ "particle left physics world, is now in unknown space -> delete");
+ vParticle.erase();
+ }
+ 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.
+ */
+
+ sequence_.doBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+ }
}
- }
- template <typename TTracking, typename TProcessList, typename TStack,
- typename TStackView>
- ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::decay(
- TStackView& view) {
- CORSIKA_LOG_DEBUG("decay");
- InverseTimeType const actual_decay_time = sequence_.getInverseLifetime(view.parent());
+ template <typename TTracking, typename TProcessList, typename TStack,
+ typename TStackView>
+ ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::decay(TStackView &
+ view) {
+ CORSIKA_LOG_DEBUG("decay");
+ InverseTimeType const actual_decay_time =
+ sequence_.getInverseLifetime(view.parent());
- corsika::UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
- const auto sample_process = uniDist(rng_);
+ UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
+ const auto sample_process = uniDist(rng_);
- auto const returnCode = sequence_.selectDecay(view, sample_process);
- if (returnCode != ProcessReturn::Decayed) {
- CORSIKA_LOG_WARN("Particle did not decay!");
+ auto const returnCode = sequence_.selectDecay(view, sample_process);
+ if (returnCode != ProcessReturn::Decayed) {
+ CORSIKA_LOG_WARN("Particle did not decay!");
+ }
+ setEventType(view, history::EventType::Decay);
+ return returnCode;
}
- setEventType(view, history::EventType::Decay);
- return returnCode;
- }
- template <typename TTracking, typename TProcessList, typename TStack,
- typename TStackView>
- ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::interaction(
- TStackView& view) {
- CORSIKA_LOG_DEBUG("collide");
+ template <typename TTracking, typename TProcessList, typename TStack,
+ typename TStackView>
+ ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::interaction(
+ TStackView & view) {
+ CORSIKA_LOG_DEBUG("collide");
- InverseGrammageType const current_inv_length =
- sequence_.getInverseInteractionLength(view.parent());
+ InverseGrammageType const current_inv_length =
+ sequence_.getInverseInteractionLength(view.parent());
- corsika::UniformRealDistribution<InverseGrammageType> uniDist(current_inv_length);
+ UniformRealDistribution<InverseGrammageType> uniDist(current_inv_length);
- const auto sample_process = uniDist(rng_);
- auto const returnCode = sequence_.selectInteraction(view, sample_process);
- if (returnCode != ProcessReturn::Interacted) {
- CORSIKA_LOG_WARN("Particle did not interace!");
+ const auto sample_process = uniDist(rng_);
+ auto const returnCode = sequence_.selectInteraction(view, sample_process);
+ if (returnCode != ProcessReturn::Interacted) {
+ CORSIKA_LOG_WARN("Particle did not interace!");
+ }
+ setEventType(view, history::EventType::Interaction);
+ return returnCode;
}
- setEventType(view, history::EventType::Interaction);
- return returnCode;
- }
- template <typename TTracking, typename TProcessList, typename TStack,
- typename TStackView>
- void Cascade<TTracking, TProcessList, TStack, TStackView>::setEventType(
- TStackView& view, [[maybe_unused]] history::EventType eventType) {
- if constexpr (TStackView::has_event) {
- for (auto&& sec : view) { sec.getEvent()->setEventType(eventType); }
+ template <typename TTracking, typename TProcessList, typename TStack,
+ typename TStackView>
+ void Cascade<TTracking, TProcessList, TStack, TStackView>::setNodes() {
+ std::for_each(stack_.begin(), stack_.end(), [&](auto& p) {
+ auto const* numericalNode =
+ environment_.getUniverse()->getContainingNode(p.getPosition());
+ p.setNode(numericalNode);
+ });
+ }
+
+ template <typename TTracking, typename TProcessList, typename TStack,
+ typename TStackView>
+ void Cascade<TTracking, TProcessList, TStack, TStackView>::setEventType(
+ TStackView & view, [[maybe_unused]] history::EventType eventType) {
+ if constexpr (TStackView::has_event) {
+ for (auto&& sec : view) { sec.getEvent()->setEventType(eventType); }
+ }
}
- }
-} // namespace corsika
+ } // namespace corsika
diff --git a/corsika/detail/framework/geometry/Helix.inl b/corsika/detail/framework/geometry/Helix.inl
index c20a8f21fccfed478815f21f5bf354ad3ad2a93a..21b1e5a61b3d37298939aea2ceba0ba9b04b8f8f 100644
--- a/corsika/detail/framework/geometry/Helix.inl
+++ b/corsika/detail/framework/geometry/Helix.inl
@@ -22,16 +22,19 @@ namespace corsika {
(vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
omegaC_;
}
+ inline VelocityVector Helix::getVelocity(TimeType const t) const {
+ return vPar_ + (vPerp_ * (cos(omegaC_ * t) - 1) + uPerp_ * sin(omegaC_ * t));
+ }
- Point Helix::getPositionFromArclength(LengthType const l) const {
+ inline Point Helix::getPositionFromArclength(LengthType const l) const {
return getPosition(getTimeFromArclength(l));
}
-
- LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
+
+ inline LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
return (vPar_ + vPerp_).getNorm() * (t2 - t1);
}
- TimeType Helix::getTimeFromArclength(LengthType const l) const {
+ inline TimeType Helix::getTimeFromArclength(LengthType const l) const {
return l / (vPar_ + vPerp_).getNorm();
}
diff --git a/corsika/detail/framework/geometry/Line.inl b/corsika/detail/framework/geometry/Line.inl
index e3a5ef6ee65db71e8c446e0704f4e2f25915b5e8..aa99662c5af00eda0a2394c21a5d92baa6def84a 100644
--- a/corsika/detail/framework/geometry/Line.inl
+++ b/corsika/detail/framework/geometry/Line.inl
@@ -14,22 +14,24 @@
namespace corsika {
- Point Line::getPosition(TimeType const t) const { return start_point_ + velocity_ * t; }
+ inline Point Line::getPosition(TimeType const t) const { return start_point_ + velocity_ * t; }
- Point Line::getPositionFromArclength(LengthType const l) const {
+ inline VelocityVector const& Line::getVelocity(TimeType const) const { return velocity_; }
+
+ inline Point Line::getPositionFromArclength(LengthType const l) const {
return start_point_ + velocity_.normalized() * l;
}
- LengthType Line::getArcLength(TimeType const t1, TimeType const t2) const {
+ inline LengthType Line::getArcLength(TimeType const t1, TimeType const t2) const {
return velocity_.getNorm() * (t2 - t1);
}
- TimeType Line::getTimeFromArclength(LengthType const t) const {
+ inline TimeType Line::getTimeFromArclength(LengthType const t) const {
return t / velocity_.getNorm();
}
- Point const& Line::getStartPoint() const { return start_point_; }
+ inline Point const& Line::getStartPoint() const { return start_point_; }
- Line::VelocityVec const& Line::getVelocity() const { return velocity_; }
+ inline VelocityVector const& Line::getVelocity() const { return velocity_; }
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Sphere.inl b/corsika/detail/framework/geometry/Sphere.inl
index 9a6688eb506fab5159bdc009d2b72512544888f4..a793984dd1b2bb5a46ba386887062875221fa136 100644
--- a/corsika/detail/framework/geometry/Sphere.inl
+++ b/corsika/detail/framework/geometry/Sphere.inl
@@ -13,12 +13,16 @@
namespace corsika {
- bool Sphere::isInside(Point const& p) const {
+ inline bool Sphere::isInside(Point const& p) const {
return radius_ * radius_ > (center_ - p).getSquaredNorm();
}
- Point const& Sphere::getCenter() const { return center_; }
+ inline Point const& Sphere::getCenter() const { return center_; }
- LengthType Sphere::getRadius() const { return radius_; }
+ inline void Sphere::setCenter(Point const& p) { center_ = p; }
+
+ inline LengthType Sphere::getRadius() const { return radius_; }
+
+ inline void Sphere::setRadius(LengthType const r) { radius_ = r; }
} // namespace corsika
diff --git a/corsika/detail/framework/geometry/Trajectory.inl b/corsika/detail/framework/geometry/Trajectory.inl
index fd8aa4f5c0d3387b5d3fb0a40dbbcba5e66423a0..ab7779c83a701e82a5f05e8dc9f9b3da9fe14661 100644
--- a/corsika/detail/framework/geometry/Trajectory.inl
+++ b/corsika/detail/framework/geometry/Trajectory.inl
@@ -8,43 +8,149 @@
#pragma once
+#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
namespace corsika {
- template <typename TType>
- Point Trajectory<TType>::getPosition(double const u) const {
- return TType::getPosition(timeLength_ * u);
+ inline VelocityVector LineTrajectory::getVelocity(double const u) const {
+ return initialVelocity_ * (1 - u) + finalVelocity_ * u;
+ }
+
+ inline TimeType LineTrajectory::getDuration(double const u) const {
+ return u * timeStep_;
+ }
+
+ inline LengthType LineTrajectory::getLength(double const u) const {
+ if (timeLength_ == 0_s) return 0_m;
+ if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s)
+ return std::numeric_limits<LengthType::value_type>::infinity() * 1_m;
+ return getDistance(u) * timeStep_ / timeLength_;
+ }
+
+ ///! set new duration along potentially bend trajectory.
+ inline void LineTrajectory::setLength(LengthType const limit) {
+ setDuration(line_.getTimeFromArclength(limit));
+ }
+
+ ///! set new duration along potentially bend trajectory.
+ // Scale other properties by "limit/timeLength_"
+ inline void LineTrajectory::setDuration(TimeType const limit) {
+ if (timeStep_ == 0_s) {
+ timeLength_ = 0_s;
+ setFinalVelocity(getVelocity(0));
+ timeStep_ = limit;
+ } else {
+ // for infinite steps there can't be a difference between
+ // curved and straight trajectory, this is fundamentally
+ // undefined: assume they are the same (which, i.e. is always correct for a
+ // straight line trajectory).
+ //
+ // Final note: only straight-line trajectories should have
+ // infinite steps! Everything else is ill-defined.
+ if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s ||
+ timeLength_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s) {
+ timeLength_ = limit;
+ timeStep_ = limit;
+ // ...and don't touch velocity
+ } else {
+ const double scale = limit / timeStep_;
+ timeLength_ *= scale;
+ setFinalVelocity(getVelocity(scale));
+ timeStep_ = limit;
+ }
+ }
+ }
+
+ inline LengthType LineTrajectory::getDistance(double const u) const {
+ assert(u <= 1);
+ assert(u >= 0);
+ return line_.getArcLength(0 * second, u * timeLength_);
}
- template <typename TType>
- TimeType Trajectory<TType>::getDuration() const {
- return timeLength_;
+ inline Line const LeapFrogTrajectory::getLine() const {
+ auto D = getPosition(1) - getPosition(0);
+ auto d = D.getNorm();
+ 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).getNorm();
+ }
+ return Line(getPosition(0), v);
}
- template <typename TType>
- LengthType Trajectory<TType>::getLength() const {
- return getDistance(timeLength_);
+ inline Point LeapFrogTrajectory::getPosition(double const u) const {
+ Point position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
+ VelocityVector velocity =
+ initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+ return position + velocity * timeStep_ * u / 2;
}
- template <typename TType>
- LengthType Trajectory<TType>::getDistance(TimeType const t) const {
- assert(t <= timeLength_);
- assert(t >= 0 * second);
- return TType::getArcLength(0 * second, t);
+ inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
+ return initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
}
- template <typename TType>
- void Trajectory<TType>::getLimitEndTo(LengthType const limit) {
- timeLength_ = TType::getTimeFromArclength(limit);
+ inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
+ return getVelocity(u).normalized();
}
- template <typename TType>
- auto Trajectory<TType>::getNormalizedDirection() const {
- static_assert(std::is_same_v<TType, corsika::Line>);
- return TType::getVelocity().normalized();
+ ///! duration along potentially bend trajectory
+ inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
+ return u * timeStep_ *
+ (double(getVelocity(u).getNorm() / initialVelocity_.getNorm()) + 1.0) / 2;
}
+ ///! total length along potentially bend trajectory
+ inline LengthType LeapFrogTrajectory::getLength(double const u) const {
+ return timeStep_ * initialVelocity_.getNorm() * u;
+ }
+
+ ///! set new duration along potentially bend trajectory.
+ inline void LeapFrogTrajectory::setLength(LengthType const limit) {
+ if (initialVelocity_.getNorm() == 0_m / 1_s) setDuration(0_s);
+ setDuration(limit / initialVelocity_.getNorm());
+ }
+
+ ///! set new duration along potentially bend trajectory.
+ // Scale other properties by "limit/timeLength_"
+ inline void LeapFrogTrajectory::setDuration(TimeType const limit) { timeStep_ = limit; }
+
+ /*
+ template <typename TType>
+ Point Trajectory<TType>::getPosition(double const u) const {
+ return TType::getPosition(timeLength_ * u);
+ }
+
+ template <typename TType>
+ TimeType Trajectory<TType>::getDuration() const {
+ return timeLength_;
+ }
+
+ template <typename TType>
+ LengthType Trajectory<TType>::getLength() const {
+ return getDistance(timeLength_);
+ }
+
+ template <typename TType>
+ LengthType Trajectory<TType>::getDistance(TimeType const t) const {
+ assert(t <= timeLength_);
+ assert(t >= 0 * second);
+ return TType::getArcLength(0 * second, t);
+ }
+
+ template <typename TType>
+ void Trajectory<TType>::getLimitEndTo(LengthType const limit) {
+ timeLength_ = TType::getTimeFromArclength(limit);
+ }
+
+ template <typename TType>
+ auto Trajectory<TType>::getNormalizedDirection() const {
+ static_assert(std::is_same_v<TType, corsika::Line>);
+ return TType::getVelocity().normalized();
+ }
+ */
+
} // namespace corsika
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 4cfeb1ab5aa7fd189711954505dc00cc0cec74d5..71cc577582bb078a566603b604a900c3c8025584 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -159,7 +159,7 @@ namespace corsika {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
- if (lambda_inv_select < lambda_inv_sum) {
+ if (lambda_inv_select <= lambda_inv_sum) {
A_.doInteraction(view);
return ProcessReturn::Interacted;
}
@@ -174,7 +174,7 @@ namespace corsika {
lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
// soon as SecondaryView::parent() is migrated!
// check if we should execute THIS process and then EXIT
- if (lambda_inv_select < lambda_inv_sum) {
+ if (lambda_inv_select <= lambda_inv_sum) {
B_.doInteraction(view);
return ProcessReturn::Interacted;
}
@@ -218,7 +218,7 @@ namespace corsika {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A_.getInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
- if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select <
+ if (decay_inv_select <= decay_inv_sum) { // more pedagogical: rndm_select <
// decay_inv_sum / decay_inv_tot
A_.doDecay(view);
return ProcessReturn::Decayed;
@@ -232,7 +232,7 @@ namespace corsika {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B_.getInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
- if (decay_inv_select < decay_inv_sum) {
+ if (decay_inv_select <= decay_inv_sum) {
B_.doDecay(view);
return ProcessReturn::Decayed;
}
diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
new file mode 100644
index 0000000000000000000000000000000000000000..cc43ec7963de8db08bf40599d501fea2e2f44d88
--- /dev/null
+++ b/corsika/detail/framework/utility/QuarticSolver.inl
@@ -0,0 +1,147 @@
+/***************************************************************************
+ * Copyright (C) 2016 by Саша Миленковић *
+ * sasa.milenkovic.xyz@gmail.com *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * ( http://www.gnu.org/licenses/gpl-3.0.en.html ) *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#pragma once
+
+#include <cmath>
+
+namespace corsika::quartic_solver {
+
+ //---------------------------------------------------------------------------
+ // solve cubic equation x^3 + a*x^2 + b*x + c = 0
+ // x - array of size 3
+ // In case 3 real roots: => x[0], x[1], x[2], return 3
+ // 2 real roots: x[0], x[1], return 2
+ // 1 real root : x[0], x[1] ± i*x[2], return 1
+ unsigned int solveP3(double* x, double a, double b, double c) {
+ double a2 = a * a;
+ double q = (a2 - 3 * b) / 9;
+ double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
+ double r2 = r * r;
+ double q3 = q * q * q;
+ double A, B;
+ if (r2 < q3) {
+ double t = r / sqrt(q3);
+ if (t < -1) t = -1;
+ if (t > 1) t = 1;
+ t = acos(t);
+ a /= 3;
+ q = -2 * sqrt(q);
+ x[0] = q * cos(t / 3) - a;
+ x[1] = q * cos((t + M_2PI) / 3) - a;
+ x[2] = q * cos((t - M_2PI) / 3) - a;
+ return 3;
+ } else {
+ A = -pow(fabs(r) + sqrt(r2 - q3), 1. / 3);
+ if (r < 0) A = -A;
+ B = (0 == A ? 0 : q / A);
+
+ a /= 3;
+ x[0] = (A + B) - a;
+ x[1] = -0.5 * (A + B) - a;
+ x[2] = 0.5 * sqrt(3.) * (A - B);
+ if (fabs(x[2]) < eps) {
+ x[2] = x[1];
+ return 2;
+ }
+
+ return 1;
+ }
+ }
+
+ //---------------------------------------------------------------------------
+ // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
+ // Attention - this function returns dynamically allocated array. It has to be released
+ // afterwards.
+ DComplex* solve_quartic(double a, double b, double c, double d) {
+ double a3 = -b;
+ double b3 = a * c - 4. * d;
+ double c3 = -a * a * d - c * c + 4. * b * d;
+
+ // cubic resolvent
+ // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+
+ double x3[3];
+ unsigned int iZeroes = solveP3(x3, a3, b3, c3);
+
+ double q1, q2, p1, p2, D, sqD, y;
+
+ y = x3[0];
+ // The essence - choosing Y with maximal absolute value.
+ if (iZeroes != 1) {
+ if (fabs(x3[1]) > fabs(y)) y = x3[1];
+ if (fabs(x3[2]) > fabs(y)) y = x3[2];
+ }
+
+ // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q)
+
+ D = y * y - 4 * d;
+ if (fabs(D) < eps) // in other words - D==0
+ {
+ q1 = q2 = y * 0.5;
+ // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g)
+ D = a * a - 4 * (b - y);
+ if (fabs(D) < eps) // in other words - D==0
+ p1 = p2 = a * 0.5;
+
+ else {
+ sqD = sqrt(D);
+ p1 = (a + sqD) * 0.5;
+ p2 = (a - sqD) * 0.5;
+ }
+ } else {
+ sqD = sqrt(D);
+ q1 = (y + sqD) * 0.5;
+ q2 = (y - sqD) * 0.5;
+ // g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer
+ p1 = (a * q1 - c) / (q1 - q2);
+ p2 = (c - a * q2) / (q1 - q2);
+ }
+
+ DComplex* retval = new DComplex[4];
+
+ // solving quadratic eq. - x^2 + p1*x + q1 = 0
+ D = p1 * p1 - 4 * q1;
+ if (D < 0.0) {
+ retval[0].real(-p1 * 0.5);
+ retval[0].imag(sqrt(-D) * 0.5);
+ retval[1] = std::conj(retval[0]);
+ } else {
+ sqD = sqrt(D);
+ retval[0].real((-p1 + sqD) * 0.5);
+ retval[1].real((-p1 - sqD) * 0.5);
+ }
+
+ // solving quadratic eq. - x^2 + p2*x + q2 = 0
+ D = p2 * p2 - 4 * q2;
+ if (D < 0.0) {
+ retval[2].real(-p2 * 0.5);
+ retval[2].imag(sqrt(-D) * 0.5);
+ retval[3] = std::conj(retval[2]);
+ } else {
+ sqD = sqrt(D);
+ retval[2].real((-p2 + sqD) * 0.5);
+ retval[3].real((-p2 - sqD) * 0.5);
+ }
+
+ return retval;
+ }
+
+} // namespace corsika::quartic_solver
\ No newline at end of file
diff --git a/corsika/detail/media/BaseExponential.inl b/corsika/detail/media/BaseExponential.inl
index b8482d02e9fbedd0627efd3029957171bdaffcb1..2ec244091fece18c6149fbc4cc4d9b22a4447b39 100644
--- a/corsika/detail/media/BaseExponential.inl
+++ b/corsika/detail/media/BaseExponential.inl
@@ -12,7 +12,6 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
namespace corsika {
@@ -23,8 +22,8 @@ namespace corsika {
template <typename TDerived>
GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
- Trajectory<Line> const& line, LengthType vL,
- Vector<dimensionless_d> const& axis) const {
+ setup::Trajectory const& line, LengthType vL,
+ DirectionVector const& axis) const {
if (vL == LengthType::zero()) { return GrammageType::zero(); }
auto const uDotA = line.getNormalizedDirection().dot(axis).magnitude();
@@ -39,8 +38,8 @@ namespace corsika {
template <typename TDerived>
LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
- Trajectory<Line> const& line, GrammageType grammage,
- Vector<dimensionless_d> const& axis) const {
+ setup::Trajectory const& line, GrammageType grammage,
+ DirectionVector const& axis) const {
auto const uDotA = line.getNormalizedDirection().dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(line.getStartPoint());
diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl
index e3b50839a11b3b07c04b67d5f6d7f2ec8bfe4655..63caea5556825ca09eba384d19528cd8c7be5d6b 100644
--- a/corsika/detail/media/FlatExponential.inl
+++ b/corsika/detail/media/FlatExponential.inl
@@ -11,7 +11,6 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/BaseExponential.hpp>
#include <corsika/media/NuclearComposition.hpp>
@@ -39,13 +38,13 @@ namespace corsika {
}
template <typename T>
- GrammageType FlatExponential<T>::getIntegratedGrammage(Trajectory<Line> const& line,
+ GrammageType FlatExponential<T>::getIntegratedGrammage(setup::Trajectory const& line,
LengthType to) const {
return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, to, axis_);
}
template <typename T>
- LengthType FlatExponential<T>::getArclengthFromGrammage(Trajectory<Line> const& line,
+ LengthType FlatExponential<T>::getArclengthFromGrammage(setup::Trajectory const& line,
GrammageType grammage) const {
return BaseExponential<FlatExponential<T>>::getArclengthFromGrammage(line, grammage,
axis_);
diff --git a/corsika/detail/media/HomogeneousMedium.inl b/corsika/detail/media/HomogeneousMedium.inl
index 1952d5f68f04db5877a2d807e5cc5149ad2f16ab..d73354aa680c085241c1f71c37daed4bb1221b47 100644
--- a/corsika/detail/media/HomogeneousMedium.inl
+++ b/corsika/detail/media/HomogeneousMedium.inl
@@ -11,7 +11,6 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/NuclearComposition.hpp>
namespace corsika {
@@ -32,13 +31,13 @@ namespace corsika {
}
template <typename T>
- GrammageType HomogeneousMedium<T>::getIntegratedGrammage(Trajectory<Line> const&,
+ GrammageType HomogeneousMedium<T>::getIntegratedGrammage(setup::Trajectory const&,
LengthType to) const {
return to * density_;
}
template <typename T>
- LengthType HomogeneousMedium<T>::getArclengthFromGrammage(Trajectory<Line> const&,
+ LengthType HomogeneousMedium<T>::getArclengthFromGrammage(setup::Trajectory const&,
GrammageType grammage) const {
return grammage / density_;
}
diff --git a/corsika/detail/media/InhomogeneousMedium.inl b/corsika/detail/media/InhomogeneousMedium.inl
index 33b9066c2a040bf74fa342ace617879cd15dd347..064e9610efe4ed04e7707cf626cae8f0a355c076 100644
--- a/corsika/detail/media/InhomogeneousMedium.inl
+++ b/corsika/detail/media/InhomogeneousMedium.inl
@@ -11,7 +11,6 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/NuclearComposition.hpp>
namespace corsika {
@@ -37,13 +36,13 @@ namespace corsika {
template <typename T, typename TDensityFunction>
GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
- Trajectory<Line> const& line, LengthType to) const {
+ setup::Trajectory const& line, LengthType to) const {
return densityFunction_.getIntegrateGrammage(line, to);
}
template <typename T, typename TDensityFunction>
LengthType InhomogeneousMedium<T, TDensityFunction>::getArclengthFromGrammage(
- Trajectory<Line> const& line, GrammageType grammage) const {
+ setup::Trajectory const& line, GrammageType grammage) const {
return densityFunction_.getArclengthFromGrammage(line, grammage);
}
diff --git a/corsika/detail/media/LinearApproximationIntegrator.inl b/corsika/detail/media/LinearApproximationIntegrator.inl
index 7d526ba03f28015fb9c9186478274a97bf7933e2..7452ff31d3047b235aaf292fedb78dd0b0552ba5 100644
--- a/corsika/detail/media/LinearApproximationIntegrator.inl
+++ b/corsika/detail/media/LinearApproximationIntegrator.inl
@@ -19,28 +19,28 @@ namespace corsika {
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
- Trajectory<Line> const& line, LengthType length) const {
+ setup::Trajectory const& line, LengthType length) const {
auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
- auto const c1 = getImplementation().rho_.getFirstDerivative(
- line.getPosition(0), line.getNormalizedDirection());
+ auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
+ line.getDirection(0));
return (c0 + 0.5 * c1 * length) * length;
}
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
- Trajectory<Line> const& line, GrammageType grammage) const {
+ setup::Trajectory const& line, GrammageType grammage) const {
auto const c0 = getImplementation().rho_(line.getPosition(0));
- auto const c1 = getImplementation().rho_.getFirstDerivative(
- line.getPosition(0), line.getNormalizedDirection());
+ auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
+ line.getDirection(0));
return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
}
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::getMaximumLength(
- Trajectory<Line> const& line, [[maybe_unused]] double relError) const {
+ setup::Trajectory const& line, [[maybe_unused]] double relError) const {
[[maybe_unused]] auto const c1 = getImplementation().rho_.getSecondDerivative(
- line.getPosition(0), line.getNormalizedDirection());
+ line.getPosition(0), line.getDirection(0));
// todo: provide a real, working implementation
return 1_m * std::numeric_limits<double>::infinity();
diff --git a/corsika/detail/media/ShowerAxis.inl b/corsika/detail/media/ShowerAxis.inl
index efb7dc58d304ddeec3656c4d9529831739aec3fe..1356197a881c0cd152aaa25559da667f4a66abc0 100644
--- a/corsika/detail/media/ShowerAxis.inl
+++ b/corsika/detail/media/ShowerAxis.inl
@@ -57,32 +57,42 @@ namespace corsika {
}
GrammageType ShowerAxis::getX(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 fraction = fractionalBin - lower;
+ unsigned int const upper = lower + 1;
- if (lower < 0) {
+ if (fractionalBin < 0) {
CORSIKA_LOG_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()) {
- const std::string err =
- fmt::format("shower axis too short, cannot extrapolate (l / max_length_ = {} )",
- l / max_length_);
- CORSIKA_LOG_ERROR(err);
- throw std::runtime_error(err.c_str());
+ CORSIKA_LOG_ERROR(
+ "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
+ l / max_length_);
+ if (throw_) {
+ const std::string err = fmt::format(
+ "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
+ l / max_length_);
+ throw std::runtime_error(err.c_str());
+ }
}
+ CORSIKA_LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}", fraction,
+ fractionalBin, lower, upper);
- assert(0 <= lambda && lambda <= 1.);
+ assert(0 <= fraction && fraction <= 1.);
- CORSIKA_LOG_TRACE("ShowerAxis::getX l={} m, lower={}, lambda={}, upper={}", l / 1_m,
- lower, lambda, upper);
+ CORSIKA_LOG_TRACE("ShowerAxis::getX l={} m, lower={}, fraction={}, upper={}", l / 1_m,
+ lower, fraction, upper);
// linear interpolation between getX[lower] and X[upper]
- return X_[upper] * lambda + X_[lower] * (1 - lambda);
+ return X_[upper] * fraction + X_[lower] * (1 - fraction);
}
LengthType ShowerAxis::getSteplength() const { return steplength_; }
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 5681239810d139e5be7d45c02aed0fb42beb18fd..cef313d25b2f0e9c4383b712f42c9fc3351409b8 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -12,14 +12,17 @@
namespace corsika {
- ObservationPlane::ObservationPlane(Plane const& obsPlane, std::string const& filename,
- bool deleteOnHit)
+ ObservationPlane::ObservationPlane(Plane const& obsPlane, DirectionVector const& x_axis,
+ std::string const& filename, bool deleteOnHit)
: plane_(obsPlane)
, outputStream_(filename)
, deleteOnHit_(deleteOnHit)
, energy_ground_(0_GeV)
- , count_ground_(0) {
- outputStream_ << "#PDG code, energy / eV, distance to center / m" << std::endl;
+ , count_ground_(0)
+ , xAxis_(x_axis.normalized())
+ , yAxis_(obsPlane.getNormal().cross(xAxis_)) {
+ outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m"
+ << std::endl;
}
ProcessReturn ObservationPlane::doContinuous(
@@ -37,8 +40,11 @@ namespace corsika {
}
const auto energy = particle.getEnergy();
+ auto const displacement = trajectory.getPosition(1) - plane_.getCenter();
+
outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
- << ' '
+ << ' ' << displacement.dot(xAxis_) / 1_m << ' '
+ << displacement.dot(yAxis_) / 1_m
<< (trajectory.getPosition(1) - plane_.getCenter()).getNorm() / 1_m
<< std::endl;
@@ -53,9 +59,71 @@ namespace corsika {
}
LengthType ObservationPlane::getMaxStepLength(
- corsika::setup::Stack::particle_type const&,
+ corsika::setup::Stack::particle_type const& vParticle,
corsika::setup::Trajectory const& trajectory) {
+ int chargeNumber;
+ if (is_nucleus(vParticle.getPID())) {
+ chargeNumber = vParticle.getNuclearZ();
+ } else {
+ chargeNumber = get_charge_number(vParticle.getPID());
+ }
+ auto const* currentLogicalVolumeNode = vParticle.getNode();
+ auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
+ vParticle.getPosition());
+ auto direction = trajectory.getLine().getV0().normalized();
+
+ if (chargeNumber != 0 &&
+ abs(plane_.getNormal().dot(trajectory.getLine().getV0().cross(magneticfield))) *
+ 1_s / 1_m / 1_T >
+ 1e-6) {
+ auto const* currentLogicalVolumeNode = vParticle.getNode();
+ auto magneticfield =
+ currentLogicalVolumeNode->getModelProperties().getMagneticField(
+ vParticle.getPosition());
+ auto k =
+ chargeNumber * constants::c * 1_eV / (vParticle.getMomentum().norm() * 1_V);
+
+ if (direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
+ (plane_.getNormal().dot(trajectory.getLine().getR0() - plane_.getCenter()) *
+ plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
+ 0) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+
+ LengthType MaxStepLength1 =
+ (sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
+ (plane_.getNormal().dot(trajectory.getLine().getR0() -
+ plane_.getCenter()) *
+ plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane_.getNormal()) / direction.getNorm()) /
+ (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
+
+ LengthType MaxStepLength2 =
+ (-sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
+ (plane_.getNormal().dot(trajectory.getLine().getR0() -
+ plane_.getCenter()) *
+ plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
+ direction.dot(plane_.getNormal()) / direction.getNorm()) /
+ (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
+
+ if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
+ std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
+ return MaxStepLength2 *
+ (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
+ .norm() *
+ 1.001;
+ } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
+ std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
+ return MaxStepLength1 *
+ (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
+ .norm() *
+ 1.001;
+ }
+ }
+
TimeType const timeOfIntersection =
(plane_.getCenter() - trajectory.getStartPoint()).dot(plane_.getNormal()) /
trajectory.getVelocity().dot(plane_.getNormal());
@@ -66,7 +134,8 @@ namespace corsika {
auto const pointOfIntersection = trajectory.getPosition(timeOfIntersection);
auto dist = (trajectory.getStartPoint() - pointOfIntersection).getNorm() * 1.0001;
- CORSIKA_LOG_TRACE("ObservationPlane::getMaxStepLength l={} m", dist / 1_m);
+ CORSIKA_LOG_TRACE("ObservationPlane w/o magnetic field: getMaxStepLength l={} m",
+ dist / 1_m);
return dist;
}
diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl
index bee651660b57d7c2e739bc9520cbc582d5ce6d9e..2812a0bd3734734f074820c60cad6c80bdb9e7ed 100644
--- a/corsika/detail/modules/TrackWriter.inl
+++ b/corsika/detail/modules/TrackWriter.inl
@@ -25,8 +25,9 @@ namespace corsika {
using namespace std::string_literals;
file_.open(filename_);
- file_ << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
- << '\n';
+ file_
+ << "# PID, E / eV, start coordinates / m, displacement vector to end / m, steplength / m "s
+ << '\n';
}
template <typename TParticle, typename TTrack>
@@ -43,7 +44,9 @@ namespace corsika {
<< std::setw(width_) << std::scientific << std::setprecision(precision_) << start[2] / 1_m
<< std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[0] / 1_m
<< std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[1] / 1_m
- << std::setw(width_) << std::scientific << std::setprecision(precision_) << delta[2] / 1_m << '\n';
+ << std::setw(width) << std::scientific << std::setprecision(precision) << delta[2] / 1_m
+ << std::setw(width) << std::scientific << std::setprecision(precision) << delta.norm() / 1_m
+ << '\n';
// clang-format on
return ProcessReturn::Ok;
diff --git a/corsika/detail/stack/history/HistoryObservationPlane.hpp b/corsika/detail/stack/history/HistoryObservationPlane.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0f6170393c8dc64261d514d08f96d99cbc39e912
--- /dev/null
+++ b/corsika/detail/stack/history/HistoryObservationPlane.hpp
@@ -0,0 +1,21 @@
+#pragma once
+
+#include <boost/histogram.hpp>
+
+namespace corsika::history {
+
+ namespace detail {
+ inline auto hist_factory() {
+ namespace bh = boost::histogram;
+ namespace bha = bh::axis;
+ auto h = bh::make_histogram(
+ bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11, "muon energy"},
+ bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11,
+ "projectile energy"},
+ bha::category<int, bh::use_default, bha::option::growth_t>{
+ {211, -211, 2212, -2212}, "projectile PDG"});
+ return h;
+ }
+ } // namespace detail
+
+} // namespace corsika::history
diff --git a/corsika/detail/stack/history/HistoryObservationPlane.inl b/corsika/detail/stack/history/HistoryObservationPlane.inl
new file mode 100644
index 0000000000000000000000000000000000000000..0536eaf3ada84db7382461ec95d2ff7cc6ffac1b
--- /dev/null
+++ b/corsika/detail/stack/history/HistoryObservationPlane.inl
@@ -0,0 +1,80 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/logging/Logging.h>
+#include <corsika/stack/history/HistoryObservationPlane.hpp>
+
+#include <boost/histogram/ostream.hpp>
+
+namespace corsika::history {
+
+ HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
+ Plane const& obsPlane,
+ bool deleteOnHit)
+ : stack_{stack}
+ , plane_{obsPlane}
+ , deleteOnHit_{deleteOnHit} {}
+
+ ProcessReturn HistoryObservationPlane::DoContinuous(
+ setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.getCenter() - trajectory.getR0()).dot(plane_.getNormal()) /
+ trajectory.getV0().dot(plane_.getNormal());
+
+ if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
+
+ if (plane_.isAbove(trajectory.getR0()) == plane_.isAbove(trajectory.getPosition(1))) {
+ return ProcessReturn::Ok;
+ }
+
+ CORSIKA_LOG_DEBUG(fmt::format("HistoryObservationPlane: Particle detected: pid={}",
+ particle.getPID()));
+
+ auto const pid = particle.getPID();
+ if (particles::isMuon(pid)) { fillHistoryHistogram(particle); }
+
+ if (deleteOnHit_) {
+ return ProcessReturn::ParticleAbsorbed;
+ } else {
+ return ProcessReturn::Ok;
+ }
+ }
+
+ LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
+ setup::Trajectory const& trajectory) {
+ TimeType const timeOfIntersection =
+ (plane_.getCenter() - trajectory.getR0()).dot(plane_.getNormal()) /
+ trajectory.getV0().dot(plane_.getNormal());
+
+ if (timeOfIntersection < TimeType::zero()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+
+ auto const pointOfIntersection = trajectory.getPosition(timeOfIntersection);
+ return (trajectory.getR0() - pointOfIntersection).norm() * 1.0001;
+ }
+
+ void HistoryObservationPlane::fillHistoryHistogram(
+ setup::Stack::ParticleType const& muon) {
+ double const muon_energy = muon.getEnergy() / 1_GeV;
+
+ int genctr{0};
+ Event const* event = muon.getEvent().get();
+ while (event) {
+ auto const projectile = stack_.cfirst() + event->projectileIndex();
+ if (event->eventType() == EventType::Interaction) {
+ genctr++;
+ double const projEnergy = projectile.getEnergy() / 1_GeV;
+ int const pdg = static_cast<int>(particles::getPDG(projectile.getPID()));
+
+ histogram_(muon_energy, projEnergy, pdg);
+ }
+ event = event->parentEvent().get(); // projectile.getEvent().get();
+ }
+ }
+} // namespace corsika::history
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index ba0040c887291baf5706a2772e089478f2ab15a6..fcd8eaf2778618ad4dc6e470448c4e1f36e795c7 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -18,8 +18,13 @@
#include <corsika/media/Environment.hpp>
#include <corsika/framework/logging/Logging.hpp>
+/* see Issue 161, we need to include SetupStack only because we need
+ to globally define StackView. This is clearly not nice and should
+ be changed, when possible. It might be that StackView needs to be
+ templated in Cascade, but this would be even worse... so we don't
+ do that until it is really needed.
+ */
#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
#include <cassert>
#include <cmath>
@@ -43,7 +48,7 @@ namespace corsika {
*
* <code>
* auto getTrack(Particle const& p)</auto>,
- * with the return type <code>geometry::Trajectory<corsika::Line>
+ * with the return type <code>geometry::Trajectory<Line>
* </code>
*
* <b>TProcessList</b> must be a ProcessSequence. *
@@ -53,6 +58,11 @@ namespace corsika {
*
*/
template <typename TTracking, typename TProcessList, typename TStack,
+ /*
+ TStackView is needed as explicit template parameter because
+ of issue 161 and the
+ inability of clang to understand "stack::MakeView" so far.
+ */
typename TStackView = corsika::setup::StackView>
class Cascade {
@@ -62,10 +72,19 @@ namespace corsika {
typedef typename VolumeTreeNode::IModelProperties MediumInterface;
public:
+ /**
+ * \group constructors
+ * \{
+ * Cascade class cannot be default constructed, but needs a valid
+ * list of physics processes for configuration at construct time.
+ */
Cascade() = delete;
-
- Cascade(corsika::Environment<MediumInterface> const& env, TTracking& tr,
- TProcessList& pl, TStack& stack)
+ Cascade(Cascade const&) = default;
+ Cascade(Cascade&&) = default;
+ ~Cascade() = default;
+ Cascade & operator=(Cascade const&)) = default;
+ Cascade(Environment<MediumInterface> const& env, TTracking& tr, TProcessList& pl,
+ TStack& stack)
: environment_(env)
, tracking_(tr)
, sequence_(pl)
@@ -75,12 +94,7 @@ namespace corsika {
CORSIKA_LOG_INFO(" - With full cascade HISTORY.");
}
}
-
- /**
- * The Init function is called before the actual cascade simulations.
- * All components of the Cascade simulation must be configured here.
- */
- void init();
+ //! \}
/**
* set the nodes for all particles on the stack according to their numerical
@@ -102,6 +116,16 @@ namespace corsika {
void forceInteraction();
private:
+ /**
+ * The Step function is executed for each particle from the
+ * stack. It will calcualte geometric transport of the particles,
+ * and apply continuous and stochastic processes to it, which may
+ * lead to energy losses, scattering, absorption, decays and the
+ * production of secondary particles.
+ *
+ * New particles produced in one step are subject to further
+ * processing, e.g. thinning, etc.
+ */
void step(Particle& vParticle);
ProcessReturn decay(TStackView& view);
@@ -109,12 +133,11 @@ namespace corsika {
void setEventType(TStackView& view, history::EventType);
// data members
- corsika::Environment<MediumInterface> const& environment_;
+ Environment<MediumInterface> const& environment_;
TTracking& tracking_;
TProcessList& sequence_;
TStack& stack_;
- corsika::default_prng_type& rng_ =
- corsika::RNGManager::getInstance().getRandomStream("cascade");
+ default_prng_type& rng_ = RNGManager::getInstance().getRandomStream("cascade");
unsigned int count_ = 0;
// but this here temporarily. Should go into dedicated file later:
diff --git a/corsika/framework/core/PhysicalUnits.hpp b/corsika/framework/core/PhysicalUnits.hpp
index 71d91c9a679518099d432e60e24340f5bddd7ef3..a9769c00d3fcae0027c0602a78d60255e8169f72 100644
--- a/corsika/framework/core/PhysicalUnits.hpp
+++ b/corsika/framework/core/PhysicalUnits.hpp
@@ -94,6 +94,8 @@ namespace corsika::units::si {
phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
using InverseGrammageType =
phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
+ using MagneticFluxType =
+ phys::units::quantity<phys::units::magnetic_flux_density_d, double>;
template <typename DimFrom, typename DimTo>
auto constexpr conversion_factor_HEP_to_SI() {
diff --git a/corsika/framework/geometry/Helix.hpp b/corsika/framework/geometry/Helix.hpp
index 826b4c3fc055d992fccc13f650d8d27ae307ef6b..c69d65dbb57b391218c599bef3aa67586394b475 100644
--- a/corsika/framework/geometry/Helix.hpp
+++ b/corsika/framework/geometry/Helix.hpp
@@ -48,6 +48,8 @@ namespace corsika {
inline Point getPosition(TimeType const t) const;
+ VelocityVec getVelocity(TimeType const t) const;
+
inline Point getPositionFromArclength(LengthType const l) const;
inline LengthType getArcLength(TimeType const t1, TimeType const t2) const;
diff --git a/corsika/framework/geometry/Intersections.hpp b/corsika/framework/geometry/Intersections.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4db60b859d91a9050a6438dd6bd5b8ed8d63cd8d
--- /dev/null
+++ b/corsika/framework/geometry/Intersections.hpp
@@ -0,0 +1,56 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <map> // for pair
+
+namespace corsika {
+
+ /**
+ *
+ * Container to store and return a list of intersections of a
+ * trajectory with a geometric volume objects in space.
+ *
+ **/
+
+ class Intersections {
+
+ Intersections(const Intersections&) = delete;
+ Intersections(Intersections&&) = delete;
+ Intersections& operator=(const Intersections&) = delete;
+
+ public:
+ Intersections()
+ : has_intersections_(false) {}
+
+ Intersections(TimeType&& t1, TimeType&& t2)
+ : has_intersections_(true)
+ , intersections_(std::make_pair(t1, t2)) {}
+
+ Intersections(TimeType&& t)
+ : has_intersections_(true)
+ , intersections_(
+ std::make_pair(t, std::numeric_limits<TimeType::value_type>::infinity() * second)) {}
+
+ bool hasIntersections() const { return has_intersections_; }
+
+ ///! where did the trajectory currently enter the volume
+ TimeType getEntry() const { return intersections_.first; }
+
+ ///! where did the trajectory currently exit the volume
+ TimeType getExit() const { return intersections_.second; }
+
+ private:
+ bool has_intersections_;
+ std::pair<TimeType, TimeType> intersections_;
+ };
+
+} // namespace corsika
diff --git a/corsika/framework/geometry/Line.hpp b/corsika/framework/geometry/Line.hpp
index ea99139db79179c188dc856171ccb0f421928955..ba7381a4a9398f219dd30be910405cab887841b6 100644
--- a/corsika/framework/geometry/Line.hpp
+++ b/corsika/framework/geometry/Line.hpp
@@ -11,26 +11,30 @@ n/*
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
namespace corsika {
/**
- * Describes a straight line in space
*
- */
+ * A Line describes a movement in three dimensional space. It
+ * consists of a Point `$\vec{p_0}$` and and a speed-Vector
+ * `$\vec{v}$`, so that it can return GetPosition as
+ * `$\vec{p_0}*\vec{v}*t$` for any value of time `$t$`.
+ *
+ **/
class Line {
- ///! \toto move this to PhysicalUnits
- using VelocityVec = Vector<SpeedType::dimension_type>;
-
public:
- Line(Point const& pR0, VelocityVec const& pV0)
+ Line(Point const& pR0, VelocityVector const& pV0)
: start_point_(pR0)
, velocity_(pV0) {}
inline Point getPosition(TimeType const t) const;
+ inline VelocityVector const& getVelocity(TimeType const) const;
+
inline Point getPositionFromArclength(LengthType const l) const;
inline LengthType getArcLength(TimeType const t1, TimeType const t2) const;
@@ -38,14 +42,12 @@ namespace corsika {
inline TimeType getTimeFromArclength(LengthType const t) const;
inline Point const& getStartPoint() const;
- inline Point& startPoint() { return start_point_; }
- inline VelocityVec const& getVelocity() const;
- inline VelocityVec& velocity() { return velocity_; }
+ inline VelocityVector const& getVelocity() const;
private:
Point start_point_;
- VelocityVec velocity_;
+ VelocityVector velocity_;
};
} // namespace corsika
diff --git a/corsika/framework/geometry/Sphere.hpp b/corsika/framework/geometry/Sphere.hpp
index 0cf5e66d781a12b8df5e04d6d4872fda97ecdb5a..3b19d7f9b7635159dd29c8457f5816def584e8c5 100644
--- a/corsika/framework/geometry/Sphere.hpp
+++ b/corsika/framework/geometry/Sphere.hpp
@@ -10,6 +10,7 @@ n/*
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/IVolume.hpp>
namespace corsika {
@@ -27,14 +28,16 @@ namespace corsika {
, radius_(pRadius) {}
//! returns true if the Point p is within the sphere
- inline bool isInside(Point const& p) const override;
+ bool isInside(Point const& p) const override;
- inline Point const& getCenter() const;
- inline Point& getCenter() { return center_; }
+ Point const& getCenter() const;
- inline LengthType getRadius() const;
- inline LengthType& getRadius() { return radius_; }
+ void setCenter(Point const&);
+ LengthType getRadius() const;
+
+ void setRadius(LengthType const);
+
private:
Point center_;
LengthType radius_;
diff --git a/corsika/framework/geometry/Trajectory.hpp b/corsika/framework/geometry/Trajectory.hpp
index 9f5ef956f6bfe3ad9903a58b4574581981335b02..6a3700eea9bf1879520eee2ad6d4a18d22fc81f5 100644
--- a/corsika/framework/geometry/Trajectory.hpp
+++ b/corsika/framework/geometry/Trajectory.hpp
@@ -11,13 +11,182 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
namespace corsika {
+ /**
+ *
+ * A Trajectory is a description of a momvement of an object in
+ * three-dimensional space that describes the trajectory (connection
+ * between two Points in space), as well as the direction of motion
+ * at any given point.
+ *
+ * A Trajectory has a start `0` and an end `1`, where
+ * e.g. getPosition(0) returns the start point and getDirection(1)
+ * the direction of motion at the end. Values outside 0...1 are not
+ * defined.
+ *
+ * A Trajectory has a length in [m], getLength, a duration in [s], getDuration.
+ *
+ * Note: so far it is assumed that the speed (d|vec{r}|/dt) between
+ * start and end does not change and is constant for the entire
+ * Trajectory.
+ *
+ **/
+
+ class LineTrajectory {
+
+ public:
+ LineTrajectory() = delete;
+ LineTrajectory(LineTrajectory const&) = default;
+ LineTrajectory(LineTrajectory&&) = default;
+ LineTrajectory& operator=(LineTrajectory const&) = delete;
+
+ /**
+ * \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
+ */
+ LineTrajectory(Line const& theLine, TimeType timeLength)
+ : line_(theLine)
+ , timeLength_(timeLength)
+ , timeStep_(timeLength)
+ , initialVelocity_(theLine.getVelocity(TimeType::zero()))
+ , finalVelocity_(theLine.getVelocity(timeLength)) {}
+
+ /**
+ * \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
+ * start of trajectory \param finalV Final velocity vector at start of trajectory
+ */
+ LineTrajectory(Line const& theLine,
+ TimeType const timeLength, // length of theLine (straight)
+ TimeType const timeStep, // length of bend step (curved)
+ VelocityVector const& initialV, VelocityVector const& finalV)
+ : line_(theLine)
+ , timeLength_(timeLength)
+ , timeStep_(timeStep)
+ , initialVelocity_(initialV)
+ , finalVelocity_(finalV) {}
+
+ Line const& getLine() const { return line_; }
+
+ Point getPosition(double const u) const { return line_.getPosition(timeLength_ * u); }
+
+ VelocityVector getVelocity(double const u) const;
+
+ DirectionVector getDirection(double const u) const {
+ return getVelocity(u).normalized();
+ }
+
+ ///! duration along potentially bend trajectory
+ TimeType getDuration(double const u = 1) const;
+
+ ///! total length along potentially bend trajectory
+ LengthType getLength(double const u = 1) const;
+
+ ///! set new duration along potentially bend trajectory.
+ void setLength(LengthType const limit);
+
+ ///! set new duration along potentially bend trajectory.
+ // Scale other properties by "limit/timeLength_"
+ void setDuration(TimeType const limit);
+
+ protected:
+ ///! total length along straight trajectory
+ LengthType getDistance(double const u) const;
+
+ void setFinalVelocity(VelocityVector const& v) { finalVelocity_ = v; }
+
+ private:
+ Line line_;
+ TimeType timeLength_; ///! length of straight step (shortest connecting line)
+ TimeType timeStep_; ///! length of bend step (curved)
+ VelocityVector initialVelocity_;
+ VelocityVector finalVelocity_;
+ };
+
+ /**
+ * The LeapFrogTrajectory stores information on one leap-frog step.
+ *
+ * The leap-frog algorithm uses a half-step and is used in magnetic
+ * field tracking. The LeapFrogTrajectory will solve the leap-frog
+ * algorithm equation for a given constant $k$ that has to be
+ * specified during construction (essentially fixing the magnetic
+ * field). Thus, different steps (length) can be dynamically
+ * generated here. The velocity vector will correctly point into the
+ * direction as calculated by the algorithm for any steplength, or
+ * intermediate position.
+ *
+ **/
+
+ class LeapFrogTrajectory {
+
+ public:
+ LeapFrogTrajectory() = delete;
+ LeapFrogTrajectory(LeapFrogTrajectory const&) = default;
+ LeapFrogTrajectory(LeapFrogTrajectory&&) = default;
+ LeapFrogTrajectory& operator=(LeapFrogTrajectory const&) = delete;
+
+ LeapFrogTrajectory(Point const& pos, VelocityVector const& initialVelocity,
+ MagneticFieldVector const& Bfield,
+ decltype(square(meter) / (square(second) * volt)) const k,
+ TimeType const timeStep) // leap-from total length
+ : initialPosition_(pos)
+ , initialVelocity_(initialVelocity)
+ , initialDirection_(initialVelocity.normalized())
+ , magneticfield_(Bfield)
+ , k_(k)
+ , timeStep_(timeStep) {}
+
+ Line const getLine() const;
+
+ Point getPosition(double const u) const;
+
+ VelocityVector getVelocity(double const u) const;
+
+ DirectionVector getDirection(double const u) const;
+
+ ///! duration along potentially bend trajectory
+ TimeType getDuration(double const u = 1) const;
+
+ ///! total length along potentially bend trajectory
+ LengthType getLength(double const u = 1) const;
+
+ ///! set new duration along potentially bend trajectory.
+ void setLength(LengthType const limit);
+
+ ///! set new duration along potentially bend trajectory.
+ // Scale other properties by "limit/timeLength_"
+ void setDuration(TimeType const limit);
+
+ private:
+ Point initialPosition_;
+ VelocityVector initialVelocity_;
+ DirectionVector initialDirection_;
+ MagneticFieldVector magneticfield_;
+ decltype(square(meter) / (square(second) * volt)) k_;
+ TimeType timeStep_;
+ };
+
+ /*
+
template <typename TType>
class Trajectory : public TType {
public:
+
using TType::getArcLength;
using TType::getPosition;
@@ -40,6 +209,7 @@ namespace corsika {
private:
TimeType timeLength_;
};
+ */
} // namespace corsika
diff --git a/corsika/framework/utility/QuarticSolver.hpp b/corsika/framework/utility/QuarticSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..912c19d8d56012c42e4be9faf65969d26c4cc9b3
--- /dev/null
+++ b/corsika/framework/utility/QuarticSolver.hpp
@@ -0,0 +1,73 @@
+/***************************************************************************
+ * Copyright (C) 2016 by Саша Миленковић *
+ * sasa.milenkovic.xyz@gmail.com *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * ( http://www.gnu.org/licenses/gpl-3.0.en.html ) *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#pragma once
+
+#include <complex>
+
+/**
+ * \todo convert to class
+ */
+
+namespace corsika::quartic_solver {
+
+ const double PI = 3.141592653589793238463L;
+ const double M_2PI = 2 * PI;
+ const double eps = 1e-12;
+
+ typedef std::complex<double> DComplex;
+
+ //---------------------------------------------------------------------------
+ // useful for testing
+ inline DComplex polinom_2(DComplex x, double a, double b) {
+ // Horner's scheme for x*x + a*x + b
+ return x * (x + a) + b;
+ }
+
+ //---------------------------------------------------------------------------
+ // useful for testing
+ inline DComplex polinom_3(DComplex x, double a, double b, double c) {
+ // Horner's scheme for x*x*x + a*x*x + b*x + c;
+ return x * (x * (x + a) + b) + c;
+ }
+
+ //---------------------------------------------------------------------------
+ // useful for testing
+ inline DComplex polinom_4(DComplex x, double a, double b, double c, double d) {
+ // Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
+ return x * (x * (x * (x + a) + b) + c) + d;
+ }
+
+ //---------------------------------------------------------------------------
+ // x - array of size 3
+ // In case 3 real roots: => x[0], x[1], x[2], return 3
+ // 2 real roots: x[0], x[1], return 2
+ // 1 real root : x[0], x[1] ± i*x[2], return 1
+ unsigned int solveP3(double* x, double a, double b, double c);
+
+ //---------------------------------------------------------------------------
+ // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
+ // Attention - this function returns dynamically allocated array. It has to be released
+ // afterwards.
+ DComplex* solve_quartic(double a, double b, double c, double d);
+
+} // namespace corsika::quartic_solver
+
+#include <corsika/detail/framework/utility/QuarticSolver.inl>
diff --git a/corsika/media/BaseExponential.hpp b/corsika/media/BaseExponential.hpp
index b0afa123a8595864f77d3a15651402ed09f1b3af..b480b8e2c4207fcf3577489996d8cb0a4254740f 100644
--- a/corsika/media/BaseExponential.hpp
+++ b/corsika/media/BaseExponential.hpp
@@ -12,7 +12,7 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
#include <limits>
@@ -25,10 +25,6 @@ namespace corsika {
template <typename TDerived>
class BaseExponential {
protected:
- units::si::MassDensityType const rho0_;
- units::si::LengthType const lambda_;
- units::si::InverseLengthType const invLambda_;
- Point const point_;
auto const& getImplementation() const;
@@ -46,9 +42,8 @@ namespace corsika {
* \f]
*/
// clang-format on
- units::si::GrammageType getIntegratedGrammage(
- Trajectory<Line> const& line, units::si::LengthType vL,
- Vector<units::si::dimensionless_d> const& axis) const;
+ GrammageType getIntegratedGrammage(setup::Trajectory const& line, LengthType vL,
+ DirectionVector const& axis) const;
// clang-format off
/**
@@ -68,14 +63,19 @@ namespace corsika {
* \f]
*/
// clang-format on
- units::si::LengthType getArclengthFromGrammage(
- Trajectory<Line> const& line, units::si::GrammageType grammage,
- Vector<units::si::dimensionless_d> const& axis) const;
+ LengthType getArclengthFromGrammage(setup::Trajectory const& line,
+ GrammageType grammage,
+ DirectionVector const& axis) const;
public:
- BaseExponential(Point const& point, units::si::MassDensityType rho0,
- units::si::LengthType lambda);
+ BaseExponential(Point const& point, MassDensityType rho0, LengthType lambda);
+ private:
+ MassDensityType const rho0_;
+ LengthType const lambda_;
+ InverseLengthType const invLambda_;
+ Point const point_;
+
}; // class BaseExponential
} // namespace corsika
diff --git a/corsika/media/DensityFunction.hpp b/corsika/media/DensityFunction.hpp
index b22b54961a0aaa821227cc42edf5f803ac923866..6242b3d787ec507a0ac7db7a2e3e0d3b225d6813 100644
--- a/corsika/media/DensityFunction.hpp
+++ b/corsika/media/DensityFunction.hpp
@@ -10,7 +10,6 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/LinearApproximationIntegrator.hpp>
namespace corsika {
@@ -27,7 +26,7 @@ namespace corsika {
DensityFunction(TDerivableRho rho)
: rho_(rho) {}
- MassDensityType evaluateAt(corsika::Point const& p) const { return rho_(p); }
+ MassDensityType evaluateAt(Point const& p) const { return rho_(p); }
};
} // namespace corsika
diff --git a/corsika/media/FlatExponential.hpp b/corsika/media/FlatExponential.hpp
index 01f16105cc26a1eb38e8e75f3efe027eed4f7ccd..2e21772d4cd1a104f133ece3e73b5fb55a936d60 100644
--- a/corsika/media/FlatExponential.hpp
+++ b/corsika/media/FlatExponential.hpp
@@ -11,9 +11,9 @@ n/*
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/BaseExponential.hpp>
#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
namespace corsika {
@@ -30,10 +30,7 @@ namespace corsika {
// clang-format on
template <typename T>
class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
- Vector<dimensionless_d> const axis_;
- NuclearComposition const nuclComp_;
-
- using Base = BaseExponential<FlatExponential<T>>;
+ using base_type = BaseExponential<FlatExponential<T>>;
public:
FlatExponential(Point const& point, Vector<dimensionless_d> const& axis,
@@ -44,10 +41,15 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
- GrammageType getIntegratedGrammage(Trajectory<Line> const& line, LengthType to) const;
+ GrammageType getIntegratedGrammage(setup::Trajectory const& line,
+ LengthType to) const;
- LengthType getArclengthFromGrammage(Trajectory<Line> const& line,
+ LengthType getArclengthFromGrammage(setup::Trajectory const& line,
GrammageType grammage) const;
+
+ private:
+ DirectionVector const axis_;
+ NuclearComposition const nuclComp_;
};
} // namespace corsika
diff --git a/corsika/media/HomogeneousMedium.hpp b/corsika/media/HomogeneousMedium.hpp
index c0e0258166accdd65c6f806f02243f37599f88e2..a659be2288aef59492c9ab84922e6495c01a46d6 100644
--- a/corsika/media/HomogeneousMedium.hpp
+++ b/corsika/media/HomogeneousMedium.hpp
@@ -11,9 +11,10 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
/**
* a homogeneous medium
*/
@@ -30,10 +31,10 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
- GrammageType getIntegratedGrammage(Trajectory<Line> const&,
+ GrammageType getIntegratedGrammage(setup::Trajectory const&,
LengthType to) const override;
- LengthType getArclengthFromGrammage(Trajectory<Line> const&,
+ LengthType getArclengthFromGrammage(setup::Trajectory const&,
GrammageType grammage) const override;
private:
diff --git a/corsika/media/IEmpty.hpp b/corsika/media/IEmpty.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1912f90fc5e09ce403ef9c4e33c163e63f96309f
--- /dev/null
+++ b/corsika/media/IEmpty.hpp
@@ -0,0 +1,42 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+namespace corsika {
+
+ /**
+ * Intended for usage as default template argument for environments
+ * with no properties. For now, the ArclengthFromGrammage is
+ * mandatory, since it is used even in the most simple Cascade code.
+ *
+ * - IEmpty is the interface definition.
+ * - Empty<IEmpty> is a possible model implementation
+ *
+ **/
+
+ class IEmpty {
+ public:
+ virtual LengthType getArclengthFromGrammage(setup::Trajectory const&,
+ GrammageType) const = 0;
+
+ virtual ~IEmpty() {}
+ };
+
+ template <typename TModel = IEmpty>
+ class Empty : public TModel {
+ public:
+ LengthType getArclengthFromGrammage(setup::Trajectory const&, GrammageType) const {
+ return 0. * meter;
+ }
+ };
+
+} // namespace corsika
diff --git a/corsika/media/IMediumModel.hpp b/corsika/media/IMediumModel.hpp
index 7fcaafd12a220376a9810a26960970349c00205f..93ad9ccf58520a2a9151e43335a829ad3e4e213b 100644
--- a/corsika/media/IMediumModel.hpp
+++ b/corsika/media/IMediumModel.hpp
@@ -11,8 +11,9 @@
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
namespace corsika {
@@ -24,10 +25,10 @@ namespace corsika {
// todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
// approach; for now, only lines are supported
- virtual GrammageType getIntegratedGrammage(Trajectory<Line> const&,
+ virtual GrammageType getIntegratedGrammage(setup::Trajectory const&,
LengthType) const = 0;
- virtual LengthType getArclengthFromGrammage(Trajectory<Line> const&,
+ virtual LengthType getArclengthFromGrammage(setup::Trajectory const&,
GrammageType) const = 0;
virtual NuclearComposition const& getNuclearComposition() const = 0;
diff --git a/corsika/media/InhomogeneousMedium.hpp b/corsika/media/InhomogeneousMedium.hpp
index 615fe891065aa132cbd09ea97812d3737a0eaecf..7c9f664aaf589a311de7207d79a214125c8cdd15 100644
--- a/corsika/media/InhomogeneousMedium.hpp
+++ b/corsika/media/InhomogeneousMedium.hpp
@@ -11,8 +11,8 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
/**
* A general inhomogeneous medium. The mass density distribution TDensityFunction must be
@@ -32,10 +32,10 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
- GrammageType getIntegratedGrammage(Trajectory<Line> const& line,
+ GrammageType getIntegratedGrammage(setup::Trajectory const& line,
LengthType to) const override;
- LengthType getArclengthFromGrammage(Trajectory<Line> const& pLine,
+ LengthType getArclengthFromGrammage(setup::Trajectory const& pLine,
GrammageType grammage) const override;
private:
diff --git a/corsika/media/LinearApproximationIntegrator.hpp b/corsika/media/LinearApproximationIntegrator.hpp
index e14e0e9b38f72d05e16fe0af6c5f0da489717421..dfd9f0037476ea9d4e8dc8ed5ae3be6acccf8943 100644
--- a/corsika/media/LinearApproximationIntegrator.hpp
+++ b/corsika/media/LinearApproximationIntegrator.hpp
@@ -11,7 +11,7 @@
#include <limits>
#include <corsika/framework/geometry/Line.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
namespace corsika {
@@ -20,12 +20,12 @@ namespace corsika {
auto const& getImplementation() const;
public:
- auto getIntegrateGrammage(Trajectory<Line> const& line, LengthType length) const;
+ auto getIntegrateGrammage(setup::Trajectory const& line, LengthType length) const;
- auto getArclengthFromGrammage(Trajectory<Line> const& line,
+ auto getArclengthFromGrammage(setup::Trajectory const& line,
GrammageType grammage) const;
- auto getMaximumLength(Trajectory<Line> const& line,
+ auto getMaximumLength(setup::Trajectory const& line,
[[maybe_unused]] double relError) const;
};
diff --git a/corsika/media/ShowerAxis.hpp b/corsika/media/ShowerAxis.hpp
index a9cc58b086bbaec36743590b9063a5927f8622f5..6b2b49751b7f5aeb1e6a5b64ff6fc2699df9f2c4 100644
--- a/corsika/media/ShowerAxis.hpp
+++ b/corsika/media/ShowerAxis.hpp
@@ -45,7 +45,7 @@ namespace corsika {
public:
template <typename TEnvModel>
ShowerAxis(Point const& pStart, Vector<length_d> const& length,
- Environment<TEnvModel> const& env, int steps = 10'000);
+ Environment<TEnvModel> const& env, bool doThrow = false, int steps = 10'000);
LengthType getSteplength() const;
@@ -57,15 +57,16 @@ namespace corsika {
GrammageType getX(LengthType) const;
- Vector<dimensionless_d> const& getDirection() const;
+ DirectionVector const& getDirection() const;
Point const& getStart() const;
private:
Point const pointStart_;
Vector<length_d> const length_;
+ bool throw_ = false;
LengthType const max_length_, steplength_;
- Vector<dimensionless_d> const axis_normalized_;
+ DirectionVector const axis_normalized_;
std::vector<GrammageType> X_;
// for storing the lengths corresponding to equidistant X values
diff --git a/corsika/media/SlidingPlanarExponential.hpp b/corsika/media/SlidingPlanarExponential.hpp
index 66e51eea7d7c2659dc401d2a8ba4cd005598af1e..9060dc8765abc90e05147503462f7a9b1de5afe9 100644
--- a/corsika/media/SlidingPlanarExponential.hpp
+++ b/corsika/media/SlidingPlanarExponential.hpp
@@ -12,10 +12,10 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
namespace corsika {
@@ -46,10 +46,10 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
- GrammageType getIntegratedGrammage(Trajectory<Line> const& line,
+ GrammageType getIntegratedGrammage(setup::Trajectory const& line,
LengthType l) const override;
- LengthType getArclengthFromGrammage(Trajectory<Line> const& line,
+ LengthType getArclengthFromGrammage(setup::Trajectory const& line,
GrammageType grammage) const override;
private:
diff --git a/corsika/media/VolumeTreeNode.hpp b/corsika/media/VolumeTreeNode.hpp
index 98ba1d4422a42c471a40a45d6756591e18b8a7f1..110d17458e170442853c298d74801cecd32bf0e3 100644
--- a/corsika/media/VolumeTreeNode.hpp
+++ b/corsika/media/VolumeTreeNode.hpp
@@ -9,21 +9,19 @@
#pragma once
#include <corsika/framework/geometry/IVolume.hpp>
-#include <corsika/media/IMediumModel.hpp>
+#include <corsika/environment/IEmpty.hpp>
#include <memory>
#include <vector>
namespace corsika {
- class Empty {}; //<! intended for usage as default template argument
-
- template <typename TModelProperties = Empty>
+ template <typename TModelProperties = IEmpty>
class VolumeTreeNode {
public:
using IModelProperties = TModelProperties;
using VTN_type = VolumeTreeNode<IModelProperties>;
- using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
+ using VTNUPtr = std::unique_ptr<VTN_type>;
using IMPSharedPtr = std::shared_ptr<IModelProperties>;
using VolUPtr = std::unique_ptr<IVolume>;
@@ -53,7 +51,7 @@ namespace corsika {
inline void excludeOverlapWith(VTNUPtr const& pNode);
- inline auto* getParent() const { return parentNode_; };
+ inline VTN_type* getParent() const { return parentNode_; };
inline auto const& getChildNodes() const { return childNodes_; }
@@ -67,23 +65,21 @@ namespace corsika {
template <typename ModelProperties, typename... Args>
inline auto setModelProperties(Args&&... args) {
- static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
- "unusable model properties type provided");
-
+ //static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
+ // "unusable type provided");
modelProperties_ = std::make_shared<ModelProperties>(std::forward<Args>(args)...);
return modelProperties_;
}
inline void setModelProperties(IMPSharedPtr ptr) { modelProperties_ = ptr; }
- /*
- template <class MediumType, typename... Args>
- static auto createMedium(Args&&... args);
+ //template <class MediumType, typename... Args>
+ //static auto createMedium(Args&&... args);
private:
std::vector<VTNUPtr> childNodes_;
- std::vector<VolumeTreeNode<IModelProperties> const*> excludedNodes_;
- VolumeTreeNode<IModelProperties> const* parentNode_ = nullptr;
+ std::vector<VTN_type const*> excludedNodes_;
+ VTN_type const* parentNode_ = nullptr;
VolUPtr geoVolume_;
IMPSharedPtr modelProperties_;
};
diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index cd645194c34de557ac58f90039744bc719433bd2..ea87ea2e49602f2869e2ef80d60f404897e89311 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -27,8 +27,8 @@ namespace corsika {
class ObservationPlane : public ContinuousProcess<ObservationPlane> {
public:
- ObservationPlane(Plane const&, std::string const&, bool = true);
- void Init() {}
+ ObservationPlane(Plane const&, DirectionVector const&, std::string const&,
+ bool = true);
ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
corsika::setup::Trajectory& vTrajectory);
@@ -46,6 +46,9 @@ namespace corsika {
bool const deleteOnHit_;
HEPEnergyType energy_ground_;
unsigned int count_ground_;
+ DirectionVector const xAxis_;
+ DirectionVector const yAxis_;
+
};
} // namespace corsika
diff --git a/corsika/modules/TrackingLine.hpp b/corsika/modules/TrackingLine.hpp
index c49c4bbdc98d6e212b24aa01a7296e2d0a75d650..7cfb386730b706ee7e42c0915f60a96066def150 100644
--- a/corsika/modules/TrackingLine.hpp
+++ b/corsika/modules/TrackingLine.hpp
@@ -12,108 +12,94 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
+
+#include <corsika/framework/logging/Logging.hpp>
+#include <corsika/modules/tracking/Intersect.hpp>
-#include <optional>
#include <type_traits>
#include <utility>
namespace corsika::tracking_line {
- std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(Line const&,
- Sphere const&);
-
- TimeType TimeOfIntersection(Line const&, Plane const&);
+ /**
+ * \class Tracking
+ *
+ *
+ *
+ **/
- class TrackingLine {
+ class Tracking : public Intersect<Tracking> {
public:
- TrackingLine(){};
-
- template <typename Particle> // was Stack previously, and argument was
- // Stack::StackIterator
- auto getTrack(Particle const& p) {
- Vector<SpeedType::dimension_type> const velocity =
- p.getMomentum() / p.getEnergy() * constants::c;
-
- auto const currentPosition = p.getPosition();
- std::cout << "TrackingLine pid: " << p.getPID()
- << " , E = " << p.getEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine pos: " << currentPosition.getCoordinates() << std::endl;
- std::cout << "TrackingLine E: " << p.getEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine p: " << p.getMomentum().getComponents() / 1_GeV
- << " GeV " << std::endl;
- std::cout << "TrackingLine v: " << velocity.getComponents() << std::endl;
-
- // to do: include effect of magnetic field
- Line line(currentPosition, velocity);
-
- auto const* currentLogicalVolumeNode = p.getNode();
- auto const numericallyInside =
- currentLogicalVolumeNode->getVolume().isInside(currentPosition);
-
- std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
-
- auto const& children = currentLogicalVolumeNode->getChildNodes();
- auto const& excluded = currentLogicalVolumeNode->getExcludedNodes();
-
- std::vector<std::pair<TimeType, decltype(p.getNode())>> intersections;
-
- // for entering from outside
- auto addIfIntersects = [&](auto const& vtn) {
- auto const& volume = vtn.getVolume();
- auto const& sphere = dynamic_cast<Sphere const&>(
- volume); // for the moment we are a bit bold here and assume
- // everything is a sphere, crashes with exception if not
-
- if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
- auto const [t1, t2] = *opt;
- std::cout << "intersection times: " << t1 / 1_s << "; "
- << t2 / 1_s
- // << " " << vtn.getModelProperties().getName()
- << std::endl;
- if (t1.magnitude() > 0)
- intersections.emplace_back(t1, &vtn);
- else if (t2.magnitude() > 0)
- std::cout << "inside other volume" << std::endl;
- }
- };
-
- for (auto const& child : children) { addIfIntersects(*child); }
- for (auto const* ex : excluded) { addIfIntersects(*ex); }
-
- {
- auto const& sphere =
- dynamic_cast<Sphere const&>(currentLogicalVolumeNode->getVolume());
- // for the moment we are a bit bold here and assume
- // everything is a sphere, crashes with exception if not
- [[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
- [[maybe_unused]] auto dummy_t1 = t1;
- intersections.emplace_back(t2, currentLogicalVolumeNode->getParent());
- }
-
- auto const minIter = std::min_element(
- intersections.cbegin(), intersections.cend(),
- [](auto const& a, auto const& b) { return a.first < b.first; });
-
- TimeType min;
+ template <typename TParticle>
+ auto getTrack(TParticle const& particle) {
+ VelocityVector const initialVelocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+
+ auto const initialPosition = particle.getPosition();
+ CORSIKA_LOG_DEBUG(
+ "Tracking pid: {}"
+ " , E = {} GeV",
+ particle.getPID(), particle.getEnergy() / 1_GeV);
+ CORSIKA_LOG_DEBUG("Tracking pos: {}", initialPosition.getCoordinates());
+ CORSIKA_LOG_DEBUG("Tracking E: {} GeV", particle.getEnergy() / 1_GeV);
+ CORSIKA_LOG_DEBUG("Tracking p: {} GeV", particle.getMomentum().getComponents() / 1_GeV);
+ CORSIKA_LOG_DEBUG("Tracking v: {} ", initialVelocity.getComponents());
+
+ // traverse the environment volume tree and find next
+ // intersection
+ auto [minTime, minNode] = Intersect<Tracking>::nextIntersect(particle);
+
+ return std::make_tuple(LineTrajectory(Line(initialPosition, initialVelocity),
+ minTime), // trajectory
+ minNode); // next volume node
+ }
- if (minIter == intersections.cend()) {
- min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
- // to handle the numerics properly
- throw std::runtime_error("no intersection with anything!");
- } else {
- min = minIter->first;
+ template <typename TParticle, typename TMedium>
+ static Intersections intersect(TParticle const& particle, Sphere const& sphere,
+ TMedium const&) {
+ auto const delta = particle.getPosition() - sphere.getCenter();
+ auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
+ auto const vSqNorm = velocity.squaredNorm();
+ auto const R = sphere.getRadius();
+
+ auto const vDotDelta = velocity.dot(delta);
+ auto const discriminant =
+ vDotDelta * vDotDelta - vSqNorm * (delta.getSquaredNorm() - R * R);
+
+ if (discriminant.magnitude() > 0) {
+ auto const sqDisc = sqrt(discriminant);
+ auto const invDenom = 1 / vSqNorm;
+ return Intersections((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
}
+ return Intersections();
+ }
- std::cout << " t-intersect: "
- << min
- // << " " << minIter->second->getModelProperties().getName()
- << std::endl;
+ template <typename TParticle, typename TBaseNodeType>
+ static Intersections intersect(TParticle const& particle,
+ TBaseNodeType const& volumeNode) {
+ Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
+ if (sphere) {
+ return Intersect(particle, *sphere, volumeNode.getModelProperties());
+ }
+ throw std::runtime_error(
+ "The Volume type provided is not supported in Intersect(particle, node)");
+ }
- return std::make_tuple(Trajectory<Line>(line, min), velocity.getNorm() * min,
- minIter->second);
+ template <typename TParticle, typename TMedium>
+ static Intersections intersect(TParticle const& particle, Plane const& plane,
+ TMedium const&) {
+ auto const delta = plane.getCenter() - particle.getPosition();
+ auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
+ auto const n = plane.getNormal();
+ auto const c = n.dot(velocity);
+
+ return Intersections(
+ c.magnitude() == 0 ? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
+ : n.dot(delta) / c);
}
};
diff --git a/corsika/modules/tracking/Intersect.hpp b/corsika/modules/tracking/Intersect.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..73593abd0fec8389a13e3779b55c73fad4b3d681
--- /dev/null
+++ b/corsika/modules/tracking/Intersect.hpp
@@ -0,0 +1,134 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/logging/Logging.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
+
+#include <limits>
+
+namespace corsika {
+
+ /**
+ *
+ * This is a CRTP class to provide a generic volume-tree
+ * intersection for the purpose of tracking.
+ *
+ * It return the closest distance in time to the next geometric
+ * intersection, as well as a pointer to the corresponding new
+ * volume.
+ *
+ * User may provide an optional global step-length limit as
+ * parameter. This may be needd for (simpler) algorithms in magnetic
+ * fields, where tracking errors grow linearly with step-length.
+ * Obviously, in the case of the step-length limit, the returend
+ * "next" volume is just the current one.
+ *
+ **/
+
+ template <typename TDerived>
+ class Intersect {
+
+ protected:
+ template <typename TParticle>
+ auto nextIntersect(
+ const TParticle& particle,
+ TimeType step_limit = std::numeric_limits<TimeType::value_type>::infinity() *
+ second) const {
+
+ typedef
+ typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
+ node_type& volumeNode =
+ *particle.getNode(); // current "logical" node, from previous tracking step
+ CORSIKA_LOG_DEBUG("volumeNode={}, numericallyInside={} ", fmt::ptr(&volumeNode),
+ volumeNode.getVolume().isInside(particle.getPosition()));
+
+ // start values:
+ TimeType minTime = step_limit;
+ node_type* minNode = &volumeNode;
+
+ // determine the first geometric collision with any other Volume boundary
+
+ // first check, where we leave the current volume
+ // this assumes our convention, that all Volume primitives must be convex
+ // thus, the last entry is always the exit point
+ const Intersections time_intersections_curr =
+ TDerived::intersect(particle, volumeNode);
+ CORSIKA_LOG_TRACE("curr node {}, parent node {}, hasIntersections={} ",
+ fmt::ptr(&volumeNode), fmt::ptr(volumeNode.getParent()),
+ time_intersections_curr.hasIntersections());
+ if (time_intersections_curr.hasIntersections()) {
+ CORSIKA_LOG_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
+ // entirely contained by currentLogicalVolumeNode
+ for (const auto& node : volumeNode.getChildNodes()) {
+
+ const Intersections time_intersections = TDerived::intersect(particle, *node);
+ if (!time_intersections.hasIntersections()) { continue; }
+ CORSIKA_LOG_DEBUG("intersection times with child volume {} : enter {} s, exit {} s",
+ fmt::ptr(node), time_intersections.getEntry() / 1_s,
+ time_intersections.getExit() / 1_s);
+
+ const auto t_entry = time_intersections.getEntry();
+ const auto t_exit = time_intersections.getExit();
+ CORSIKA_LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
+ 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", 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();
+ }
+ }
+
+ // these are volumes from the previous tree-level that are cut-out partly from the
+ // current volume
+ for (node_type* node : volumeNode.getExcludedNodes()) {
+
+ const Intersections time_intersections = TDerived::intersect(particle, *node);
+ if (!time_intersections.hasIntersections()) { continue; }
+ CORSIKA_LOG_DEBUG("intersection times with exclusion volume {} : enter {} s, exit {} s",
+ fmt::ptr(node), time_intersections.getEntry() / 1_s,
+ time_intersections.getExit() / 1_s);
+ const auto t_entry = time_intersections.getEntry();
+ const auto t_exit = time_intersections.getExit();
+ CORSIKA_LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
+ 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
+ minTime = t_entry;
+ minNode = node;
+ }
+ }
+ CORSIKA_LOG_TRACE("t-intersect: {}, node {} ", minTime, fmt::ptr(minNode));
+ return std::make_tuple(minTime, minNode);
+ }
+ };
+
+} // namespace corsika::tracking
diff --git a/corsika/setup/SetupTrajectory.hpp b/corsika/setup/SetupTrajectory.hpp
index 8bed53ce8cde7cf5e3f43b94418b3a5eb8839984..254e8d25a9b10680a965bab5ebc2bb3bdbe97e70 100644
--- a/corsika/setup/SetupTrajectory.hpp
+++ b/corsika/setup/SetupTrajectory.hpp
@@ -8,12 +8,51 @@ n/*
#pragma once
+#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Helix.hpp>
#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/modules/TrackingLine.hpp>
+//#include <corsika/modules/TrackingCurved.hpp> // simple leap-frog implementation
+//#include <corsika/modules/TrackingLeapFrog.hpp> // more complete leap-frog
+// implementation
+
namespace corsika::setup {
+ /**
+ \file SetupTrajectory.hpp
+
+ Note/Warning: Tracking and Trajectory must fit together !
+
+ tracking_leapfrog_curved::Tracking is the result of the Bachelor
+ thesis of Andre Schmidt, KIT. This is a leap-frog algorithm with
+ an analytical, precise calculation of volume intersections. This
+ algorithm needs a LeapFrogTrajectory.
+
+ tracking_leapfrog_straight::Tracking is a more simple and direct
+ leap-frog implementation. The two halve steps are coded explicitly
+ as two straight segments. Intersections with other volumes are
+ calculate only on the straight segments. This algorithm is based
+ on LineTrajectory.
+
+ tracking_line::Tracking is a pure straight tracker. It is based on
+ LineTrajectory.
+ */
+
+ /**
+ The default tracking algorithm.
+ */
+
+ // typedef corsika::process::tracking_leapfrog_curved::Tracking Tracking;
+ // typedef corsika::process::tracking_leapfrog_straight::Tracking Tracking;
+ typedef corsika::tracking_line::Tracking Tracking;
+
+ /**
+ The default trajectory.
+ */
/// definition of Trajectory base class, to be used in tracking and cascades
- typedef corsika::Trajectory<corsika::Line> Trajectory;
+ typedef LineTrajectory Trajectory;
+ // typedef corsika::geometry::LeapFrogTrajectory Trajectory;
} // namespace corsika::setup
diff --git a/corsika/stack/history/HistoryObservationPlane.cpp b/corsika/stack/history/HistoryObservationPlane.cpp
deleted file mode 100644
index b1425f18be3ff20fe68f8ffbaf2f6ff9244d6e78..0000000000000000000000000000000000000000
--- a/corsika/stack/history/HistoryObservationPlane.cpp
+++ /dev/null
@@ -1,84 +0,0 @@
-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/logging/Logging.h>
-#include <corsika/stack/history/HistoryObservationPlane.hpp>
-
-#include <boost/histogram/ostream.hpp>
-
-#include <iomanip>
-#include <iostream>
-
-using namespace corsika::units::si;
-using namespace corsika::history;
-using namespace corsika;
-
-HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
- geometry::Plane const& obsPlane,
- bool deleteOnHit)
- : stack_{stack}
- , plane_{obsPlane}
- , deleteOnHit_{deleteOnHit} {}
-
-corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
- setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
- TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
- trajectory.GetV0().dot(plane_.GetNormal());
-
- if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
-
- if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
- return process::EProcessReturn::eOk;
- }
-
- C8LOG_DEBUG(fmt::format("HistoryObservationPlane: Particle detected: pid={}",
- particle.GetPID()));
-
- auto const pid = particle.GetPID();
- if (particles::IsMuon(pid)) { fillHistoryHistogram(particle); }
-
- if (deleteOnHit_) {
- return process::EProcessReturn::eParticleAbsorbed;
- } else {
- return process::EProcessReturn::eOk;
- }
-}
-
-LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
- setup::Trajectory const& trajectory) {
- TimeType const timeOfIntersection =
- (plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
- trajectory.GetV0().dot(plane_.GetNormal());
-
- if (timeOfIntersection < TimeType::zero()) {
- return std::numeric_limits<double>::infinity() * 1_m;
- }
-
- auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
- return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
-}
-
-void HistoryObservationPlane::fillHistoryHistogram(
- setup::Stack::ParticleType const& muon) {
- double const muon_energy = muon.GetEnergy() / 1_GeV;
-
- int genctr{0};
- Event const* event = muon.GetEvent().get();
- while (event) {
- auto const projectile = stack_.cfirst() + event->projectileIndex();
- if (event->eventType() == EventType::Interaction) {
- genctr++;
- double const projEnergy = projectile.GetEnergy() / 1_GeV;
- int const pdg = static_cast<int>(particles::GetPDG(projectile.GetPID()));
-
- histogram_(muon_energy, projEnergy, pdg);
- }
- event = event->parentEvent().get(); // projectile.GetEvent().get();
- }
-}
diff --git a/corsika/stack/history/HistoryObservationPlane.hpp b/corsika/stack/history/HistoryObservationPlane.hpp
index fb67ca7f9b8897cb6d6f2b11c6ce02101a5ea653..9c97f8722218f24449ef32538be81ea2ac5c1056 100644
--- a/corsika/stack/history/HistoryObservationPlane.hpp
+++ b/corsika/stack/history/HistoryObservationPlane.hpp
@@ -18,43 +18,33 @@
#include <functional>
+// the detail namespace: here the histrograms are defined
+//! \todo add options/parameters to 'detail::hist_factory()'
+#include <corsika/detail/stack/history/HistoryObservationPlane.hpp>
+
namespace corsika::history {
- namespace detail {
- inline auto hist_factory() {
- namespace bh = boost::histogram;
- namespace bha = bh::axis;
- auto h = bh::make_histogram(
- bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11, "muon energy"},
- bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11,
- "projectile energy"},
- bha::category<int, bh::use_default, bha::option::growth_t>{
- {211, -211, 2212, -2212}, "projectile PDG"});
- return h;
- }
- } // namespace detail
-
- class HistoryObservationPlane
- : public corsika::process::ContinuousProcess<HistoryObservationPlane> {
+
+ class HistoryObservationPlane : public ContinuousProcess<HistoryObservationPlane> {
public:
- HistoryObservationPlane(setup::Stack const&, geometry::Plane const&, bool = true);
+ HistoryObservationPlane(setup::Stack const&, Plane const&, bool = true);
- corsika::units::si::LengthType MaxStepLength(
- corsika::setup::Stack::ParticleType const&,
- corsika::setup::Trajectory const& vTrajectory);
+ LengthType getMaxStepLength(setup::Stack::particle_type const&,
+ setup::Trajectory const& vTrajectory);
- corsika::process::EProcessReturn DoContinuous(
- corsika::setup::Stack::ParticleType const& vParticle,
- corsika::setup::Trajectory const& vTrajectory);
+ ProcessReturn doContinuous(setup::Stack::particle_type const& vParticle,
+ setup::Trajectory const& vTrajectory);
auto const& histogram() const { return histogram_; }
private:
- void fillHistoryHistogram(setup::Stack::ParticleType const&);
+ void fillHistoryHistogram(setup::Stack::particle_type const&);
setup::Stack const& stack_;
- geometry::Plane const plane_;
+ Plane const plane_;
bool const deleteOnHit_;
decltype(detail::hist_factory()) histogram_ = detail::hist_factory();
};
} // namespace corsika::history
+
+#include <corsika/detail/stack/history/HistoryObservationPlane.inl>
diff --git a/tests/common/SetupEnvironment.hpp b/tests/common/SetupEnvironment.hpp
index 6c6ee7677ab3a724ccbd8204df851a411aaa119f..26399e247051b7a83b90b1cc668a27f5fdb766eb 100644
--- a/tests/common/SetupEnvironment.hpp
+++ b/tests/common/SetupEnvironment.hpp
@@ -23,8 +23,10 @@
namespace corsika::setup::testing {
inline std::tuple<std::unique_ptr<setup::Environment>, CoordinateSystem const*,
- setup::Environment::BaseNodeType const*>
- setup_environment(Code vTargetCode) {
+ setup::Environment::BaseNodeType*>
+ setupEnvironment(Code const vTargetCode,
+ MagneticFluxType const& BfieldZ = MagneticFluxType::zero()) {
+
auto env = std::make_unique<setup::Environment>();
auto& universe = *(env->getUniverse());
CoordinateSystemPtr const& cs = env->getCoordinateSystem();
@@ -32,8 +34,7 @@ namespace corsika::setup::testing {
/**
* our world is a sphere at 0,0,0 with R=infty
*/
- auto world = setup::Environment::createNode<Sphere>(
- Point{cs, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+ auto world = setup::Environment::createNode<Sphere>(Point{cs, 0_m, 0_m, 0_m}, 100_km);
/**
* construct suited environment medium model:
@@ -42,10 +43,10 @@ namespace corsika::setup::testing {
UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
world->setModelProperties<MyHomogeneousModel>(
- Medium::AirDry1Atm, Vector(cs, 0_T, 0_T, 1_T), 1_kg / (1_m * 1_m * 1_m),
+ Medium::AirDry1Atm, Vector(cs, 0_T, 0_T, BfieldZ), 1_kg / (1_m * 1_m * 1_m),
NuclearComposition(std::vector<Code>{vTargetCode}, std::vector<float>{1.}));
- setup::Environment::BaseNodeType const* nodePtr = world.get();
+ setup::Environment::BaseNodeType* nodePtr = world.get();
universe.addChild(std::move(world));
return std::make_tuple(std::move(env), &cs, nodePtr);
diff --git a/tests/common/SetupTestTrajectory.hpp b/tests/common/SetupTestTrajectory.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..edcf4a8554a7ce3bc8e861ce0413c39d35e34db1
--- /dev/null
+++ b/tests/common/SetupTestTrajectory.hpp
@@ -0,0 +1,36 @@
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Helix.hpp>
+#include <corsika/framework/geometry/Trajectory.hpp>
+
+#include <corsika/modules/TrackingLine.hpp>
+//#include <corsika/modules/TrackingCurved.hpp> // simple leap-frog implementation
+//#include <corsika/modules/TrackingLeapFrog.hpp> // more complete leap-frog
+// implementation
+
+namespace corsika::setup::testing {
+
+ template <typename TTrack>
+ TTrack make_track(Line const& line, TimeType const tEnd);
+
+ template <>
+ inline LineTrajectory make_track<LineTrajectory>(Line const& line,
+ TimeType const tEnd) {
+ return LineTrajectory(line, tEnd);
+ }
+
+ /*
+ template <>
+ inline LeapFrogTrajectory make_track<LeapFrogTrajectory>(Line const& line,
+ TimeType const tEnd) {
+
+ auto const k = square(0_m) / (square(1_s) * 1_V);
+ return LeapFrogTrajectory(
+ line.getStartPoint(), line.getVelocity(),
+ MagneticFieldVector{line.getStartPoint().getCoordinateSystem(), 0_T, 0_T, 0_T},
+ k, tEnd);
+ }
+ */
+} // namespace corsika::setup::testing
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index 67fd3ae3bb442f5a7a7e5b296db951be44ce9bd1..6bc0c607927e608c22f8bc893cde178d315ac2cb 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -118,7 +118,7 @@ TEST_CASE("Cascade", "[Cascade]") {
auto env = make_dummy_env();
auto const& rootCS = env.getCoordinateSystem();
- tracking_line::TrackingLine tracking;
+ tracking_line::Tracking tracking;
StackInspector<TestCascadeStack> stackInspect(1, true, E0);
NullModel nullModel;
@@ -133,7 +133,7 @@ TEST_CASE("Cascade", "[Cascade]") {
MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
- Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack,
+ Cascade<tracking_line::Tracking, decltype(sequence), TestCascadeStack,
TestCascadeStackView>
EAS(env, tracking, sequence, stack);
diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp
index d1bc4bdc5984e265a0ebd87f31b27314ded6506e..dd860f2df535e75871c9e1bba0b0dfc5db20905b 100644
--- a/tests/framework/testGeometry.cpp
+++ b/tests/framework/testGeometry.cpp
@@ -309,8 +309,7 @@ TEST_CASE("Trajectories") {
Point r0(rootCS, {0_m, 0_m, 0_m});
SECTION("Line") {
- Vector<SpeedType::dimension_type> v0(rootCS,
- {3_m / second, 0_m / second, 0_m / second});
+ VelocityVector v0(rootCS, {3_m / second, 0_m / second, 0_m / second});
Line const line(r0, v0);
CHECK(
@@ -329,22 +328,18 @@ TEST_CASE("Trajectories") {
.magnitude() == Approx(0).margin(absMargin));
auto const t = 1_s;
- Trajectory<Line> base(line, t);
+ LineTrajectory base(line, t);
CHECK(line.getPosition(t).getCoordinates() == base.getPosition(1.).getCoordinates());
- CHECK(base.getArcLength(1_s, 2_s) / 1_m == Approx(3));
-
- CHECK((base.getNormalizedDirection().getComponents(rootCS) -
+ CHECK((base.getDirection(0).getComponents(rootCS) -
QuantityVector<dimensionless_d>{1, 0, 0})
.getNorm() == Approx(0).margin(absMargin));
}
SECTION("Helix") {
- Vector<SpeedType::dimension_type> const vPar(
- rootCS, {0_m / second, 0_m / second, 4_m / second});
+ VelocityVector const vPar(rootCS, {0_m / second, 0_m / second, 4_m / second});
- Vector<SpeedType::dimension_type> const vPerp(
- rootCS, {3_m / second, 0_m / second, 0_m / second});
+ VelocityVector const vPerp(rootCS, {3_m / second, 0_m / second, 0_m / second});
auto const T = 1_s;
auto const omegaC = 2 * M_PI / T;
@@ -365,11 +360,5 @@ TEST_CASE("Trajectories") {
helix.getPositionFromArclength(helix.getArcLength(0_s, 7_s)))
.getNorm()
.magnitude() == Approx(0).margin(absMargin));
-
- auto const t = 1234_s;
- Trajectory<Helix> const base(helix, t);
- CHECK(helix.getPosition(t).getCoordinates() == base.getPosition(1.).getCoordinates());
-
- CHECK(base.getArcLength(0_s, 1_s) / 1_m == Approx(5));
}
}