diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 464ce2402377ddacbf94a9147aa958ebbd9bdf3d..9ebe0555b9b6f4745c25e3e1bf09a3c1723a05eb 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -170,7 +170,7 @@ build-clang-8:
script:
- cd build
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
rules:
- if: $CI_MERGE_REQUEST_ID
when: manual
@@ -185,8 +185,6 @@ build-clang-8:
reports:
junit:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
- paths:
- - ${CI_PROJECT_DIR}/build/test.log.gz
cache:
paths:
- ${CI_PROJECT_DIR}/build/
@@ -227,7 +225,7 @@ test-clang-8:
- cd build
- cmake --build . -- -j4
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
rules:
- if: '$CI_MERGE_REQUEST_ID && $CI_MERGE_REQUEST_TITLE =~ /^Draft:/'
allow_failure: false
@@ -241,8 +239,6 @@ test-clang-8:
reports:
junit:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
- paths:
- - ${CI_PROJECT_DIR}/build/test.log.gz
cache:
paths:
- ${CI_PROJECT_DIR}/build/
@@ -338,7 +334,7 @@ example-clang-8:
- cd build
- cmake --build . -- -j4
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
- make -j4 run_examples | gzip -v -9 > examples.log.gz
rules:
- if: '$CI_MERGE_REQUEST_ID && $CI_MERGE_REQUEST_TITLE =~ /^Draft:/'
@@ -357,7 +353,6 @@ example-clang-8:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
paths:
- ${CI_PROJECT_DIR}/build/examples.log.gz
- - ${CI_PROJECT_DIR}/build/test.log.gz
cache:
paths:
- ${CI_PROJECT_DIR}/build/
@@ -450,7 +445,7 @@ install-clang-8:
- cmake .. -DCMAKE_BUILD_TYPE=Release
- cmake --build . -- -j4
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
- make -j4 run_examples | gzip -v -9 > examples.log.gz
rules:
- if: '$CI_MERGE_REQUEST_LABELS =~ /Ready for code review/' # run on merge requests, if label 'Ready for code review' is set
@@ -475,7 +470,6 @@ install-clang-8:
junit:
- ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml
paths:
- - ${CI_PROJECT_DIR}/build/test.log.gz
- ${CI_PROJECT_DIR}/build/examples.log.gz
# release for gcc
@@ -513,7 +507,7 @@ coverage:
- cd build
- cmake .. -DCMAKE_BUILD_TYPE=Coverage
- cmake --build . -- -j4
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4
- cmake --build . --target coverage
- tar czf coverage-report.tar.gz coverage-report
coverage: '/^.*functions\.+:\s(.*\%)\s/'
@@ -530,7 +524,6 @@ coverage:
when: always
expire_in: 1 year
paths:
- - ${CI_PROJECT_DIR}/build/test.log.gz
- ${CI_PROJECT_DIR}/build/coverage-report.tar.gz
cache:
paths:
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 46ecfb3cb0afa503e7ed174c003ab13b45f26a2e..3231ab4a8a643b682a50723e9b64dc3837e1c761 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -19,6 +19,7 @@ set (
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
+ IEmpty.hpp
IMediumModel.h
NuclearComposition.h
HomogeneousMedium.h
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 774b7492761469affbfc2d6bf326cac5501dbbf2..1f807de45bc698718c74f93cbc7021c8329936b2 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -8,7 +8,6 @@
#pragma once
-#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 8cad83555833331e3d3c57912379f33ca3ac8eee..ad3306810d4efccbee8905d5d281785e6180880d 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -8,16 +8,14 @@
#pragma once
-#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/IEmpty.hpp>
#include <corsika/geometry/Volume.h>
#include <memory>
#include <vector>
namespace corsika::environment {
- class Empty {}; //<! intended for usage as default template argument
-
- template <typename TModelProperties = Empty>
+ template <typename TModelProperties = IEmpty>
class VolumeTreeNode {
public:
using IModelProperties = TModelProperties;
@@ -116,14 +114,15 @@ namespace corsika::environment {
void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
+ /*
template <class MediumType, typename... Args>
static auto CreateMedium(Args&&... args) {
- static_assert(std::is_base_of_v<IMediumModel, MediumType>,
+ static_assert(std::is_base_of_v<, MediumType>,
"unusable type provided, needs to be derived from \"IMediumModel\"");
return std::make_shared<MediumType>(std::forward<Args>(args)...);
}
-
+ */
private:
std::vector<VTNUPtr> fChildNodes;
std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 3e81ef973e333f1f3d264810f0f1a88138e5acfa..6dc602de24ba9b03c475c721e2e178aa99296412 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -19,8 +19,6 @@
#include <corsika/stack/history/EventType.hpp>
#include <corsika/stack/history/HistorySecondaryProducer.hpp>
-#include <corsika/setup/SetupTrajectory.h>
-
/* see Issue 161, we need to include SetupStack only because we need
to globally define StackView. This is clearly not nice and should
be changed, when possible. It might be that StackView needs to be
@@ -174,8 +172,9 @@ namespace corsika::cascade {
// assert that particle stays outside void Universe if it has no
// model properties set
- assert(currentLogicalNode != &*environment_.GetUniverse() ||
- environment_.GetUniverse()->HasModelProperties());
+ assert((currentLogicalNode != &*environment_.GetUniverse() ||
+ environment_.GetUniverse()->HasModelProperties()) &&
+ "FATAL: The environment model has no valid properties set!");
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime =
@@ -237,9 +236,9 @@ namespace corsika::cascade {
}
C8LOG_DEBUG("sth. happening before geometric limit ? {}",
- ((min_distance < geomMaxLength) ? "yes" : "no"));
+ ((min_distance <= geomMaxLength) ? "yes" : "no"));
- if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
+ if (min_distance <= geomMaxLength) { // interaction to happen within geometric limit
// check whether decay or interaction limits this step the
// outcome of decay or interaction MAY be a) new particles in
@@ -333,7 +332,7 @@ namespace corsika::cascade {
const auto sample_process = uniDist(rng_);
auto const returnCode = process_sequence_.SelectInteraction(view, sample_process);
if (returnCode != process::EProcessReturn::eInteracted) {
- C8LOG_WARN("Particle did not interace!");
+ C8LOG_WARN("Particle did not interact!");
}
SetEventType(view, history::EventType::Interaction);
return returnCode;
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 1018bb002b63885244a067a76746de34b5a02d3c..2ca4dd7f88d8fd5e97aad9ddab1f2750ad96377b 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -13,8 +13,6 @@
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/NullModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/Tracking.h>
-//#include <corsika/process/tracking_bfield/Tracking.h>
#include <corsika/particles/ParticleProperties.h>
@@ -36,44 +34,61 @@ using namespace corsika::geometry;
#include <limits>
using namespace std;
-/*
- The dummy env must support GetMagneticField(), and a density model
+/**
+ * testCascade implements an e.m. Heitler model with energy splitting
+ * and a critical energy.
+ *
+ * It resembles one of the most simple cascades you can simulate with CORSIKA8.
+ **/
+/*
+ The dummy env (here) doesn't need to have any propoerties
*/
auto MakeDummyEnv() {
TestEnvironmentType env; // dummy environment
auto& universe = *(env.GetUniverse());
- const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
- auto theMedium = TestEnvironmentType::CreateNode<Sphere>(
+ auto world = TestEnvironmentType::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_m * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel = environment::UniformMagneticField<
- environment::HomogeneousMedium<TestEnvironmentInterface>>;
+ using MyEmptyModel = environment::Empty<environment::IEmpty>;
+ world->SetModelProperties<MyEmptyModel>();
- theMedium->SetModelProperties<MyHomogeneousModel>(
- geometry::Vector(cs, 0_T, 0_T, 0_T), 1_g / (1_cm * 1_cm * 1_cm),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
-
- universe.AddChild(std::move(theMedium));
+ universe.AddChild(std::move(world));
return env;
}
+/**
+ * \class DummyTracking
+ *
+ * For the Heitler model we don't need particle transport.
+ **/
+class DummyTracking {
+
+public:
+ template <typename TParticle>
+ auto GetTrack(TParticle const& particle) {
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ geometry::Vector<SpeedType::dimension_type> const initialVelocity =
+ particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
+ return std::make_tuple(
+ geometry::LineTrajectory(geometry::Line(particle.GetPosition(), initialVelocity),
+ 0_s), // trajectory
+ particle.GetNode()); // next volume node
+ }
+};
+
class ProcessSplit : public process::InteractionProcess<ProcessSplit> {
int fCalls = 0;
- GrammageType fX0;
public:
- ProcessSplit(GrammageType const X0)
- : fX0(X0) {}
-
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle const&) const {
- return fX0;
+ return 0_g / square(1_cm);
}
template <typename TSecondaryView>
@@ -126,6 +141,8 @@ public:
TEST_CASE("Cascade", "[Cascade]") {
+ logging::SetLevel(logging::level::trace);
+
HEPEnergyType E0 = 100_GeV;
random::RNGManager& rmng = random::RNGManager::GetInstance();
@@ -133,14 +150,12 @@ TEST_CASE("Cascade", "[Cascade]") {
auto env = MakeDummyEnv();
auto const& rootCS = env.GetCoordinateSystem();
- tracking_line::Tracking tracking;
stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0);
process::NullModel nullModel;
- const GrammageType X0 = 20_g / square(1_cm);
const HEPEnergyType Ecrit = 85_MeV;
- ProcessSplit split(X0);
+ ProcessSplit split;
ProcessCut cut(Ecrit);
auto sequence = process::sequence(nullModel, stackInspect, split, cut);
TestCascadeStack stack;
@@ -153,7 +168,8 @@ TEST_CASE("Cascade", "[Cascade]") {
particles::GetMass(particles::Code::Electron)))}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
- cascade::Cascade<tracking_line::Tracking, decltype(sequence), TestCascadeStack,
+ DummyTracking tracking;
+ cascade::Cascade<DummyTracking, decltype(sequence), TestCascadeStack,
TestCascadeStackView>
EAS(env, tracking, sequence, stack);
@@ -161,7 +177,7 @@ TEST_CASE("Cascade", "[Cascade]") {
EAS.Run();
CHECK(cut.GetCount() == 2048);
- CHECK(cut.GetCalls() == 2047);
+ CHECK(cut.GetCalls() == 2047); // final particle is still on stack and not yet deleted
CHECK(split.GetCalls() == 2047);
}
diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h
index 848c597ee3985e15aa2cfdad956b1dc621e2ca21..e057441433a58f24e374ee925497fc445ab3cae9 100644
--- a/Framework/Cascade/testCascade.h
+++ b/Framework/Cascade/testCascade.h
@@ -9,15 +9,14 @@
#pragma once
#include <corsika/environment/Environment.h>
-#include <corsika/environment/IMagneticFieldModel.h>
+#include <corsika/environment/IEmpty.hpp>
#include <corsika/stack/CombinedStack.h>
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/node/GeometryNodeStackExtension.h>
#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
-using TestEnvironmentInterface =
- corsika::environment::IMagneticFieldModel<corsika::environment::IMediumModel>;
+using TestEnvironmentInterface = corsika::environment::IEmpty;
using TestEnvironmentType = corsika::environment::Environment<TestEnvironmentInterface>;
template <typename T>
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 32d3025940e78be21d15d58f5593f32d4a8f514e..02a3d062ea09792f9c9e4f40d030bd7fc3c02414 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -215,7 +215,7 @@ namespace corsika::process {
const auto particle = view.parent();
lambda_inv_sum += A_.GetInverseInteractionLength(particle);
// check if we should execute THIS process and then EXIT
- if (lambda_inv_select < lambda_inv_sum) {
+ if (lambda_inv_select <= lambda_inv_sum) {
A_.DoInteraction(view);
return EProcessReturn::eInteracted;
}
@@ -229,7 +229,7 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B_.GetInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
- if (lambda_inv_select < lambda_inv_sum) {
+ if (lambda_inv_select <= lambda_inv_sum) {
B_.DoInteraction(view);
return EProcessReturn::eInteracted;
}
@@ -279,8 +279,8 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A_.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
- if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select <
- // decay_inv_sum / decay_inv_tot
+ if (decay_inv_select <= decay_inv_sum) { // more pedagogical: rndm_select <
+ // decay_inv_sum / decay_inv_tot
A_.DoDecay(view);
return EProcessReturn::eDecayed;
}
@@ -294,7 +294,7 @@ namespace corsika::process {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B_.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
- if (decay_inv_select < decay_inv_sum) {
+ if (decay_inv_select <= decay_inv_sum) {
B_.DoDecay(view);
return EProcessReturn::eDecayed;
}
diff --git a/Processes/Tracking/Intersect.hpp b/Processes/Tracking/Intersect.hpp
index 277b66f0252931d82b468c781509938eba141d36..277c5f6336478ce8dd2cbf0c84290cd72be8b6a0 100644
--- a/Processes/Tracking/Intersect.hpp
+++ b/Processes/Tracking/Intersect.hpp
@@ -14,12 +14,25 @@
#include <corsika/logging/Logging.h>
#include <corsika/geometry/Intersections.hpp>
+#include <limits>
+
namespace corsika::process::tracking {
/**
* \class Intersect
*
+ * This is a CRTP class to provide a generic volume-tree
+ * intersection for the purpose of tracking.
+ *
+ * It return the closest distance in time to the next geometric
+ * intersection, as well as a pointer to the corresponding new
+ * volume.
*
+ * User may provide an optional global step-length limit as
+ * parameter. This may be needd for (simpler) algorithms in magnetic
+ * fields, where tracking errors grow linearly with step-length.
+ * Obviously, in the case of the step-length limit, the returend
+ * "next" volume is just the current one.
*
**/
@@ -28,45 +41,23 @@ namespace corsika::process::tracking {
protected:
template <typename TParticle>
- auto nextIntersect(const TParticle& particle) const {
+ auto nextIntersect(
+ const TParticle& particle,
+ corsika::units::si::TimeType step_limit =
+ std::numeric_limits<corsika::units::si::TimeType::value_type>::infinity() *
+ corsika::units::si::second) const {
using namespace corsika::units::si;
using namespace corsika::geometry;
- const Point& initialPosition = particle.GetPosition();
-
typedef
typename std::remove_reference<decltype(*particle.GetNode())>::type node_type;
- node_type& volumeNode = *particle.GetNode();
+ node_type& volumeNode =
+ *particle.GetNode(); // current "logical" node, from previous tracking step
C8LOG_DEBUG("volumeNode={}, numericallyInside={} ", fmt::ptr(&volumeNode),
- volumeNode.GetVolume().Contains(initialPosition));
-
- auto const velocity =
- particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c;
-
- // 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
- const auto& magneticfield =
- volumeNode.GetModelProperties().GetMagneticField(initialPosition);
- const auto magnitudeB = magneticfield.norm();
- const int chargeNumber = particle.GetChargeNumber();
- auto const momentumVerticalMag =
- particle.GetMomentum() -
- particle.GetMomentum().parallelProjectionOnto(magneticfield);
- LengthType const gyroradius =
- (chargeNumber == 0 || magnitudeB == 0_T
- ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
- : momentumVerticalMag.norm() * 1_V /
- (corsika::units::constants::c * abs(chargeNumber) * magnitudeB *
- 1_eV));
- const double maxRadians = 0.01;
- const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
- C8LOG_DEBUG("gyroradius {}, steplimit: {} m = {} s", gyroradius, steplimit,
- steplimit / velocity.norm());
+ volumeNode.GetVolume().Contains(particle.GetPosition()));
// start values:
- TimeType minTime = steplimit / velocity.norm();
+ TimeType minTime = step_limit;
node_type* minNode = &volumeNode;
// determine the first geometric collision with any other Volume boundary
diff --git a/Processes/TrackingLeapFrogCurved/Tracking.h b/Processes/TrackingLeapFrogCurved/Tracking.h
index ee0ba6e7f950e5be656b2d2c3b42dc30e1c7dc5b..4ffa438289aa093946adf320a16e75a288dbe5eb 100644
--- a/Processes/TrackingLeapFrogCurved/Tracking.h
+++ b/Processes/TrackingLeapFrogCurved/Tracking.h
@@ -107,13 +107,34 @@ namespace corsika::process {
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
+ const auto& magneticfield =
+ volumeNode.GetModelProperties().GetMagneticField(position);
+ const auto magnitudeB = magneticfield.norm();
+ const int chargeNumber = particle.GetChargeNumber();
+ auto const momentumVerticalMag =
+ particle.GetMomentum() -
+ particle.GetMomentum().parallelProjectionOnto(magneticfield);
+ LengthType const gyroradius =
+ (chargeNumber == 0 || magnitudeB == 0_T
+ ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m
+ : momentumVerticalMag.norm() * 1_V /
+ (corsika::units::constants::c * abs(chargeNumber) * magnitudeB *
+ 1_eV));
+ const double maxRadians = 0.01;
+ const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
+ const TimeType steplimit_time = steplimit / initialVelocity.norm();
+ C8LOG_DEBUG("gyroradius {}, steplimit: {} m = {} s", gyroradius, steplimit,
+ steplimit_time);
+
// traverse the environment volume tree and find next
// intersection
- auto [minTime, minNode] = tracking::Intersect<Tracking>::nextIntersect(particle);
+ auto [minTime, minNode] =
+ tracking::Intersect<Tracking>::nextIntersect(particle, steplimit_time);
- const int chargeNumber = particle.GetChargeNumber();
- const auto& magneticfield =
- volumeNode.GetModelProperties().GetMagneticField(position);
const auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV /
(particle.GetEnergy() * 1_V);
return std::make_tuple(
diff --git a/Processes/TrackingLeapFrogStraight/Tracking.h b/Processes/TrackingLeapFrogStraight/Tracking.h
index fd15b89c87fbff03977ee9074686fa78704c7909..d9e7c57de709866bf410a8050016d70bdf62ebd3 100644
--- a/Processes/TrackingLeapFrogStraight/Tracking.h
+++ b/Processes/TrackingLeapFrogStraight/Tracking.h
@@ -39,14 +39,14 @@ namespace corsika::process {
* tracking_line::Tracking and adds a (two-step) Leap-Frog
* algorithms with two halve-steps and magnetic deflection.
*
- * The two halve steps are implemented as two
+ * The two halve steps are implemented as two straight explicit
* `tracking_line::Tracking`s and all geometry intersections are,
* thus, based on those two straight line elements.
*
* As a precaution for numerical instability, the steplength is
* limited to correspond to a straight line distance to the next
* volume intersection. In typical situations this leads to about
- * one full leap-frog step to the next volume boundary.
+ * (at least) one full leap-frog step to the next volume boundary.
*
**/