From 04211f1252d843df4d8d525ff76014d28a4e7752 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 5 Jan 2021 18:05:44 +0100
Subject: [PATCH] various improvements
---
corsika/corsika.hpp | 32 ++-
corsika/detail/framework/geometry/Sphere.inl | 2 +-
.../detail/framework/geometry/Trajectory.inl | 156 -----------
corsika/detail/media/Universe.inl | 4 +-
corsika/detail/media/VolumeTreeNode.inl | 15 +-
.../detail/modules/tracking/TrackingLine.inl | 59 ----
.../modules/tracking/TrackingStraight.inl | 100 +++++++
.../detail/stack/NuclearStackExtension.inl | 2 +-
corsika/detail/stack/SuperStupidStack.inl | 93 -------
corsika/framework/core/ParticleProperties.hpp | 5 +
corsika/framework/geometry/IVolume.hpp | 2 +-
corsika/framework/geometry/Sphere.hpp | 2 +-
corsika/framework/geometry/Trajectory.hpp | 216 ---------------
corsika/framework/stack/Stack.hpp | 24 ++
corsika/media/Universe.hpp | 2 +-
corsika/media/VolumeTreeNode.hpp | 4 +-
corsika/modules/tracking/Intersect.hpp | 104 +------
.../tracking/TrackingLeapFrogCurved.hpp | 225 +--------------
.../tracking/TrackingLeapFrogStraight.hpp | 167 +-----------
corsika/modules/tracking/TrackingStraight.hpp | 76 +-----
corsika/setup/SetupTrajectory.hpp | 5 +-
corsika/stack/NuclearStackExtension.hpp | 8 +-
corsika/stack/SuperStupidStack.hpp | 176 ------------
examples/geometry_example.cpp | 8 +-
examples/stack_example.cpp | 8 +-
tests/common/SetupTestTrajectory.hpp | 9 +-
tests/common/testTestStack.hpp | 19 +-
tests/framework/testCascade.cpp | 2 +-
tests/framework/testGeometry.cpp | 11 +-
tests/framework/testRandom.cpp | 12 +-
tests/modules/CMakeLists.txt | 4 +-
tests/modules/testInteractionCounter.cpp | 32 +--
tests/modules/testQGSJetII.cpp | 44 +--
tests/modules/testStackInspector.cpp | 2 +-
tests/modules/testSwitchProcess.cpp | 258 ------------------
tests/stack/CMakeLists.txt | 2 +-
tests/stack/testNuclearStackExtension.cpp | 6 +-
...perStupidStack.cpp => testSimpleStack.cpp} | 8 +-
38 files changed, 303 insertions(+), 1601 deletions(-)
delete mode 100644 corsika/detail/framework/geometry/Trajectory.inl
delete mode 100644 corsika/detail/modules/tracking/TrackingLine.inl
create mode 100644 corsika/detail/modules/tracking/TrackingStraight.inl
delete mode 100644 corsika/detail/stack/SuperStupidStack.inl
delete mode 100644 corsika/framework/geometry/Trajectory.hpp
delete mode 100644 corsika/stack/SuperStupidStack.hpp
delete mode 100644 tests/modules/testSwitchProcess.cpp
rename tests/stack/{testSuperStupidStack.cpp => testSimpleStack.cpp} (91%)
diff --git a/corsika/corsika.hpp b/corsika/corsika.hpp
index 137cccc6b..36fa90b63 100644
--- a/corsika/corsika.hpp
+++ b/corsika/corsika.hpp
@@ -8,10 +8,32 @@
#pragma once
-// Usage:
-// to get the version X.YY.Z,
-// set CORSIKA_VERSION X0YY0Z
-//
+/**
+ @file corsika.hpp
+
+ The CORSIKA 8 air shower simulation framework.
+
+ @mainpage Technical documentation of the CORSIKA 8 software framework
+
+ Software documentatin and reference guide for the CORSIKA 8
+ software framework for air shower simulations. CORSIKA 8 is developed
+ at <a href="https://gitlab.ikp.kit.edu/AirShowerPhysics">https://gitlab.ikp.kit.edu</a>.
+ If you got the code from somewhere else, consider to switch to the official development
+ repository. If you want to report bugs, or want to suggest features or future
+ development, please submit an "issue" on this gitlab server. We only accept Issues and
+ discussion via our central development server https://gitlab.ikp.kit.edu.
+
+ Write to corsika-devel@lists.kit.edu, or even register yourself at
+ https://www.lists.kit.edu/sympa/info/corsika-devel to get in contact
+ with other developers.
+
+ For more information about the project, see @ref md_README.
+ **/
+
+/*! Usage:
+ * to get the version X.YY.Z,
+ * set CORSIKA_VERSION X0YY0Z
+ */
#define CORSIKA_VERSION 800000
/*! \def CORSIKA_MAJOR_VERSION
@@ -30,4 +52,4 @@
* \brief The preprocessor macro \p CORSIKA_PATCH_NUMBER encodes the
* patch number of the CORSIKA library.
*/
-#define CORSIKA_PATCH_NUMBER 0
\ No newline at end of file
+#define CORSIKA_PATCH_NUMBER 0
diff --git a/corsika/detail/framework/geometry/Sphere.inl b/corsika/detail/framework/geometry/Sphere.inl
index a793984dd..be1ba2774 100644
--- a/corsika/detail/framework/geometry/Sphere.inl
+++ b/corsika/detail/framework/geometry/Sphere.inl
@@ -13,7 +13,7 @@
namespace corsika {
- inline bool Sphere::isInside(Point const& p) const {
+ inline bool Sphere::contains(Point const& p) const {
return radius_ * radius_ > (center_ - p).getSquaredNorm();
}
diff --git a/corsika/detail/framework/geometry/Trajectory.inl b/corsika/detail/framework/geometry/Trajectory.inl
deleted file mode 100644
index ab7779c83..000000000
--- a/corsika/detail/framework/geometry/Trajectory.inl
+++ /dev/null
@@ -1,156 +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.
- */
-
-#pragma once
-
-#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 {
-
- 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_);
- }
-
- 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);
- }
-
- 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;
- }
-
- inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
- return initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
- }
-
- inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
- return getVelocity(u).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/media/Universe.inl b/corsika/detail/media/Universe.inl
index 6a2bdd6e9..b3f4d2541 100644
--- a/corsika/detail/media/Universe.inl
+++ b/corsika/detail/media/Universe.inl
@@ -16,6 +16,6 @@ namespace corsika {
: corsika::Sphere(Point{pCS, 0 * meter, 0 * meter, 0 * meter},
meter * std::numeric_limits<double>::infinity()) {}
- bool Universe::isInside(corsika::Point const&) const { return true; }
+ bool Universe::contains(corsika::Point const&) const { return true; }
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/detail/media/VolumeTreeNode.inl b/corsika/detail/media/VolumeTreeNode.inl
index 6f978f86b..bdae7cc12 100644
--- a/corsika/detail/media/VolumeTreeNode.inl
+++ b/corsika/detail/media/VolumeTreeNode.inl
@@ -13,18 +13,17 @@
namespace corsika {
- //! convenience function equivalent to Volume::isInside
template <typename IModelProperties>
- bool VolumeTreeNode<IModelProperties>::isInside(Point const& p) const {
- return geoVolume_->isInside(p);
+ bool VolumeTreeNode<IModelProperties>::contains(Point const& p) const {
+ return geoVolume_->contains(p);
}
template <typename IModelProperties>
inline VolumeTreeNode<IModelProperties> const*
- VolumeTreeNode<IModelProperties>::isExcluded(Point const& p) const {
+ VolumeTreeNode<IModelProperties>::excludes(Point const& p) const {
auto exclContainsIter =
std::find_if(excludedNodes_.cbegin(), excludedNodes_.cend(),
- [&](auto const& s) { return bool(s->isInside(p)); });
+ [&](auto const& s) { return bool(s->contains(p)); });
return exclContainsIter != excludedNodes_.cend() ? *exclContainsIter : nullptr;
}
@@ -35,14 +34,14 @@ namespace corsika {
template <typename IModelProperties>
VolumeTreeNode<IModelProperties> const*
VolumeTreeNode<IModelProperties>::getContainingNode(Point const& p) const {
- if (!isInside(p)) { return nullptr; }
+ if (!contains(p)) { return nullptr; }
if (auto const childContainsIter =
std::find_if(childNodes_.cbegin(), childNodes_.cend(),
- [&](auto const& s) { return bool(s->isInside(p)); });
+ [&](auto const& s) { return bool(s->contains(p)); });
childContainsIter == childNodes_.cend()) // not contained in any of the children
{
- if (auto const exclContainsIter = isExcluded(p)) // contained in any excluded nodes
+ if (auto const exclContainsIter = excludes(p)) // contained in any excluded nodes
{
return exclContainsIter->getContainingNode(p);
} else {
diff --git a/corsika/detail/modules/tracking/TrackingLine.inl b/corsika/detail/modules/tracking/TrackingLine.inl
deleted file mode 100644
index 52200c42d..000000000
--- a/corsika/detail/modules/tracking/TrackingLine.inl
+++ /dev/null
@@ -1,59 +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.
- */
-
-#pragma once
-
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/QuantityVector.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/media/Environment.hpp>
-
-#include <limits>
-#include <stdexcept>
-#include <utility>
-
-namespace corsika::tracking_line {
-
- std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(
- corsika::Line const& line, corsika::Sphere const& sphere) {
- auto const delta = line.getStartPoint() - sphere.getCenter();
- auto const v = line.getVelocity();
- auto const vSqNorm =
- v.getSquaredNorm(); // todo: get rid of this by having V0 normalized always
- auto const R = sphere.getRadius();
-
- auto const vDotDelta = v.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 std::make_pair((-vDotDelta - sqDisc) * invDenom,
- (-vDotDelta + sqDisc) * invDenom);
- } else {
- return {};
- }
- }
-
- TimeType getTimeOfIntersection(Line const& vLine, Plane const& vPlane) {
-
- auto const delta = vPlane.getCenter() - vLine.getStartPoint();
- auto const v = vLine.getVelocity();
- auto const n = vPlane.getNormal();
- auto const c = n.dot(v);
-
- if (c.magnitude() == 0) {
- return std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
- } else {
- return n.dot(delta) / c;
- }
- }
-
-} // namespace corsika::tracking_line
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
new file mode 100644
index 000000000..33ad57ede
--- /dev/null
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -0,0 +1,100 @@
+/*
+ * (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/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/geometry/StraightTrajectory.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
+#include <corsika/framework/logging/Logging.hpp>
+#include <corsika/modules/tracking/Intersect.hpp>
+
+#include <type_traits>
+#include <utility>
+
+namespace corsika::tracking_line {
+
+ template <typename TParticle>
+ inline auto Tracking::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] = nextIntersect(particle);
+
+ return std::make_tuple(StraightTrajectory(Line(initialPosition, initialVelocity),
+ minTime), // trajectory
+ minNode); // next volume node
+ }
+
+ template <typename TParticle, typename TMedium>
+ Intersections Tracking::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.getSquaredNorm();
+ 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();
+ }
+
+ template <typename TParticle, typename TBaseNodeType>
+ Intersections Tracking::intersect(TParticle const& particle,
+ TBaseNodeType const& volumeNode) {
+ Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
+ if (sphere) {
+ typedef typename std::remove_const_t<
+ std::remove_reference_t<decltype(volumeNode.getModelProperties())>>
+ medium_type;
+ return Tracking::intersect<TParticle, medium_type>(particle, *sphere,
+ volumeNode.getModelProperties());
+ }
+ throw std::runtime_error(
+ "The Volume type provided is not supported in Intersect(particle, node)");
+ }
+
+ template <typename TParticle, typename TMedium>
+ Intersections Tracking::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);
+ }
+
+} // namespace corsika::tracking_line
diff --git a/corsika/detail/stack/NuclearStackExtension.inl b/corsika/detail/stack/NuclearStackExtension.inl
index 2d10e28ee..b1fe9ba1c 100644
--- a/corsika/detail/stack/NuclearStackExtension.inl
+++ b/corsika/detail/stack/NuclearStackExtension.inl
@@ -13,7 +13,7 @@
#include <corsika/framework/stack/Stack.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/stack/SuperStupidStack.hpp>
+#include <corsika/stack/SimpleStack.hpp>
#include <algorithm>
#include <tuple>
diff --git a/corsika/detail/stack/SuperStupidStack.inl b/corsika/detail/stack/SuperStupidStack.inl
deleted file mode 100644
index 35d00a6d1..000000000
--- a/corsika/detail/stack/SuperStupidStack.inl
+++ /dev/null
@@ -1,93 +0,0 @@
-/*
- * (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/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/stack/Stack.hpp>
-
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-
-#include <string>
-#include <tuple>
-#include <vector>
-
-namespace corsika {
-
- template <typename StackIteratorInterface>
- void ParticleInterface<StackIteratorInterface>::setParticleData(
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v) {
- this->setPID(std::get<0>(v));
- this->setEnergy(std::get<1>(v));
- this->setMomentum(std::get<2>(v));
- this->setPosition(std::get<3>(v));
- this->setTime(std::get<4>(v));
- }
-
- template <typename StackIteratorInterface>
- void ParticleInterface<StackIteratorInterface>::setParticleData(
- ParticleInterface<StackIteratorInterface> const&,
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v) {
- this->setPID(std::get<0>(v));
- this->setEnergy(std::get<1>(v));
- this->setMomentum(std::get<2>(v));
- this->setPosition(std::get<3>(v));
- this->setTime(std::get<4>(v));
- }
-
- inline void SuperStupidStackImpl::clear() {
- dataPID_.clear();
- dataE_.clear();
- momentum_.clear();
- position_.clear();
- time_.clear();
- }
-
- inline void SuperStupidStackImpl::copy(size_t i1, size_t i2) {
- dataPID_[i2] = dataPID_[i1];
- dataE_[i2] = dataE_[i1];
- momentum_[i2] = momentum_[i1];
- position_[i2] = position_[i1];
- time_[i2] = time_[i1];
- }
-
- inline void SuperStupidStackImpl::swap(size_t i1, size_t i2) {
- std::swap(dataPID_[i2], dataPID_[i1]);
- std::swap(dataE_[i2], dataE_[i1]);
- std::swap(momentum_[i2], momentum_[i1]);
- std::swap(position_[i2], position_[i1]);
- std::swap(time_[i2], time_[i1]);
- }
-
- inline void SuperStupidStackImpl::incrementSize() {
- dataPID_.push_back(Code::Unknown);
- dataE_.push_back(0 * electronvolt);
-
- CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
-
- momentum_.push_back(
- MomentumVector(dummyCS, {0 * electronvolt, 0 * electronvolt, 0 * electronvolt}));
-
- position_.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
- time_.push_back(0 * second);
- }
-
- inline void SuperStupidStackImpl::decrementSize() {
- if (dataE_.size() > 0) {
- dataPID_.pop_back();
- dataE_.pop_back();
- momentum_.pop_back();
- position_.pop_back();
- time_.pop_back();
- }
- }
-
-} // namespace corsika
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index bec9a2700..2deb03a86 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -10,6 +10,11 @@
@file ParticleProperties.hpp
Interface to particle properties
+
+ The properties of all particles are saved in static and flat
+ arrays. There is a enum corsika::Code to identify each
+ particles, and each individual particles has its own static class,
+ which can be used to retrieve its physical properties.
*/
#pragma once
diff --git a/corsika/framework/geometry/IVolume.hpp b/corsika/framework/geometry/IVolume.hpp
index 59b9acb55..1960dfdf2 100644
--- a/corsika/framework/geometry/IVolume.hpp
+++ b/corsika/framework/geometry/IVolume.hpp
@@ -16,7 +16,7 @@ namespace corsika {
public:
//! returns true if the Point p is within the volume
- virtual bool isInside(Point const& p) const = 0;
+ virtual bool contains(Point const& p) const = 0;
virtual ~IVolume() = default;
};
diff --git a/corsika/framework/geometry/Sphere.hpp b/corsika/framework/geometry/Sphere.hpp
index cea283b3d..81e5e0b15 100644
--- a/corsika/framework/geometry/Sphere.hpp
+++ b/corsika/framework/geometry/Sphere.hpp
@@ -28,7 +28,7 @@ namespace corsika {
, radius_(pRadius) {}
//! returns true if the Point p is within the sphere
- bool isInside(Point const& p) const override;
+ bool contains(Point const& p) const override;
Point const& getCenter() const;
diff --git a/corsika/framework/geometry/Trajectory.hpp b/corsika/framework/geometry/Trajectory.hpp
deleted file mode 100644
index 6a3700eea..000000000
--- a/corsika/framework/geometry/Trajectory.hpp
+++ /dev/null
@@ -1,216 +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.
- */
-
-#pragma once
-
-#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;
-
- Trajectory(TType const& theT, TimeType timeLength)
- : TType(theT)
- , timeLength_(timeLength) {}
-
- Point getPosition(double const u) const;
-
- TimeType getDuration() const;
-
- LengthType getLength() const;
-
- LengthType getDistance(TimeType const t) const;
-
- void getLimitEndTo(LengthType const limit);
-
- auto getNormalizedDirection() const;
-
- private:
- TimeType timeLength_;
- };
- */
-
-} // namespace corsika
-
-#include <corsika/detail/framework/geometry/Trajectory.inl>
diff --git a/corsika/framework/stack/Stack.hpp b/corsika/framework/stack/Stack.hpp
index dab5bb2d9..4a6f09056 100644
--- a/corsika/framework/stack/Stack.hpp
+++ b/corsika/framework/stack/Stack.hpp
@@ -17,6 +17,30 @@
#include <utility>
#include <type_traits>
+/**
+ @file Stack.hpp
+
+ Description of particle stacks
+
+ In the CORSIKA 8 framework particle data is always stored in
+ particle stacks. A particle is, thus, always a reference (iterator)
+ to an entry on a stack, e.g.
+
+ \verbatim
+ ModelStack stack;
+ stack.begin(); // returns an iterator: StackIterator, ConstStackIterator
+
+ *stack.begin(); // return a reference to ParticleInterfaceType, which is the class
+provided by the user to read/write particle properties
+
+ \endverbatim
+
+ All functionality and algorithms for stack data access is located in the namespace
+corsika::stack
+
+ The minimal example of how to use this is shown in stack_example.cc
+**/
+
namespace corsika {
/**
diff --git a/corsika/media/Universe.hpp b/corsika/media/Universe.hpp
index af1ead8cb..cdd5c2bd2 100644
--- a/corsika/media/Universe.hpp
+++ b/corsika/media/Universe.hpp
@@ -15,7 +15,7 @@ namespace corsika {
struct Universe : public corsika::Sphere {
inline Universe(corsika::CoordinateSystemPtr const& pCS);
- inline bool isInside(corsika::Point const&) const override;
+ inline bool contains(corsika::Point const&) const override;
};
} // namespace corsika
diff --git a/corsika/media/VolumeTreeNode.hpp b/corsika/media/VolumeTreeNode.hpp
index 2647b4def..fba07bc58 100644
--- a/corsika/media/VolumeTreeNode.hpp
+++ b/corsika/media/VolumeTreeNode.hpp
@@ -29,9 +29,9 @@ namespace corsika {
: geoVolume_(std::move(pVolume)) {}
//! convenience function equivalent to Volume::isInside
- inline bool isInside(Point const& p) const;
+ inline bool contains(Point const& p) const;
- inline VolumeTreeNode<IModelProperties> const* isExcluded(Point const& p) const;
+ inline VolumeTreeNode<IModelProperties> const* excludes(Point const& p) const;
/** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
* \class Point \p p, or nullptr iff \p p is not contained in this volume.
diff --git a/corsika/modules/tracking/Intersect.hpp b/corsika/modules/tracking/Intersect.hpp
index 14ba805b7..459cef850 100644
--- a/corsika/modules/tracking/Intersect.hpp
+++ b/corsika/modules/tracking/Intersect.hpp
@@ -37,100 +37,18 @@ namespace corsika {
class Intersect {
protected:
+ /**
+ * Determines next intersection with any of the geometry volumes.
+ *
+ *
+ */
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);
- }
+ auto nextIntersect(TParticle const& particle,
+ TimeType const step_limit =
+ std::numeric_limits<TimeType::value_type>::infinity() *
+ second) const;
};
} // namespace corsika
+
+#include <corsika/detail/modules/tracking/Intersect.inl>
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
index bb75daf62..883f96dbc 100644
--- a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -12,7 +12,7 @@
#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/LeapFrogTrajectory.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/geometry/Intersections.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
@@ -29,6 +29,7 @@ namespace corsika {
namespace tracking_leapfrog_curved {
/**
+ * \file TrackingLeapFrogCurved.hpp
* \function LeapFrogStep
*
* Performs one leap-frog step consistent of two halve-steps with steplength/2
@@ -36,32 +37,7 @@ namespace corsika {
*boundary.
**/
template <typename TParticle>
- auto LeapFrogStep(TParticle const& particle, LengthType steplength) {
- if (particle.getMomentum().norm() == 0_GeV) {
- return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
- double(0));
- } // charge of the particle
- int const chargeNumber = particle.getChargeNumber();
- auto const* currentLogicalVolumeNode = particle.getNode();
- MagneticFieldVector const& magneticfield =
- currentLogicalVolumeNode->getModelProperties().getMagneticField(
- particle.getPosition());
- VelocityVector velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- decltype(meter / (second * volt)) k =
- chargeNumber * constants::cSquared * 1_eV /
- (velocity.getNorm() * particle.getEnergy() * 1_V);
- DirectionVector direction = velocity.normalized();
- auto position = particle.getPosition(); // First Movement
- // assuming magnetic field does not change during movement
- position =
- position + direction * steplength / 2; // Change of direction by magnetic field
- direction =
- direction + direction.cross(magneticfield) * steplength * k; // Second Movement
- position = position + direction * steplength / 2;
- auto steplength_true = steplength * (1.0 + (double)direction.getNorm()) / 2;
- return std::make_tuple(position, direction.normalized(), steplength_true);
- }
+ auto make_LeapFrogStep(TParticle const& particle, LengthType steplength);
/**
*
@@ -82,185 +58,15 @@ namespace corsika {
: straightTracking_{tracking_line::Tracking()} {}
template <typename TParticle>
- auto getTrack(TParticle const& particle) {
- VelocityVector const initialVelocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
-
- auto const position = particle.getPosition();
- CORSIKA_LOG_DEBUG(
- "Tracking pid: {}"
- " , E = {} GeV",
- particle.getPID(), particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("Tracking pos: {}", position.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());
-
- typedef
- typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
- node_type& volumeNode = *particle.getNode();
-
- // for the event of magnetic fields and curved trajectories, we need to limit
- // maximum step-length since we need to follow curved
- // trajectories segment-wise -- at least if we don't employ concepts as "Helix
- // Trajectories" or similar
- MagneticFieldVector const& magneticfield =
- volumeNode.getModelProperties().getMagneticField(position);
- MagneticFluxType const magnitudeB = magneticfield.getNorm();
- int const chargeNumber = particle.getChargeNumber();
- bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
-
- if (no_deflection) { return getLinearTrajectory(particle); }
-
- HEPMomentumType const pAlongB_delta =
- (particle.getMomentum() -
- particle.getMomentum().getParallelProjectionOnto(magneticfield))
- .getNorm();
-
- if (pAlongB_delta == 0_GeV) {
- // particle travel along, parallel to magnetic field. Rg is
- // "0", but for purpose of step limit we return infinity here.
- CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
- return getLinearTrajectory(particle);
- }
-
- LengthType const gyroradius =
- (pAlongB_delta * 1_V /
- (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
-
- const double maxRadians = 0.01;
- const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
- const TimeType steplimit_time = steplimit / initialVelocity.getNorm();
- CORSIKA_LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
- steplimit_time);
-
- // traverse the environment volume tree and find next
- // intersection
- auto [minTime, minNode] = nextIntersect(particle, steplimit_time);
-
- const auto k =
- chargeNumber * constants::cSquared * 1_eV / (particle.getEnergy() * 1_V);
- return std::make_tuple(
- LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
- minTime), // trajectory
- minNode); // next volume node
- }
+ auto getTrack(TParticle const& particle);
template <typename TParticle, typename TMedium>
static Intersections intersect(const TParticle& particle, const Sphere& sphere,
- const TMedium& medium) {
-
- if (sphere.getRadius() == 1_km * std::numeric_limits<double>::infinity()) {
- return Intersections();
- }
-
- int const chargeNumber = particle.getChargeNumber();
- auto const& position = particle.getPosition();
- MagneticFieldVector const& magneticfield = medium.getMagneticField(position);
-
- VelocityVector const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- DirectionVector const directionBefore =
- velocity.normalized(); // determine steplength to next volume
-
- auto const projectedDirection = directionBefore.cross(magneticfield);
- auto const projectedDirectionSqrNorm = projectedDirection.getSquaredNorm();
- bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
-
- if (chargeNumber == 0 || magneticfield.getNorm() == 0_T || isParallel) {
- return tracking_line::Tracking::intersect(particle, sphere, medium);
- }
-
- bool const numericallyInside = sphere.isInside(particle.getPosition());
-
- const auto absVelocity = velocity.getNorm();
- auto energy = particle.getEnergy();
- auto k = chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V);
-
- auto const denom =
- (directionBefore.cross(magneticfield)).getSquaredNorm() * k * k;
- const double a =
- ((directionBefore.cross(magneticfield)).dot(position - sphere.getCenter()) *
- k +
- 1) *
- 4 / (1_m * 1_m * denom);
- const double b = directionBefore.dot(position - sphere.getCenter()) * 8 /
- (denom * 1_m * 1_m * 1_m);
- const double c = ((position - sphere.getCenter()).getSquaredNorm() -
- (sphere.getRadius() * sphere.getRadius())) *
- 4 / (denom * 1_m * 1_m * 1_m * 1_m);
- CORSIKA_LOG_TRACE("denom={}, a={}, b={}, c={}", denom, a, b, c);
- std::complex<double>* solutions = quartic_solver::solve_quartic(0, a, b, c);
- LengthType d_enter, d_exit;
- int first = 0, first_entry = 0, first_exit = 0;
- for (int i = 0; i < 4; i++) {
- if (solutions[i].imag() == 0) {
- LengthType const dist = solutions[i].real() * 1_m;
- CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
- if (numericallyInside) {
- // there must be an entry (negative) and exit (positive) solution
- if (dist < -0.0001_m) { // security margin to assure transfer to next
- // logical volume
- if (first_entry == 0) {
- d_enter = dist;
- } else {
- d_enter = std::max(d_enter, dist); // closest negative to zero (-1e-4) m
- }
- first_entry++;
-
- } else { // thus, dist >= -0.0001_m
-
- if (first_exit == 0) {
- d_exit = dist;
- } else {
- d_exit = std::min(d_exit, dist); // closest positive to zero (-1e-4) m
- }
- first_exit++;
- }
- first = int(first_exit > 0) + int(first_entry > 0);
-
- } else { // thus, numericallyInside == false
-
- // both physical solutions (entry, exit) must be positive, and as small as
- // possible
- if (dist < -0.0001_m) { // need small numerical margin, to assure transport
- // into next logical volume
- continue;
- }
- if (first == 0) {
- d_enter = dist;
- } else {
- if (dist < d_enter) {
- d_exit = d_enter;
- d_enter = dist;
- } else {
- d_exit = dist;
- }
- }
- first++;
- }
- } // loop over solutions
- }
- delete[] solutions;
-
- if (first != 2) { // entry and exit points found
- CORSIKA_LOG_DEBUG("no intersection! count={}", first);
- return Intersections();
- }
- return Intersections(d_enter / absVelocity, d_exit / absVelocity);
- }
+ const TMedium& medium);
template <typename TParticle, typename TBaseNodeType>
static Intersections intersect(const TParticle& particle,
- const TBaseNodeType& volumeNode) {
- const Sphere* sphere = dynamic_cast<const Sphere*>(&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)");
- }
+ const TBaseNodeType& volumeNode);
protected:
/**
@@ -270,22 +76,7 @@ namespace corsika {
*
*/
template <typename TParticle>
- auto getLinearTrajectory(TParticle& particle) {
-
- // perform simple linear tracking
- auto [straightTrajectory, minNode] = straightTracking_.getTrack(particle);
-
- // return as leap-frog trajectory
- return std::make_tuple(
- LeapFrogTrajectory(
- straightTrajectory.getLine().getStartPoint(),
- straightTrajectory.getLine().getVelocity(),
- MagneticFieldVector(particle.getPosition().getCoordinateSystem(), 0_T,
- 0_T, 0_T),
- square(0_m) / (square(1_s) * 1_V),
- straightTrajectory.getDuration()), // trajectory
- minNode); // next volume node
- }
+ auto getLinearTrajectory(TParticle& particle);
protected:
tracking_line::Tracking
@@ -296,3 +87,5 @@ namespace corsika {
} // namespace tracking_leapfrog_curved
} // namespace corsika
+
+#include <corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl>
diff --git a/corsika/modules/tracking/TrackingLeapFrogStraight.hpp b/corsika/modules/tracking/TrackingLeapFrogStraight.hpp
index 0bc38b93d..f54e355ba 100644
--- a/corsika/modules/tracking/TrackingLeapFrogStraight.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogStraight.hpp
@@ -11,7 +11,7 @@
#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/StraightTrajectory.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
@@ -50,171 +50,16 @@ namespace corsika {
* relative to full linear step to next volume boundary. This
* should not be less than 0.5, otherwise you risk that
* particles will never travel from one volume to the next
- * one. A cross should be possible (even likely). If
- * firstFraction is too big (~1) the resulting calculated error
- * will be largest.
+ * one. A crossing must be possible (even likely). However, if
+ * firstFraction is too big (~1) the resulting calculated
+ * numerical error will be largest.
*
*/
Tracking(double const firstFraction = 0.55)
: firstFraction_(firstFraction) {}
template <typename Particle>
- auto getTrack(Particle& particle) {
- VelocityVector initialVelocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
-
- const Point initialPosition = particle.getPosition();
- CORSIKA_LOG_DEBUG(
- "TrackingB pid: {}"
- " , E = {} GeV",
- particle.getPID(), particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB pos: {}", initialPosition.getCoordinates());
- CORSIKA_LOG_DEBUG("TrackingB E: {} GeV", particle.getEnergy() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB p: {} GeV",
- particle.getMomentum().getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("TrackingB v: {} ", initialVelocity.getComponents());
-
- typedef decltype(particle.getNode()) node_type;
- node_type const volumeNode = particle.getNode();
-
- // check if particle is moving at all
- auto const absVelocity = initialVelocity.getNorm();
- if (absVelocity * 1_s == 0_m) {
- return std::make_tuple(
- LineTrajectory(Line(initialPosition, initialVelocity), 0_s), volumeNode);
- }
-
- // charge of the particle, and magnetic field
- const int chargeNumber = particle.getChargeNumber();
- auto magneticfield =
- volumeNode->getModelProperties().getMagneticField(initialPosition);
- const auto magnitudeB = magneticfield.getNorm();
- CORSIKA_LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
- magneticfield.getComponents() / 1_uT, chargeNumber,
- magnitudeB / 1_T);
- bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
-
- // check, where the first halve-step direction has geometric intersections
- const auto [initialTrack, initialTrackNextVolume] =
- tracking_line::Tracking::getTrack(particle);
- { [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
- const auto initialTrackLength = initialTrack.getLength(1);
-
- CORSIKA_LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
- initialTrack.getPosition(0).getCoordinates(),
- initialTrack.getPosition(1).getCoordinates(),
- initialTrackLength);
-
- // if particle is non-deflectable, we are done:
- if (no_deflection) {
- CORSIKA_LOG_DEBUG("no deflection. tracking finished");
- return std::make_tuple(initialTrack, initialTrackNextVolume);
- }
-
- HEPMomentumType const pAlongB_delta =
- (particle.getMomentum() -
- particle.getMomentum().getParallelProjectionOnto(magneticfield))
- .getNorm();
-
- if (pAlongB_delta == 0_GeV) {
- // particle travel along, parallel to magnetic field. Rg is
- // "0", but for purpose of step limit we return infinity here.
- CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
- return std::make_tuple(initialTrack, initialTrackNextVolume);
- }
-
- LengthType const gyroradius =
- (pAlongB_delta * 1_V /
- (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
-
- // we need to limit maximum step-length since we absolutely
- // need to follow strongly curved trajectories segment-wise,
- // at least if we don't employ concepts as "Helix
- // Trajectories" or similar
- const double maxRadians = 0.01;
- const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
- CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
-
- // calculate first halve step for "steplimit"
- auto const initialMomentum = particle.getMomentum();
- auto const absMomentum = initialMomentum.getNorm();
- DirectionVector const direction = initialVelocity.normalized();
-
- // avoid any intersections within first halve steplength
- LengthType const firstHalveSteplength =
- std::min(steplimit, initialTrackLength * firstFraction_);
-
- CORSIKA_LOG_DEBUG(
- "first halve step length {}, steplimit={}, initialTrackLength={}",
- firstHalveSteplength, steplimit, initialTrackLength);
- // perform the first halve-step
- const Point position_mid = initialPosition + direction * firstHalveSteplength;
- const auto k =
- chargeNumber * constants::c * 1_eV / (particle.getMomentum().getNorm() * 1_V);
- const auto new_direction =
- direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
- const auto new_direction_norm = new_direction.getNorm(); // by design this is >1
- CORSIKA_LOG_DEBUG(
- "position_mid={}, new_direction={}, (new_direction_norm)={}, deflection={}",
- position_mid.getCoordinates(), new_direction.getComponents(),
- new_direction_norm,
- acos(std::min(1.0, direction.dot(new_direction) / new_direction_norm)) * 180 /
- M_PI);
-
- // check, where the second halve-step direction has geometric intersections
- particle.setPosition(position_mid);
- particle.setMomentum(new_direction * absMomentum);
- const auto [finalTrack, finalTrackNextVolume] =
- tracking_line::Tracking::getTrack(particle);
- particle.setPosition(initialPosition); // this is not nice...
- particle.setMomentum(initialMomentum); // this is not nice...
-
- LengthType const finalTrackLength = finalTrack.getLength(1);
- LengthType const secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
-
- // check if volume transition is obvious, OR
- // for numerical reasons, particles slighly bend "away" from a
- // volume boundary have a very hard time to cross the border,
- // thus, if secondLeapFrogLength is just slighly shorter (1e-4m) than
- // finalTrackLength we better just [extend the
- // secondLeapFrogLength slightly and] force the volume
- // crossing:
- bool const switch_volume = finalTrackLength - 0.0001_m <= secondLeapFrogLength;
- LengthType const secondHalveStepLength =
- std::min(secondLeapFrogLength, finalTrackLength);
-
- CORSIKA_LOG_DEBUG(
- "finalTrack(0)={}, finalTrack(1)={}, finalTrackLength={}, "
- "secondLeapFrogLength={}, secondHalveStepLength={}, "
- "secondLeapFrogLength-finalTrackLength={}, "
- "secondHalveStepLength-finalTrackLength={}, "
- "nextVol={}, transition={}",
- finalTrack.getPosition(0).getCoordinates(),
- finalTrack.getPosition(1).getCoordinates(), finalTrackLength,
- secondLeapFrogLength, secondHalveStepLength,
- secondLeapFrogLength - finalTrackLength,
- secondHalveStepLength - finalTrackLength, fmt::ptr(finalTrackNextVolume),
- switch_volume);
-
- // perform the second halve-step
- auto const new_direction_normalized = new_direction.normalized();
- const Point finalPosition =
- position_mid + new_direction_normalized * secondHalveStepLength;
-
- const LengthType totalStep = firstHalveSteplength + secondHalveStepLength;
- const auto delta_pos = finalPosition - initialPosition;
- const auto distance = delta_pos.getNorm();
-
- return std::make_tuple(
- LineTrajectory(Line(initialPosition,
- (distance == 0_m ? initialVelocity
- : delta_pos.normalized() * absVelocity)),
- distance / absVelocity, // straight distance
- totalStep / absVelocity, // bend distance
- initialVelocity,
- new_direction_normalized * absVelocity), // trajectory
- (switch_volume ? finalTrackNextVolume : volumeNode));
- }
+ auto getTrack(Particle& particle);
protected:
double firstFraction_;
@@ -223,3 +68,5 @@ namespace corsika {
} // namespace tracking_leapfrog_straight
} // namespace corsika
+
+#include <corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl>
diff --git a/corsika/modules/tracking/TrackingStraight.hpp b/corsika/modules/tracking/TrackingStraight.hpp
index f563766d5..c595a75a5 100644
--- a/corsika/modules/tracking/TrackingStraight.hpp
+++ b/corsika/modules/tracking/TrackingStraight.hpp
@@ -13,7 +13,7 @@
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/framework/geometry/StraightTrajectory.hpp>
#include <corsika/framework/geometry/Intersections.hpp>
#include <corsika/framework/logging/Logging.hpp>
#include <corsika/modules/tracking/Intersect.hpp>
@@ -24,9 +24,10 @@
namespace corsika::tracking_line {
/**
- * \class Tracking
- *
+ * Tracking of particles without charge or in no magnetic fields.
*
+ * This will simply follow traight lines, but intersect them with
+ * the geometry volume model of the environment.
*
**/
@@ -35,77 +36,26 @@ namespace corsika::tracking_line {
using Intersect<Tracking>::nextIntersect;
public:
+ //! Determine track of particle
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] = nextIntersect(particle);
-
- return std::make_tuple(LineTrajectory(Line(initialPosition, initialVelocity),
- minTime), // trajectory
- minNode); // next volume node
- }
+ auto getTrack(TParticle const& particle);
+ //! find intersection of Sphere with Track
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.getSquaredNorm();
- 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();
- }
+ TMedium const&);
+ //! find intersection of Volume with Track
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)");
- }
+ TBaseNodeType const& volumeNode);
+ //! find intersection of Plane with Track
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);
- }
+ TMedium const&);
};
} // namespace corsika::tracking_line
-// #include <corsika/detail/modules/tracking/TrackingLine.inl>
+#include <corsika/detail/modules/tracking/TrackingStraight.inl>
diff --git a/corsika/setup/SetupTrajectory.hpp b/corsika/setup/SetupTrajectory.hpp
index a20d1db3d..fdbfcbeaa 100644
--- a/corsika/setup/SetupTrajectory.hpp
+++ b/corsika/setup/SetupTrajectory.hpp
@@ -8,7 +8,8 @@
#pragma once
-#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/framework/geometry/StraightTrajectory.hpp>
+#include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
#include <corsika/modules/Tracking.hpp>
namespace corsika::setup {
@@ -45,7 +46,7 @@ namespace corsika::setup {
The default trajectory.
*/
/// definition of Trajectory base class, to be used in tracking and cascades
- typedef LineTrajectory Trajectory;
+ typedef StraightTrajectory Trajectory;
// typedef corsika::geometry::LeapFrogTrajectory Trajectory;
} // namespace corsika::setup
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index c9559a753..bcd187c05 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -14,7 +14,7 @@
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
-#include <corsika/stack/SuperStupidStack.hpp>
+#include <corsika/stack/SimpleStack.hpp>
#include <algorithm>
#include <tuple>
@@ -27,7 +27,7 @@ namespace corsika::nuclear_stack {
* Define ParticleInterface for NuclearStackExtension Stack derived from
* ParticleInterface of Inner stack class
*
- * Add A and Z data to existing stack (currently SuperStupidStack) of particle
+ * Add A and Z data to existing stack (currently SimpleStack) of particle
* properties. This is done via inheritance, not via CombinedStack since the nuclear
* data is stored ONLY when needed (for nuclei) and not for all particles. Thus, this is
* a new, derived Stack object.
@@ -210,11 +210,11 @@ namespace corsika::nuclear_stack {
//
template <typename TStackIter>
using ExtendedParticleInterfaceType =
- NuclearParticleInterface<SuperStupidStack::pi_type, TStackIter>;
+ NuclearParticleInterface<SimpleStack::pi_type, TStackIter>;
// the particle data stack with extra nuclear information:
using ParticleDataStack =
- NuclearStackExtension<SuperStupidStack, ExtendedParticleInterfaceType>;
+ NuclearStackExtension<SimpleStack, ExtendedParticleInterfaceType>;
} // namespace corsika::nuclear_stack
diff --git a/corsika/stack/SuperStupidStack.hpp b/corsika/stack/SuperStupidStack.hpp
deleted file mode 100644
index 81e44a49c..000000000
--- a/corsika/stack/SuperStupidStack.hpp
+++ /dev/null
@@ -1,176 +0,0 @@
-/*
- * (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/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/stack/Stack.hpp>
-
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
-
-#include <string>
-#include <tuple>
-#include <vector>
-
-namespace corsika {
-
- /**
- * Example of a particle object on the stack.
- */
-
- template <typename StackIteratorInterface>
- struct ParticleInterface : public ParticleBase<StackIteratorInterface> {
-
- private:
- typedef ParticleBase<StackIteratorInterface> super_type;
-
- public:
- std::string asString() const {
- return fmt::format("particle: i={}, PID={}, E={}GeV", super_type::getIndex(),
- get_name(this->getPID()), this->getEnergy() / 1_GeV);
- }
-
- void setParticleData(
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v);
-
- void setParticleData(
- ParticleInterface<StackIteratorInterface> const&,
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v);
-
- /// individual setters
- void setPID(Code const id) {
- super_type::getStackData().setPID(super_type::getIndex(), id);
- }
- void setEnergy(HEPEnergyType const& e) {
- super_type::getStackData().setEnergy(super_type::getIndex(), e);
- }
- void setMomentum(MomentumVector const& v) {
- super_type::getStackData().setMomentum(super_type::getIndex(), v);
- }
- void setPosition(Point const& v) {
- super_type::getStackData().setPosition(super_type::getIndex(), v);
- }
- void setTime(TimeType const& v) {
- super_type::getStackData().setTime(super_type::getIndex(), v);
- }
-
- /// individual getters
- Code getPID() const {
- return super_type::getStackData().getPID(super_type::getIndex());
- }
- HEPEnergyType getEnergy() const {
- return super_type::getStackData().getEnergy(super_type::getIndex());
- }
- MomentumVector getMomentum() const {
- return super_type::getStackData().getMomentum(super_type::getIndex());
- }
- Point getPosition() const {
- return super_type::getStackData().getPosition(super_type::getIndex());
- }
- TimeType getTime() const {
- return super_type::getStackData().getTime(super_type::getIndex());
- }
- /**
- * @name derived quantities
- *
- * @{
- */
- DirectionVector getDirection() const {
- return this->getMomentum() / this->getEnergy();
- }
-
- HEPMassType getMass() const { return get_mass(this->getPID()); }
-
- int16_t getChargeNumber() const { return get_charge_number(this->getPID()); }
- ///@}
- };
-
- /**
- * Memory implementation of the most simple (stupid) particle stack object.
- *
- */
-
- class SuperStupidStackImpl {
-
- public:
- typedef std::vector<Code> code_vector_type;
- typedef std::vector<HEPEnergyType> energy_vector_type;
- typedef std::vector<Point> point_vector_type;
- typedef std::vector<TimeType> time_vector_type;
- typedef std::vector<MomentumVector> momentum_vector_type;
-
- SuperStupidStackImpl() = default;
-
- SuperStupidStackImpl(SuperStupidStackImpl const& other) = default;
-
- SuperStupidStackImpl(SuperStupidStackImpl&& other) = default;
-
- SuperStupidStackImpl& operator=(SuperStupidStackImpl const& other) = default;
-
- SuperStupidStackImpl& operator=(SuperStupidStackImpl&& other) = default;
-
- void dump() const {}
-
- inline void clear();
-
- unsigned int getSize() const { return dataPID_.size(); }
- unsigned int getCapacity() const { return dataPID_.size(); }
-
- void setPID(size_t i, Code const id) { dataPID_[i] = id; }
- void setEnergy(size_t i, HEPEnergyType const& e) { dataE_[i] = e; }
- void setMomentum(size_t i, MomentumVector const& v) { momentum_[i] = v; }
- void setPosition(size_t i, Point const& v) { position_[i] = v; }
- void setTime(size_t i, TimeType const& v) { time_[i] = v; }
-
- Code getPID(size_t i) const { return dataPID_[i]; }
-
- HEPEnergyType getEnergy(size_t i) const { return dataE_[i]; }
-
- MomentumVector getMomentum(size_t i) const { return momentum_[i]; }
-
- Point getPosition(size_t i) const { return position_[i]; }
- TimeType getTime(size_t i) const { return time_[i]; }
-
- HEPEnergyType getDataE(size_t i) const { return dataE_[i]; }
-
- void setDataE(size_t i, HEPEnergyType const& dataE) { dataE_[i] = dataE; }
-
- Code getDataPid(size_t i) const { return dataPID_[i]; }
-
- void setDataPid(size_t i, Code const& dataPid) { dataPID_[i] = dataPid; }
- /**
- * Function to copy particle at location i2 in stack to i1
- */
- inline void copy(size_t i1, size_t i2);
-
- /**
- * Function to copy particle at location i2 in stack to i1
- */
- inline void swap(size_t i1, size_t i2);
-
- inline void incrementSize();
- inline void decrementSize();
-
- private:
- /// the actual memory to store particle data
- code_vector_type dataPID_;
- energy_vector_type dataE_;
- momentum_vector_type momentum_;
- point_vector_type position_;
- time_vector_type time_;
-
- }; // end class SuperStupidStackImpl
-
- typedef Stack<SuperStupidStackImpl, ParticleInterface> SuperStupidStack;
-
-} // namespace corsika
-
-#include <corsika/detail/stack/SuperStupidStack.inl>
diff --git a/examples/geometry_example.cpp b/examples/geometry_example.cpp
index cca5656e3..8c471cb60 100644
--- a/examples/geometry_example.cpp
+++ b/examples/geometry_example.cpp
@@ -51,12 +51,12 @@ int main() {
assert(norm == 1 * meter * meter);
Sphere s(p1, 10_m); // define a sphere around a point with a radius
- CORSIKA_LOG_INFO("p1 inside s:{} ", s.isInside(p2));
- assert(s.isInside(p2) == 1);
+ CORSIKA_LOG_INFO("p1 inside s:{} ", s.contains(p2));
+ assert(s.contains(p2) == 1);
Sphere s2(p1, 3_um); // another sphere
- CORSIKA_LOG_INFO("p1 inside s2: {}", s2.isInside(p2));
- assert(s2.isInside(p2) == 0);
+ CORSIKA_LOG_INFO("p1 inside s2: {}", s2.contains(p2));
+ assert(s2.contains(p2) == 0);
// let's try parallel projections:
auto const v1 = Vector<length_d>(root, {1_m, 1_m, 0_m});
diff --git a/examples/stack_example.cpp b/examples/stack_example.cpp
index 5e5f3ba70..196efa5c9 100644
--- a/examples/stack_example.cpp
+++ b/examples/stack_example.cpp
@@ -7,7 +7,7 @@
*/
#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/stack/SuperStupidStack.hpp>
+#include <corsika/stack/SimpleStack.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
@@ -19,7 +19,7 @@
using namespace corsika;
using namespace std;
-void fill(SuperStupidStack& s) {
+void fill(SimpleStack& s) {
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
for (int i = 0; i < 11; ++i) {
s.addParticle(std::make_tuple(Code::Electron, 1.5_GeV * i,
@@ -28,7 +28,7 @@ void fill(SuperStupidStack& s) {
}
}
-void read(SuperStupidStack& s) {
+void read(SimpleStack& s) {
assert(s.getEntries() == 11); // stack has 11 particles
HEPEnergyType total_energy;
@@ -44,7 +44,7 @@ void read(SuperStupidStack& s) {
int main() {
std::cout << "stack_example" << std::endl;
- SuperStupidStack s;
+ SimpleStack s;
fill(s);
read(s);
return 0;
diff --git a/tests/common/SetupTestTrajectory.hpp b/tests/common/SetupTestTrajectory.hpp
index dd028b0c7..0c06d7f39 100644
--- a/tests/common/SetupTestTrajectory.hpp
+++ b/tests/common/SetupTestTrajectory.hpp
@@ -11,8 +11,9 @@
#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/framework/geometry/StraightTrajectory.hpp>
+// #include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
// #include <corsika/modules/TrackingLine.hpp>
// #include <corsika/modules/TrackingCurved.hpp> // simple leap-frog implementation
// #include <corsika/modules/TrackingLeapFrog.hpp> // more complete leap-frog
@@ -24,9 +25,9 @@ namespace corsika::setup::testing {
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);
+ inline StraightTrajectory make_track<StraightTrajectory>(Line const& line,
+ TimeType const tEnd) {
+ return StraightTrajectory(line, tEnd);
}
/*
diff --git a/tests/common/testTestStack.hpp b/tests/common/testTestStack.hpp
index a905231fb..70f4d5f14 100644
--- a/tests/common/testTestStack.hpp
+++ b/tests/common/testTestStack.hpp
@@ -36,9 +36,8 @@ public:
void setData(const unsigned int i, const double v) { data_[i] = v; }
double getData(const unsigned int i) const { return data_[i]; }
-
// these functions are also needed by the Stack interface
- void incrementSize() { data_.push_back(0.); }
+ void incrementSize() { data_.push_back(0.); }
void decrementSize() {
if (data_.size() > 0) { data_.pop_back(); }
}
@@ -60,13 +59,11 @@ private:
*
*/
template <typename StackIteratorInterface>
-class TestParticleInterface
- : public corsika::ParticleBase<StackIteratorInterface> {
-
+class TestParticleInterface : public corsika::ParticleBase<StackIteratorInterface> {
+
typedef corsika::ParticleBase<StackIteratorInterface> super_type;
-
-public:
+public:
/*
The SetParticleData methods are called for creating new entries
on the stack. You can specifiy various parametric versions to
@@ -81,6 +78,10 @@ public:
}
// here are the fundamental methods for access to TestStackData data
- void setData(const double v) { super_type::getStackData().setData(super_type::getIndex(), v); }
- double getData() const { return super_type::getStackData().getData(super_type::getIndex()); }
+ void setData(const double v) {
+ super_type::getStackData().setData(super_type::getIndex(), v);
+ }
+ double getData() const {
+ return super_type::getStackData().getData(super_type::getIndex());
+ }
};
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index 6664f0acb..27880dfb8 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -68,7 +68,7 @@ public:
VelocityVector const initialVelocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
return std::make_tuple(
- LineTrajectory(
+ StraightTrajectory(
Line(particle.getPosition(), initialVelocity),
std::numeric_limits<TimeType::value_type>::infinity() * 1_s), // trajectory,
// just
diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp
index 8f9ef0b85..3472cf2dd 100644
--- a/tests/framework/testGeometry.cpp
+++ b/tests/framework/testGeometry.cpp
@@ -16,7 +16,8 @@
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/framework/geometry/StraightTrajectory.hpp>
+#include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
#include <PhysicalUnitsCatch2.hpp> // namespace corsike::testing
@@ -228,7 +229,7 @@ TEST_CASE("CoordinateSystem hirarchy") {
}
TEST_CASE("Sphere") {
- CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
+ CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
Point center(rootCS, {0_m, 3_m, 4_m});
Sphere sphere(center, 5_m);
@@ -241,8 +242,8 @@ TEST_CASE("Sphere") {
}
SECTION("isInside") {
- CHECK_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m})));
- CHECK(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m})));
+ CHECK_FALSE(sphere.contains(Point(rootCS, {100_m, 0_m, 0_m})));
+ CHECK(sphere.contains(Point(rootCS, {2_m, 3_m, 4_m})));
}
}
@@ -270,7 +271,7 @@ TEST_CASE("Trajectories") {
.magnitude() == Approx(0).margin(absMargin));
auto const t = 1_s;
- LineTrajectory base(line, t);
+ StraightTrajectory base(line, t);
CHECK(line.getPosition(t).getCoordinates() == base.getPosition(1.).getCoordinates());
CHECK((base.getDirection(0).getComponents(rootCS) -
diff --git a/tests/framework/testRandom.cpp b/tests/framework/testRandom.cpp
index f31a2d9dc..17c9c7c7a 100644
--- a/tests/framework/testRandom.cpp
+++ b/tests/framework/testRandom.cpp
@@ -23,25 +23,25 @@ SCENARIO("random-number streams can be registered and retrieved") {
RNGManager& rngManager = RNGManager::getInstance();
WHEN("the sequence name is not registered") {
- REQUIRE(rngManager.isRegistered("stream_A") == false);
+ CHECK(rngManager.isRegistered("stream_A") == false);
THEN("a sequence is registered by name") {
rngManager.registerRandomStream("stream_A");
THEN("the sequence can be retrieved") {
- REQUIRE_NOTHROW(rngManager.getRandomStream("stream_A"));
+ CHECK_NOTHROW(rngManager.getRandomStream("stream_A"));
THEN("we can check that the sequence exists") {
- REQUIRE_NOTHROW(rngManager.getRandomStream("stream_A"));
+ CHECK_NOTHROW(rngManager.getRandomStream("stream_A"));
THEN("an unknown sequence cannot be retrieved") {
- REQUIRE(rngManager.isRegistered("stream_A") == true);
+ CHECK(rngManager.isRegistered("stream_A") == true);
THEN("an unknown sequence cannot be retrieved") {
- REQUIRE_THROWS(rngManager.getRandomStream("stream_UNKNOWN"));
+ CHECK_THROWS(rngManager.getRandomStream("stream_UNKNOWN"));
THEN("an unknown sequence is not registered") {
- REQUIRE(rngManager.isRegistered("stream_UNKNOWN") == false);
+ CHECK(rngManager.isRegistered("stream_UNKNOWN") == false);
}
}
}
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 01f18af02..d72de5b53 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -2,17 +2,15 @@ set (test_modules_sources
TestMain.cpp
testStackInspector.cpp
testTracking.cpp
- # testExecTime.cpp
+ # testExecTime.cpp --> needs to be fixed, see #326
testObservationPlane.cpp
testQGSJetII.cpp
- # testSwitchProcess.cpp -> gone
testPythia8.cpp
testUrQMD.cpp
testCONEX.cpp
testOnShellCheck.cpp
testParticleCut.cpp
testSibyll.cpp
- # testNullModel.cpp -> gone, framework now
)
CORSIKA_ADD_TEST (testModules SOURCES ${test_modules_sources})
diff --git a/tests/modules/testInteractionCounter.cpp b/tests/modules/testInteractionCounter.cpp
index 62a8c17df..6a147372f 100644
--- a/tests/modules/testInteractionCounter.cpp
+++ b/tests/modules/testInteractionCounter.cpp
@@ -1,5 +1,5 @@
/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2021 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
@@ -46,7 +46,7 @@ TEST_CASE("InteractionCounter", "[process]") {
InteractionCounter countedProcess(d);
SECTION("GetInteractionLength") {
- REQUIRE(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
+ CHECK(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm);
}
auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
@@ -56,19 +56,19 @@ TEST_CASE("InteractionCounter", "[process]") {
unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
Z, 105_TeV, nodePtr, *csPtr);
- REQUIRE(stackPtr->getEntries() == 1);
- REQUIRE(secViewPtr->getEntries() == 0);
+ CHECK(stackPtr->getEntries() == 1);
+ CHECK(secViewPtr->getEntries() == 0);
auto const ret = countedProcess.DoInteraction(*secViewPtr);
- REQUIRE(ret == nullptr);
+ CHECK(ret == nullptr);
auto const& h = countedProcess.GetHistogram().labHist();
- REQUIRE(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
- REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
+ CHECK(h.at(h.axis(0).index(1'000'070'140), h.axis(1).index(1.05e14)) == 1);
+ CHECK(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
auto const& h2 = countedProcess.GetHistogram().CMSHist();
- REQUIRE(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
- REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
+ CHECK(h2.at(h2.axis(0).index(1'000'070'140), h2.axis(1).index(1.6e12)) == 1);
+ CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
countedProcess.GetHistogram().saveLab("testInteractionCounter_file1.npz",
utl::SaveMode::overwrite);
@@ -80,19 +80,19 @@ TEST_CASE("InteractionCounter", "[process]") {
auto constexpr code = particles::Code::Lambda0;
auto [stackPtr, secViewPtr] =
setup::testing::setupStack(code, 0, 0, 105_TeV, nodePtr, *csPtr);
- REQUIRE(stackPtr->getEntries() == 1);
- REQUIRE(secViewPtr->getEntries() == 0);
+ CHECK(stackPtr->getEntries() == 1);
+ CHECK(secViewPtr->getEntries() == 0);
auto const ret = countedProcess.DoInteraction(*secViewPtr);
- REQUIRE(ret == nullptr);
+ CHECK(ret == nullptr);
auto const& h = countedProcess.GetHistogram().labHist();
- REQUIRE(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
- REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
+ CHECK(h.at(h.axis(0).index(3122), h.axis(1).index(1.05e14)) == 1);
+ CHECK(std::accumulate(h.cbegin(), h.cend(), 0) == 1);
auto const& h2 = countedProcess.GetHistogram().CMSHist();
- REQUIRE(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
- REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
+ CHECK(h2.at(h2.axis(0).index(3122), h2.axis(1).index(1.6e12)) == 1);
+ CHECK(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1);
}
}
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 00f0d8013..eeb6d95a7 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -42,9 +42,9 @@ TEST_CASE("CORSIKA_DATA", "[processes]") {
SECTION("check CORSIKA_DATA") {
const char* data = std::getenv("CORSIKA_DATA");
- // these REQUIRES are needed:
- REQUIRE(data != 0);
- REQUIRE(std::experimental::filesystem::is_directory(
+ // these CHECKS are needed:
+ CHECK(data != 0);
+ CHECK(std::experimental::filesystem::is_directory(
std::experimental::filesystem::path(std::string(data) + "/QGSJetII")));
std::cout << "data: " << data << " isDir: "
<< std::experimental::filesystem::is_directory(std::string(data) +
@@ -62,37 +62,37 @@ TEST_CASE("QgsjetII", "[processes]") {
}
SECTION("QgsjetII -> Corsika") {
- REQUIRE(Code::PiPlus == corsika::qgsjetII::convertFromQgsjetII(
- corsika::qgsjetII::QgsjetIICode::PiPlus));
+ CHECK(Code::PiPlus == corsika::qgsjetII::convertFromQgsjetII(
+ corsika::qgsjetII::QgsjetIICode::PiPlus));
}
SECTION("Corsika -> QgsjetII") {
- REQUIRE(corsika::qgsjetII::convertToQgsjetII(Code::PiMinus) ==
- corsika::qgsjetII::QgsjetIICode::PiMinus);
- REQUIRE(corsika::qgsjetII::convertToQgsjetIIRaw(Code::Proton) == 2);
+ CHECK(corsika::qgsjetII::convertToQgsjetII(Code::PiMinus) ==
+ corsika::qgsjetII::QgsjetIICode::PiMinus);
+ CHECK(corsika::qgsjetII::convertToQgsjetIIRaw(Code::Proton) == 2);
}
SECTION("canInteractInQgsjetII") {
- REQUIRE(corsika::qgsjetII::canInteract(Code::Proton));
- REQUIRE(corsika::qgsjetII::canInteract(Code::KPlus));
- REQUIRE(corsika::qgsjetII::canInteract(Code::Nucleus));
- // REQUIRE(corsika::qgsjetII::canInteract(Helium::getCode()));
+ CHECK(corsika::qgsjetII::canInteract(Code::Proton));
+ CHECK(corsika::qgsjetII::canInteract(Code::KPlus));
+ CHECK(corsika::qgsjetII::canInteract(Code::Nucleus));
+ // CHECK(corsika::qgsjetII::canInteract(Helium::getCode()));
- REQUIRE_FALSE(corsika::qgsjetII::canInteract(Code::EtaC));
- REQUIRE_FALSE(corsika::qgsjetII::canInteract(Code::SigmaC0));
+ CHECK_FALSE(corsika::qgsjetII::canInteract(Code::EtaC));
+ CHECK_FALSE(corsika::qgsjetII::canInteract(Code::SigmaC0));
}
SECTION("cross-section type") {
- REQUIRE(corsika::qgsjetII::getQgsjetIIXSCode(Code::Neutron) ==
- corsika::qgsjetII::QgsjetIIXSClass::Baryons);
- REQUIRE(corsika::qgsjetII::getQgsjetIIXSCode(Code::K0Long) ==
- corsika::qgsjetII::QgsjetIIXSClass::Kaons);
- REQUIRE(corsika::qgsjetII::getQgsjetIIXSCode(Code::Proton) ==
- corsika::qgsjetII::QgsjetIIXSClass::Baryons);
- REQUIRE(corsika::qgsjetII::getQgsjetIIXSCode(Code::PiMinus) ==
- corsika::qgsjetII::QgsjetIIXSClass::LightMesons);
+ CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::Neutron) ==
+ corsika::qgsjetII::QgsjetIIXSClass::Baryons);
+ CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::K0Long) ==
+ corsika::qgsjetII::QgsjetIIXSClass::Kaons);
+ CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::Proton) ==
+ corsika::qgsjetII::QgsjetIIXSClass::Baryons);
+ CHECK(corsika::qgsjetII::getQgsjetIIXSCode(Code::PiMinus) ==
+ corsika::qgsjetII::QgsjetIIXSClass::LightMesons);
}
}
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
index a28d274f3..31f3765ae 100644
--- a/tests/modules/testStackInspector.cpp
+++ b/tests/modules/testStackInspector.cpp
@@ -27,7 +27,7 @@ TEST_CASE("StackInspector", "[processes]") {
Point const origin(rootCS, {0_m, 0_m, 0_m});
VelocityVector v(rootCS, 0_m / second, 0_m / second, 1_m / second);
Line line(origin, v);
- LineTrajectory track(line, 10_s);
+ StraightTrajectory track(line, 10_s);
TestCascadeStack stack;
stack.clear();
diff --git a/tests/modules/testSwitchProcess.cpp b/tests/modules/testSwitchProcess.cpp
deleted file mode 100644
index 296315768..000000000
--- a/tests/modules/testSwitchProcess.cpp
+++ /dev/null
@@ -1,258 +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/modules/SwitchProcess.hpp>
-
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/stack/SecondaryView.hpp>
-#include <corsika/framework/stack/Stack.hpp>
-
-#include <catch2/catch.hpp>
-
-#include <algorithm>
-#include <random>
-
-using namespace corsika;
-using namespace corsika::units::si;
-
-class TestStackData {
-
-public:
- // these functions are needed for the Stack interface
- void Init() {}
- void Clear() { fData.clear(); }
- unsigned int GetSize() const { return fData.size(); }
- unsigned int GetCapacity() const { return fData.size(); }
- void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
- void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
-
- // custom data access function
- void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
- HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
-
- // these functions are also needed by the Stack interface
- void IncrementSize() { fData.resize(fData.size() + 1); }
- void DecrementSize() {
- if (fData.size() > 0) { fData.pop_back(); }
- }
-
- // custom private data section
-private:
- std::vector<HEPEnergyType> fData;
-};
-
-/**
- * From static_cast of a StackIterator over entries in the
- * TestStackData class you get and object of type
- * TestParticleInterface defined here
- *
- * It provides Get/Set methods to read and write data to the "Data"
- * storage of TestStackData obtained via
- * "StackIteratorInterface::GetStackData()", given the index of the
- * iterator "StackIteratorInterface::GetIndex()"
- *
- */
-template <typename StackIteratorInterface>
-class TestParticleInterface : public corsika::ParticleBase<StackIteratorInterface> {
-
-public:
- using corsika::ParticleBase<StackIteratorInterface>::GetStackData;
- using corsika::ParticleBase<StackIteratorInterface>::GetIndex;
-
- /*
- The SetParticleData methods are called for creating new entries
- on the stack. You can specifiy various parametric versions to
- perform this task:
- */
-
- // default version for particle-creation from input data
- void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
- void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
- std::tuple<HEPEnergyType> v) {
- SetEnergy(std::get<0>(v));
- }
-
- // here are the fundamental methods for access to TestStackData data
- void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
- HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
-};
-
-using SimpleStack = corsika::Stack<TestStackData, TestParticleInterface>;
-
-// see issue 161
-#if defined(__clang__)
-using StackTestView = corsika::SecondaryView<TestStackData, TestParticleInterface>;
-#elif defined(__GNUC__) || defined(__GNUG__)
-using StackTestView = corsika::MakeView<SimpleStack>::type;
-#endif
-
-auto constexpr kgMSq = 1_kg / (1_m * 1_m);
-
-template <int N>
-struct DummyProcess : InteractionProcess<DummyProcess<N>> {
- void Init() {}
-
- template <typename TParticle>
- corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
- return N * kgMSq;
- }
-
- template <typename TSecondaries>
- void doInteraction(TSecondaries&) {
- // to figure out which process was selected in the end, we produce N
- // secondaries for DummyProcess<N>
-
- for (int i = 0; i < N; ++i) {
- // vSec.AddSecondary(std::make_tuple(vSec.GetEnergy() / N)); // <-- FIXME, when
- // SwitchProcess is removedalter
- }
- }
-};
-
-using DummyLowEnergyProcess = DummyProcess<1>;
-using DummyHighEnergyProcess = DummyProcess<2>;
-using DummyAdditionalProcess = DummyProcess<3>;
-
-TEST_CASE("SwitchProcess from InteractionProcess") {
- DummyLowEnergyProcess low;
- DummyHighEnergyProcess high;
- DummyAdditionalProcess proc;
-
- switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
- auto seq = corsika::make_sequence(switchProcess, proc);
-
- SimpleStack stack;
-
- SECTION("low energy") {
- stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
- auto p = stack.GetNextParticle();
-
- // low energy process returns 1 kg/m²
- SECTION("interaction length") {
- REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
- REQUIRE(seq.getInteractionLength(p) / kgMSq == Approx(3. / 4));
- }
-
- // low energy process creates 1 secondary
- //~ SECTION("SelectInteraction") {
- //~ typename SimpleStack::ParticleType theParticle =
- //~ stack.GetNextParticle(); // as in corsika::Cascade
- //~ StackTestView view(theParticle);
- //~ auto projectile = view.GetProjectile();
-
- //~ InverseGrammageType invLambda = 0 / kgMSq;
- //~ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
-
- //~ REQUIRE(view.GetSize() == 1);
- //~ }
- }
-
- SECTION("high energy") {
- stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
- auto p = stack.GetNextParticle();
-
- // high energy process returns 2 kg/m²
- SECTION("interaction length") {
- REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
- REQUIRE(seq.getInteractionLength(p) / kgMSq == Approx(6. / 5));
- }
-
- // high energy process creates 2 secondaries
- SECTION("SelectInteraction") {
- typename SimpleStack::ParticleType theParticle =
- stack.GetNextParticle(); // as in corsika::Cascade
- StackTestView view(theParticle);
- auto projectile = view.GetProjectile();
-
- InverseGrammageType invLambda = 0 / kgMSq;
- switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
-
- REQUIRE(view.GetSize() == 2);
- }
- }
-}
-
-TEST_CASE("SwitchProcess from ProcessSequence") {
- DummyProcess<1> innerA;
- DummyProcess<2> innerB;
- DummyProcess<3> outer;
- DummyProcess<4> additional;
-
- auto seq = corsika::make_sequence(innerA, innerB);
-
- switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
- auto completeSeq = corsika::make_sequence(switchProcess, additional);
-
- SimpleStack stack;
-
- SECTION("low energy") {
- stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
- auto p = stack.GetNextParticle();
-
- SECTION("interaction length") {
- REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
- REQUIRE(completeSeq.getInteractionLength(p) / kgMSq == Approx(4. / 7));
- }
-
- SECTION("SelectInteraction") {
- std::vector<int> numberOfSecondaries;
-
- for (int i = 0; i < 1000; ++i) {
- typename SimpleStack::ParticleType theParticle =
- stack.GetNextParticle(); // as in corsika::Cascade
- StackTestView view(theParticle);
-
- double r = i / 1000.;
- InverseGrammageType invLambda = r * 7. / 4 / kgMSq;
-
- InverseGrammageType accumulator = 0 / kgMSq;
- completeSeq.selectInteraction(view, invLambda, accumulator);
-
- numberOfSecondaries.push_back(view.GetSize());
- }
-
- auto const mean =
- std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
- numberOfSecondaries.size();
- REQUIRE(mean == Approx(12. / 7.).margin(0.01));
- }
- }
-
- SECTION("high energy") {
- stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
- auto p = stack.GetNextParticle();
-
- SECTION("interaction length") {
- REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
- REQUIRE(completeSeq.getInteractionLength(p) / kgMSq == Approx(12. / 7.));
- }
-
- SECTION("SelectInteraction") {
- std::vector<int> numberOfSecondaries;
-
- for (int i = 0; i < 1000; ++i) {
- typename SimpleStack::ParticleType theParticle =
- stack.GetNextParticle(); // as in corsika::Cascade
- StackTestView view(theParticle);
-
- double r = i / 1000.;
- InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
-
- InverseGrammageType accumulator = 0 / kgMSq;
- completeSeq.selectInteraction(view, invLambda, accumulator);
-
- numberOfSecondaries.push_back(view.GetSize());
- }
-
- auto const mean =
- std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
- numberOfSecondaries.size();
- REQUIRE(mean == Approx(24. / 7.).margin(0.01));
- }
- }
-}
diff --git a/tests/stack/CMakeLists.txt b/tests/stack/CMakeLists.txt
index d559a5e13..ffb635624 100644
--- a/tests/stack/CMakeLists.txt
+++ b/tests/stack/CMakeLists.txt
@@ -4,7 +4,7 @@ set (test_stack_sources
testHistoryView.cpp
testGeometryNodeStackExtension.cpp
testDummyStack.cpp
- testSuperStupidStack.cpp
+ testSimpleStack.cpp
testNuclearStackExtension.cpp
)
diff --git a/tests/stack/testNuclearStackExtension.cpp b/tests/stack/testNuclearStackExtension.cpp
index d64ae6da9..13a3ec18c 100644
--- a/tests/stack/testNuclearStackExtension.cpp
+++ b/tests/stack/testNuclearStackExtension.cpp
@@ -22,7 +22,7 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
SECTION("write non nucleus") {
- nuclear_stack::NuclearStackExtension<SuperStupidStack,
+ nuclear_stack::NuclearStackExtension<SimpleStack,
nuclear_stack::ExtendedParticleInterfaceType>
s;
s.addParticle(std::make_tuple(
@@ -33,7 +33,7 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
SECTION("write nucleus") {
- nuclear_stack::NuclearStackExtension<SuperStupidStack,
+ nuclear_stack::NuclearStackExtension<SimpleStack,
nuclear_stack::ExtendedParticleInterfaceType>
s;
s.addParticle(std::make_tuple(
@@ -212,7 +212,7 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
SECTION("not allowed") {
- nuclear_stack::NuclearStackExtension<SuperStupidStack,
+ nuclear_stack::NuclearStackExtension<SimpleStack,
nuclear_stack::ExtendedParticleInterfaceType>
s;
diff --git a/tests/stack/testSuperStupidStack.cpp b/tests/stack/testSimpleStack.cpp
similarity index 91%
rename from tests/stack/testSuperStupidStack.cpp
rename to tests/stack/testSimpleStack.cpp
index a038d22c2..7aa764359 100644
--- a/tests/stack/testSuperStupidStack.cpp
+++ b/tests/stack/testSimpleStack.cpp
@@ -9,7 +9,7 @@
#define protected public // to also test the internal state of objects
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
-#include <corsika/stack/SuperStupidStack.hpp>
+#include <corsika/stack/SimpleStack.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <catch2/catch.hpp>
@@ -17,13 +17,13 @@
using namespace corsika;
using namespace std;
-TEST_CASE("SuperStupidStack", "[stack]") {
+TEST_CASE("SimpleStack", "[stack]") {
const CoordinateSystemPtr& dummyCS = get_root_CoordinateSystem();
SECTION("read+write") {
- SuperStupidStack s;
+ SimpleStack s;
s.addParticle(std::make_tuple(
Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
@@ -39,7 +39,7 @@ TEST_CASE("SuperStupidStack", "[stack]") {
SECTION("write+delete") {
- SuperStupidStack s;
+ SimpleStack s;
for (int i = 0; i < 99; ++i)
s.addParticle(std::make_tuple(
--
GitLab