From ad686f52f54efb12f02de6835c26567871ba8161 Mon Sep 17 00:00:00 2001
From: Ralf Ulrich <ralf.ulrich@kit.edu>
Date: Mon, 18 Jan 2021 08:22:27 +0100
Subject: [PATCH] first improve log output
---
corsika/detail/modules/ObservationPlane.inl | 133 ++++++------------
corsika/framework/core/Logging.hpp | 4 +-
.../framework/process/ContinuousProcess.hpp | 2 +-
corsika/framework/process/ProcessSequence.hpp | 36 +++--
corsika/modules/ObservationPlane.hpp | 3 +-
examples/vertical_EAS.cpp | 4 +-
tests/framework/testProcessSequence.cpp | 18 ++-
7 files changed, 91 insertions(+), 109 deletions(-)
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index e266a357f..526fb273a 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -7,6 +7,8 @@
*/
#include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
#include <fstream>
@@ -27,115 +29,64 @@ namespace corsika {
ProcessReturn ObservationPlane::doContinuous(
corsika::setup::Stack::particle_type& particle,
- corsika::setup::Trajectory& trajectory) {
- TimeType const timeOfIntersection =
- (plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
- trajectory.getVelocity(0).dot(plane_.getNormal());
+ corsika::setup::Trajectory& trajectory,
+ bool const stepLimit) {
- if (timeOfIntersection < TimeType::zero()) { return ProcessReturn::Ok; }
+ auto const* volumeNode = particle.getNode();
- if (plane_.isAbove(trajectory.getPosition(0)) ==
- plane_.isAbove(trajectory.getPosition(1))) {
- return ProcessReturn::Ok;
- }
-
- const auto energy = particle.getEnergy();
- auto const displacement = trajectory.getPosition(1) - plane_.getCenter();
-
- outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
- << ' ' << displacement.dot(xAxis_) / 1_m << ' '
- << displacement.dot(yAxis_) / 1_m
- << (trajectory.getPosition(1) - plane_.getCenter()).getNorm() / 1_m
- << std::endl;
+ if (!stepLimit) {
+ return ProcessReturn::Ok;
+ }
- if (deleteOnHit_) {
- count_ground_++;
- energy_ground_ += energy;
- particle.erase();
- return ProcessReturn::ParticleAbsorbed;
- } else {
- return ProcessReturn::Ok;
+ HEPEnergyType const energy = particle.getEnergy();
+ TimeType const timeOfIntersection = particle.getTime();
+ Point const pointOfIntersection = particle.getPosition();
+ Vector const displacement = pointOfIntersection - plane_.getCenter();
+
+ outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
+ << ' ' << displacement.dot(xAxis_) / 1_m << ' '
+ << displacement.dot(yAxis_) / 1_m
+ << (pointOfIntersection - plane_.getCenter()).getNorm() / 1_m
+ << '\n';
+
+ if (deleteOnHit_) {
+ count_ground_++;
+ energy_ground_ += energy;
+ particle.erase();
+ return ProcessReturn::ParticleAbsorbed;
+ } else {
+ return ProcessReturn::Ok;
}
}
LengthType ObservationPlane::getMaxStepLength(
- corsika::setup::Stack::particle_type const& vParticle,
+ corsika::setup::Stack::particle_type const& particle,
corsika::setup::Trajectory const& trajectory) {
- int chargeNumber;
- if (is_nucleus(vParticle.getPID())) {
- chargeNumber = vParticle.getNuclearZ();
- } else {
- chargeNumber = get_charge_number(vParticle.getPID());
- }
- auto const* currentLogicalVolumeNode = vParticle.getNode();
- auto magneticfield = currentLogicalVolumeNode->getModelProperties().getMagneticField(
- vParticle.getPosition());
- auto direction = trajectory.getVelocity(0).normalized();
-
- if (chargeNumber != 0 &&
- abs(plane_.getNormal().dot(
- trajectory.getLine().getVelocity().cross(magneticfield))) *
- 1_s / 1_m / 1_T >
- 1e-6) {
- auto const* currentLogicalVolumeNode = vParticle.getNode();
- auto magneticfield =
- currentLogicalVolumeNode->getModelProperties().getMagneticField(
- vParticle.getPosition());
- auto k =
- chargeNumber * constants::c * 1_eV / (vParticle.getMomentum().getNorm() * 1_V);
-
- if (direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k) <
- 0) {
- return std::numeric_limits<double>::infinity() * 1_m;
- }
+ auto const& volumeNode = particle.getNode();
- LengthType MaxStepLength1 =
- (sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane_.getNormal()) / direction.getNorm()) /
- (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
- LengthType MaxStepLength2 =
- (-sqrt(direction.dot(plane_.getNormal()) * direction.dot(plane_.getNormal()) -
- (plane_.getNormal().dot(trajectory.getPosition(0) - plane_.getCenter()) *
- plane_.getNormal().dot(direction.cross(magneticfield)) * 2 * k)) -
- direction.dot(plane_.getNormal()) / direction.getNorm()) /
- (plane_.getNormal().dot(direction.cross(magneticfield)) * k);
-
- if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
- return std::numeric_limits<double>::infinity() * 1_m;
- } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
- std::cout << " steplength to obs plane 2: " << MaxStepLength2 << std::endl;
- return MaxStepLength2 *
- (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
- .getNorm() *
- 1.001;
- } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
- std::cout << " steplength to obs plane 1: " << MaxStepLength1 << std::endl;
- return MaxStepLength1 *
- (direction + direction.cross(magneticfield) * MaxStepLength2 * k / 2)
- .getNorm() *
- 1.001;
- }
- }
+ typedef typename std::remove_const_t<
+ std::remove_reference_t<decltype(volumeNode->getModelProperties())>>
+ medium_type;
+
+ Intersections const intersection = setup::Tracking::intersect<corsika::setup::Stack::particle_type, medium_type>(particle, plane_,
+ volumeNode->getModelProperties());
- TimeType const timeOfIntersection =
- (plane_.getCenter() - trajectory.getPosition(0)).dot(plane_.getNormal()) /
- trajectory.getVelocity(0).dot(plane_.getNormal());
- if (timeOfIntersection < TimeType::zero()) {
+ TimeType const timeOfIntersection = intersection.getEntry();
+
+ if (timeOfIntersection < TimeType::zero()) {
+ return std::numeric_limits<double>::infinity() * 1_m;
+ }
+ if (timeOfIntersection > trajectory.getDuration()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
- auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm() * 1.0001;
- CORSIKA_LOG_TRACE("ObservationPlane w/o magnetic field: getMaxStepLength l={} m",
+ auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm() * 1.001; // make max. step length a bit longer, to assure crossing
+ CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m",
dist / 1_m);
return dist;
}
diff --git a/corsika/framework/core/Logging.hpp b/corsika/framework/core/Logging.hpp
index 4d4c43a7b..29ef3e13c 100644
--- a/corsika/framework/core/Logging.hpp
+++ b/corsika/framework/core/Logging.hpp
@@ -35,9 +35,9 @@
// if this is a Debug build, include debug messages in objects
#ifdef DEBUG
// trace is the highest level of logging (ALL messages will be printed)
-#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_DEBUG
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_TRACE
#else // otherwise, remove everything but "error" and worse messages
-#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_INFO
+#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_DEBUG
#endif
#include <spdlog/fmt/ostr.h> // will output whenerver a streaming operator is found
diff --git a/corsika/framework/process/ContinuousProcess.hpp b/corsika/framework/process/ContinuousProcess.hpp
index 1343c0822..b5ba72f98 100644
--- a/corsika/framework/process/ContinuousProcess.hpp
+++ b/corsika/framework/process/ContinuousProcess.hpp
@@ -30,7 +30,7 @@ namespace corsika {
// here starts the interface part
// -> enforce TDerived to implement DoContinuous...
template <typename TParticle, typename TTrack>
- ProcessReturn doContinuous(TParticle&, TTrack const&) const;
+ ProcessReturn doContinuous(TParticle&, TTrack const&, bool const stepLimit) const;
// -> enforce TDerived to implement MaxStepLength...
template <typename TParticle, typename TTrack>
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 70337e109..1758f23d1 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -26,6 +26,21 @@
namespace corsika {
+ namespace detail {
+ template <typename TProcess, int N=0, typename Enable=void>
+ struct CountContinuous {
+ enum { count = N };
+ };
+
+ template <typename TProcess, int N=0>
+ struct CountContinuous<TProcess, N,
+ typename std::enable_if_t<std::is_base_of_v<ContinuousProcess<typename std::decay_t<TProcess>>>> {
+ enum { count = N+1 };
+ };
+ }
+
+
+
/**
*
* Definition of a static process list/sequence
@@ -39,13 +54,12 @@ namespace corsika {
* they are just classes. This allows us to handle both, rvalue as
* well as lvalue Processes in the ProcessSequence.
*
- * The sequence, and the processes use CRTP.
+ * (The sequence, and the processes use the CRTP, curiously recurring template pattern).
*
- * \todo There are several FIXME's in the ProcessSequence.inl due to
- * outstanding migration of SecondaryView::parent()
**/
- template <typename TProcess1, typename TProcess2 = NullModel>
+ template <typename TProcess1, typename TProcess2 = NullModel,
+ int NContinuous = detail::CountContinuous<TProcess1, detail::CountContinous<TProcess1>::count>::count>
class ProcessSequence : public BaseProcess<ProcessSequence<TProcess1, TProcess2>> {
using process1_type = typename std::decay_t<TProcess1>;
@@ -57,6 +71,9 @@ namespace corsika {
static bool constexpr t1SwitchProcSeq = is_switch_process_sequence_v<process1_type>;
static bool constexpr t2SwitchProcSeq = is_switch_process_sequence_v<process2_type>;
+ // we want to count continuous processes, to each one has a unique index
+ if constexpr (std::is_base_of_v<ContinuousProcess<typename std::decay_t<TProcess1>> )
+
// make sure only BaseProcess types TProcess1/2 are passed
static_assert(std::is_base_of_v<BaseProcess<process1_type>, process1_type>,
"can only use process derived from BaseProcess in "
@@ -65,9 +82,6 @@ namespace corsika {
"can only use process derived from BaseProcess in "
"ProcessSequence, for Process 2");
- TProcess1 A_; /// process/list A, this is a reference, if possible
- TProcess2 B_; /// process/list B, this is a reference, if possible
-
public:
// resource management
ProcessSequence() = delete; // only initialized objects
@@ -118,7 +132,7 @@ namespace corsika {
inline void doStack(TStack& stack);
template <typename TParticle, typename TTrack>
- inline LengthType getMaxStepLength(TParticle& particle, TTrack& vTrack);
+ inline std::pair<LengthType,void*> getMaxStepLength(TParticle& particle, TTrack& vTrack);
template <typename TParticle>
inline GrammageType getInteractionLength(TParticle&& particle) {
@@ -147,6 +161,12 @@ namespace corsika {
inline ProcessReturn selectDecay(
TSecondaryView& view, [[maybe_unused]] InverseTimeType decay_inv_select,
[[maybe_unused]] InverseTimeType decay_inv_sum = InverseTimeType::zero());
+
+ private:
+ int countContinuous_ = Ncontinuous;
+ TProcess1 A_; /// process/list A, this is a reference, if possible
+ TProcess2 B_; /// process/list B, this is a reference, if possible
+
};
/**
diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index 7effcd93b..8276c0d12 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -31,7 +31,8 @@ namespace corsika {
bool = true);
ProcessReturn doContinuous(corsika::setup::Stack::particle_type& vParticle,
- corsika::setup::Trajectory& vTrajectory);
+ corsika::setup::Trajectory& vTrajectory,
+ bool const stepLimit);
LengthType getMaxStepLength(corsika::setup::Stack::particle_type const&,
corsika::setup::Trajectory const& vTrajectory);
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 02ec3a52b..4ebffab2e 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -91,8 +91,8 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
- logging::set_level(logging::level::info);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
+ logging::set_level(logging::level::trace);
CORSIKA_LOG_INFO("vertical_EAS");
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index ac80847ca..cf414de5f 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -247,10 +247,10 @@ struct DummyView {
DummyData& parent() { return p_; }
};
-TEST_CASE("Process Sequence", "[Process Sequence]") {
+TEST_CASE("Process Sequence General", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
SECTION("Check construction") {
globalCount = 0;
@@ -373,10 +373,10 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
}
}
-TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
+TEST_CASE("Switch Process Sequence", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
SECTION("Check construction") {
@@ -479,3 +479,13 @@ TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
CHECK(checkSec == 0);
}
}
+
+TEST_CASE("Continuous Process Indexing", "ProcessSequence") {
+
+ logging::set_level(logging::level::info);
+ corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
+
+ SECTION("Check construction") {
+
+ }
+}
--
GitLab