IAP GITLAB

Skip to content
Snippets Groups Projects
Commit c3863546 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

first try

parent 35455921
No related branches found
No related tags found
2 merge requests!91Resolve "define further classes of processes (MaintenanceProcess?)",!76Resolve "Handling of boundary crossings in geometry tree"
Pipeline #333 failed
Stage: build
Hide whitespace changes
Inline Side-by-side
Environment/VolumeTreeNode.h
+ 1
− 0
View file @ c3863546
......@@ -13,6 +13,7 @@
#define _include_VolumeTreeNode_H
#include <corsika/geometry/Volume.h>
#include <corsika/environment/IMediumModel.h>
#include <memory>
#include <vector>
......
Framework/Cascade/Cascade.h
+ 33
− 26
View file @ c3863546
......@@ -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);
}
}
......
Processes/TrackingLine/TrackingLine.cc
+ 0
− 67
View file @ c3863546
......@@ -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>
......
Processes/TrackingLine/TrackingLine.h
+ 95
− 2
View file @ c3863546
......@@ -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
......
Processes/TrackingLine/testTrackingLine.cc
+ 1
− 1
View file @ c3863546
......@@ -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)
......
Processes/TrackingLine/testTrackingLineStack.cc deleted 100644 → 0
+ 0
− 37
View file @ 35455921
/*
* (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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum