From c38635469d52640288615e1ca1265d20815de468 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 14 Feb 2019 17:13:47 +0100
Subject: [PATCH] first try
---
Environment/VolumeTreeNode.h | 1 +
Framework/Cascade/Cascade.h | 59 ++++++-----
Processes/TrackingLine/TrackingLine.cc | 67 -------------
Processes/TrackingLine/TrackingLine.h | 97 ++++++++++++++++++-
Processes/TrackingLine/testTrackingLine.cc | 2 +-
.../TrackingLine/testTrackingLineStack.cc | 37 -------
6 files changed, 130 insertions(+), 133 deletions(-)
delete mode 100644 Processes/TrackingLine/testTrackingLineStack.cc
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 3331c41bf..e125b243d 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -13,6 +13,7 @@
#define _include_VolumeTreeNode_H
#include <corsika/geometry/Volume.h>
+#include <corsika/environment/IMediumModel.h>
#include <memory>
#include <vector>
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 2d55e9030..4bddb9c57 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -14,9 +14,9 @@
#include <corsika/environment/Environment.h>
#include <corsika/process/ProcessReturn.h>
+#include <corsika/random/ExponentialDistribution.h>
#include <corsika/random/RNGManager.h>
#include <corsika/random/UniformRealDistribution.h>
-#include <corsika/random/ExponentialDistribution.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
@@ -37,7 +37,7 @@ namespace corsika::cascade {
* plugged into the cascade simulation.
*
* <b>Tracking</b> must be a class according to the
- * TrackingInterface providing the functions:
+ * TrackingInterface providing the functions:
* <code>auto GetTrack(Particle const& p)</auto>,
* with the return type <code>geometry::Trajectory<corsika::geometry::Line>
* </code>
@@ -77,19 +77,19 @@ namespace corsika::cascade {
fProcessSequence.Init();
fStack.Init();
}
-
+
/**
* set the nodes for all particles on the stack according to their numerical
* position
*/
void SetNodes() {
- std::for_each(fStack.begin(), fStack.end(), [&](auto& p) {
- auto const* numericalNode =
+ std::for_each(fStack.begin(), fStack.end(), [&](auto& p) {
+ auto const* numericalNode =
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
- p.SetNode(numericalNode);
-
- std::cout << "initial node " << p.GetNode() << std::endl;
- });
+ p.SetNode(numericalNode);
+
+ std::cout << "initial node " << p.GetNode() << std::endl;
+ });
}
/**
@@ -98,11 +98,12 @@ namespace corsika::cascade {
*/
void Run() {
SetNodes();
-
- while (!fStack.IsEmpty()) {
- while (!fStack.IsEmpty()) {
+
+ while (!fStack.IsEmpty() && countSteps < maxSteps) {
+ while (!fStack.IsEmpty() && countSteps < maxSteps) {
auto pNext = fStack.GetNextParticle();
Step(pNext);
+ countSteps++;
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
@@ -125,7 +126,7 @@ namespace corsika::cascade {
using namespace corsika::units::si;
// determine geometric tracking
- corsika::setup::Trajectory step = fTracking.GetTrack(particle);
+ auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(particle);
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
@@ -139,12 +140,13 @@ namespace corsika::cascade {
<< ", next_interact=" << next_interact << std::endl;
auto const* currentLogicalNode = particle.GetNode();
-
+
auto const* currentNumericalNode =
- fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
-
- std::cout << "nodes: " << currentLogicalNode << " " << currentNumericalNode << std::endl;
-
+ fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+
+ std::cout << "nodes: " << currentLogicalNode << " " << currentNumericalNode
+ << std::endl;
+
if (currentNumericalNode != currentLogicalNode) {
throw std::runtime_error("numerical and logical nodes don't match");
}
@@ -155,7 +157,8 @@ namespace corsika::cascade {
// convert next_step from grammage to length
LengthType const distance_interact =
- currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step, next_interact);
+ currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
+ next_interact);
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
@@ -178,7 +181,7 @@ namespace corsika::cascade {
// take minimum of geometry, interaction, decay for next step
auto const min_distance =
- std::min({distance_interact, distance_decay, distance_max});
+ std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
std::cout << " move particle by : " << min_distance << std::endl;
@@ -191,7 +194,8 @@ namespace corsika::cascade {
// particle.GetNode(); // previous VolumeNode
//~ particle.SetNode(
- //~ currentLogicalNode); // NOTE @Max : here we need to distinguish: IF particle step is
+ //~ currentLogicalNode); // NOTE @Max : here we need to distinguish: IF particle
+ //step is
// limited by tracking (via fTracking.GetTrack()), THEN we need
// to check/update VolumeNodes. In all other cases it is
// guaranteed that we are still in the same volume
@@ -208,11 +212,9 @@ namespace corsika::cascade {
}
std::cout << "sth. happening before geometric limit ? "
- << ((min_distance < distance_max) ? "yes" : "no") << std::endl;
-
- if (min_distance < distance_max) { // interaction to happen within geometric limit
- // check whether decay or interaction limits this step
+ << ((min_distance < geomMaxLength) ? "yes" : "no") << std::endl;
+ if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
@@ -225,7 +227,7 @@ namespace corsika::cascade {
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
fProcessSequence.SelectInteraction(particle, step, fStack, sample_process,
inv_lambda_count);
- } else {
+ } else if (min_distance == distance_decay) {
std::cout << "decay" << std::endl;
InverseTimeType const actual_decay_time =
fProcessSequence.GetTotalInverseLifetime(particle);
@@ -235,7 +237,12 @@ namespace corsika::cascade {
const auto sample_process = uniDist(fRNG);
InverseTimeType inv_decay_count = 0 / second;
fProcessSequence.SelectDecay(particle, fStack, sample_process, inv_decay_count);
+ } else { // step-length limitation within volume
+ std::cout << "step-length limitation" << std::endl;
}
+ } else { // boundary crossing
+ std::cout << "boundary crossing! next node = " << nextVol << std::endl;
+ particle.SetNode(nextVol);
}
}
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index 36a267a92..2c610c9c2 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -62,73 +62,6 @@ namespace corsika::process::tracking_line {
corsika::environment::Environment const& pEnv)
: fEnvironment(pEnv) {}
- template <class Stack, class Trajectory>
- Trajectory TrackingLine<Stack, Trajectory>::GetTrack(Particle const& p) {
- using std::cout;
- using std::endl;
- using namespace corsika::units::si;
- using namespace corsika::geometry;
- geometry::Vector<SpeedType::dimension_type> const velocity =
- p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
-
- auto const currentPosition = p.GetPosition();
- std::cout << "TrackingLine pid: " << p.GetPID() << " , E = " << p.GetEnergy() / 1_GeV
- << " GeV" << std::endl;
- std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates() << std::endl;
- std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
- << " GeV " << std::endl;
- std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
-
- // to do: include effect of magnetic field
- geometry::Line line(currentPosition, velocity);
-
- auto const* currentVolumeNode =
- fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
- auto const& children = currentVolumeNode->GetChildNodes();
- auto const& excluded = currentVolumeNode->GetExcludedNodes();
-
- std::vector<TimeType> intersectionTimes;
-
- auto addIfIntersects = [&](auto& vtn) {
- auto const& volume = vtn.GetVolume();
- auto const& sphere = dynamic_cast<geometry::Sphere const&>(
- volume); // for the moment we are a bit bold here and assume
- // everything is a sphere, crashes with exception if not
-
- if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
- auto const [t1, t2] = *opt;
- if (t1.magnitude() >= 0)
- intersectionTimes.push_back(t1);
- else if (t2.magnitude() >= 0)
- intersectionTimes.push_back(t2);
- }
- };
-
- for (auto const& child : children) { addIfIntersects(*child); }
-
- for (auto const* child : excluded) { addIfIntersects(*child); }
-
- addIfIntersects(*currentVolumeNode);
-
- auto const minIter =
- std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
-
- TimeType min;
-
- if (minIter == intersectionTimes.cend()) {
- min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
- // to handle the numerics properly
- //~ throw std::runtime_error("no intersection with anything!");
- } else {
- min = *minIter;
- }
-
- std::cout << " t-intersect: " << min << std::endl;
-
- return Trajectory(line, min);
- }
-
} // namespace corsika::process::tracking_line
#include <corsika/setup/SetupStack.h>
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index bd0aff7f2..a6d749d82 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -12,9 +12,12 @@
#ifndef _include_corsika_processes_TrackingLine_h_
#define _include_corsika_processes_TrackingLine_h_
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/units/PhysicalUnits.h>
-#include <map> // for std::pair
#include <optional>
+#include <type_traits>
+#include <utility>
namespace corsika::environment {
class Environment;
@@ -42,7 +45,97 @@ namespace corsika::process {
TrackingLine(corsika::environment::Environment const& pEnv);
- Trajectory GetTrack(Particle const& p);
+ auto GetTrack(Particle const& p) {
+ using std::cout;
+ using std::endl;
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ geometry::Vector<SpeedType::dimension_type> const velocity =
+ p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
+
+ auto const currentPosition = p.GetPosition();
+ std::cout << "TrackingLine pid: " << p.GetPID()
+ << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+ std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
+ << std::endl;
+ std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+ std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
+ << " GeV " << std::endl;
+ std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+
+ // to do: include effect of magnetic field
+ geometry::Line line(currentPosition, velocity);
+
+ auto const* currentLogicalVolumeNode = p.GetNode();
+ //~ auto const* currentNumericalVolumeNode =
+ //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+ auto const numericallyInside =
+ currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
+
+ auto const& children = currentLogicalVolumeNode->GetChildNodes();
+ auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
+
+ std::vector<std::pair<TimeType, environment::BaseNodeType const*>> intersections;
+
+ auto addIfIntersects = [&](auto const& vtn, auto const& nextNode) {
+ static_assert(std::is_same_v<decltype(vtn), decltype(nextNode)>);
+ auto const& volume = vtn.GetVolume();
+ auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+ volume); // for the moment we are a bit bold here and assume
+ // everything is a sphere, crashes with exception if not
+
+ if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
+ auto const [t1, t2] = *opt;
+ std::cout << "intersection times: " << t1 / 1_s << "; " << t2 / 1_s
+ << std::endl;
+ if (t1.magnitude() > 0)
+ intersections.emplace_back(t1, &nextNode);
+ else if (t2.magnitude() > 0)
+ intersections.emplace_back(t2, &nextNode);
+ }
+ };
+
+ for (auto const& child : children) { addIfIntersects(*child, *child); }
+
+ for (auto const* ex : excluded) { addIfIntersects(*ex, *ex); }
+
+ if (numericallyInside) {
+ addIfIntersects(
+ *currentLogicalVolumeNode,
+ *currentLogicalVolumeNode
+ ->GetParent()); // todo: add parent node to vector, not current!
+ } else {
+ auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+ currentLogicalVolumeNode->GetVolume());
+ // for the moment we are a bit bold here and assume
+ // everything is a sphere, crashes with exception if not
+ auto const [t1, t2] = *TimeOfIntersection(line, sphere);
+
+ if (t1 > 0_s) {
+ assert(t2 > 0_s);
+ intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
+ }
+ }
+
+ auto const minIter = std::min_element(
+ intersections.cbegin(), intersections.cend(),
+ [](auto const& a, auto const& b) { return a.first < b.first; });
+
+ TimeType min;
+
+ if (minIter == intersections.cend()) {
+ min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+ // to handle the numerics properly
+ throw std::runtime_error("no intersection with anything!");
+ } else {
+ min = minIter->first;
+ }
+
+ std::cout << " t-intersect: " << min << std::endl;
+
+ return std::make_tuple(Trajectory(line, min), velocity.norm() * min,
+ minIter->second);
+ }
};
} // namespace tracking_line
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index d0c225999..19efb14f3 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -85,7 +85,7 @@ TEST_CASE("TrackingLine") {
0_m / second, 1_m / second);
Line line(origin, v);
- auto const traj = tracking.GetTrack(p);
+ auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
.GetComponents(cs)
diff --git a/Processes/TrackingLine/testTrackingLineStack.cc b/Processes/TrackingLine/testTrackingLineStack.cc
deleted file mode 100644
index 4c2192388..000000000
--- a/Processes/TrackingLine/testTrackingLineStack.cc
+++ /dev/null
@@ -1,37 +0,0 @@
-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
-
-#ifndef _include_process_trackinling_teststack_h_
-#define _include_process_trackinling_teststack_h_
-
-typedef corsika::units::si::hepmomentum_d MOMENTUM;
-
-struct DummyParticle {
- HEPEnergyType fEnergy;
- Vector<MOMENTUM> fMomentum;
- Point fPosition;
-
- DummyParticle(HEPEnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
- : fEnergy(pEnergy)
- , fMomentum(pMomentum)
- , fPosition(pPosition) {}
-
- auto GetEnergy() const { return fEnergy; }
- auto GetMomentum() const { return fMomentum; }
- auto GetPosition() const { return fPosition; }
- auto GetPID() const { return corsika::particles::Code::Unknown; }
-};
-
-struct DummyStack {
- using ParticleType = DummyParticle;
-};
-
-#endif
--
GitLab