From 2e730f784d6d773310acf2e9ec865369a3ddb56f Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 26 Mar 2019 17:24:15 -0300
Subject: [PATCH] working prototype
---
Documentation/Examples/cascade_example.cc | 29 +++++++-----
Framework/Cascade/Cascade.h | 26 +++++++----
Processes/TrackingLine/TrackingLine.h | 55 ++++++++++++-----------
3 files changed, 64 insertions(+), 46 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index c0dcbdb7..615c2e1b 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -220,26 +220,28 @@ struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess> {
//~ environment::BaseNodeType const& fA, fB;
- MyBoundaryCrossingProcess(std::string const& filename) {fFile.open(filename);}
+ MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
template <typename Particle>
EProcessReturn DoBoundaryCrossing(Particle& p,
typename Particle::BaseNodeType const& from,
typename Particle::BaseNodeType const& to) {
- std::cout << "boundary crossing! from:" << &from << "; to: " << &to << std::endl;
+ std::cout << "boundary crossing! from: " << from.GetModelProperties().GetName()
+ << "; to: " << to.GetModelProperties().GetName() << std::endl;
//~ if ((&fA == &from && &fB == &to) || (&fA == &to && &fB == &from)) {
- //~ p.Delete();
+ p.Delete();
//~ }
-
+
auto const& name = particles::GetName(p.GetPID());
auto const start = p.GetPosition().GetCoordinates();
-
- fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' ' << start[2] / 1_m << '\n';
+
+ fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
+ << start[2] / 1_m << '\n';
return EProcessReturn::eOk;
}
-
+
std::ofstream fFile;
void Init() {}
@@ -286,7 +288,7 @@ int main() {
std::vector<float>{(float)1. - fox, fox}));
auto innerMedium = EnvType::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, -500_m}, 1000_m);
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, -500_m}, 5000_m);
innerMedium->SetModelProperties<
TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
@@ -299,6 +301,11 @@ int main() {
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
+ universe.SetModelProperties<
+ TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
+ "Universe", 0_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
@@ -322,8 +329,7 @@ int main() {
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence =
/*p0 <<*/ sibyll << sibyllNuc << decay << cut /* << hadronicElastic */
- << boundaryCrossing
- << trackWriter;
+ << boundaryCrossing << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -335,8 +341,7 @@ int main() {
(nuclA - nuclZ) * particles::Neutron::GetMass();
const HEPEnergyType E0 = nuclA * 10_TeV;
double phi = 0;
- for (double theta = 0; theta < 180; theta += 20)
- {
+ for (double theta = 0; theta < 180; theta += 20) {
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 45eaece1..a8185a0f 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -54,7 +54,8 @@ namespace corsika::cascade {
template <typename Tracking, typename ProcessList, typename Stack>
class Cascade {
using Particle = typename Stack::ParticleType;
- using VolumeTreeNode = std::remove_pointer_t<decltype(((Particle*) nullptr)->GetNode())>;
+ using VolumeTreeNode =
+ std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
using MediumInterface = typename VolumeTreeNode::IModelProperties;
// we only want fully configured objects
@@ -65,8 +66,8 @@ namespace corsika::cascade {
* Cascade class cannot be default constructed, but needs a valid
* list of physics processes for configuration at construct time.
*/
- Cascade(corsika::environment::Environment<MediumInterface> const& env, Tracking& tr, ProcessList& pl,
- Stack& stack)
+ Cascade(corsika::environment::Environment<MediumInterface> const& env, Tracking& tr,
+ ProcessList& pl, Stack& stack)
: fEnvironment(env)
, fTracking(tr)
, fProcessSequence(pl)
@@ -91,7 +92,8 @@ namespace corsika::cascade {
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
p.SetNode(numericalNode);
- std::cout << "initial node " << p.GetNode()->GetModelProperties().GetName() << std::endl;
+ std::cout << "initial node " << p.GetNode()->GetModelProperties().GetName()
+ << std::endl;
});
}
@@ -227,12 +229,18 @@ namespace corsika::cascade {
} else { // step-length limitation within volume
std::cout << "step-length limitation" << std::endl;
}
- std::cout << "nodes: " << currentLogicalNode->GetModelProperties().GetName() << " " << currentNumericalNode->GetModelProperties().GetName()
- << std::endl;
+ auto const* numericalNodeAfterStep =
+ fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+
+ std::cout << "nodes: " << currentLogicalNode->GetModelProperties().GetName()
+ << " " << numericalNodeAfterStep->GetModelProperties().GetName()
+ << std::endl;
- if (currentNumericalNode != currentLogicalNode) {
- throw std::runtime_error("numerical and logical nodes don't match");
- }
+ if (numericalNodeAfterStep != currentLogicalNode) {
+ std::cout << "position " << particle.GetPosition().GetCoordinates()
+ << std::endl;
+ throw std::runtime_error("numerical and logical nodes don't match");
+ }
} else { // boundary crossing
std::cout << "boundary crossing! next node = " << nextVol << std::endl;
particle.SetNode(nextVol);
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 00129700..4db45da2 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -12,35 +12,38 @@
#ifndef _include_corsika_processes_TrackingLine_h_
#define _include_corsika_processes_TrackingLine_h_
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/geometry/Trajectory.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
#include <optional>
#include <type_traits>
#include <utility>
namespace corsika::environment {
- template <typename IEnvironmentModel> class Environment;
- template <typename IEnvironmentModel> class VolumeTreeNode;
-}
+ template <typename IEnvironmentModel>
+ class Environment;
+ template <typename IEnvironmentModel>
+ class VolumeTreeNode;
+} // namespace corsika::environment
namespace corsika::process {
namespace tracking_line {
-
- std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
- TimeOfIntersection(corsika::geometry::Line const& line,
- geometry::Sphere const& sphere);
+
+ std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+ TimeOfIntersection(corsika::geometry::Line const& line,
+ geometry::Sphere const& sphere);
class TrackingLine {
public:
- TrackingLine() {};
+ TrackingLine(){};
- template <typename Particle> // was Stack previously, and argument was Stack::StackIterator
- auto GetTrack(Particle const& p) {
+ template <typename Particle> // was Stack previously, and argument was
+ // Stack::StackIterator
+ auto GetTrack(Particle const& p) {
using namespace corsika::units::si;
using namespace corsika::geometry;
geometry::Vector<SpeedType::dimension_type> const velocity =
@@ -49,8 +52,8 @@ namespace corsika::process {
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 pos: " << currentPosition.GetCoordinates() << " ["
+ << p.GetNode()->GetModelProperties().GetName() << "]\n";
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
@@ -64,15 +67,15 @@ namespace corsika::process {
//~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
auto const numericallyInside =
currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
-
- std::cout << "numericallyInside = " << (numericallyInside ? "true":"false");
+
+ std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
auto const& children = currentLogicalVolumeNode->GetChildNodes();
auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
- // for entering from outside
+ // for entering from outside
auto addIfIntersects = [&](auto const& vtn) {
auto const& volume = vtn.GetVolume();
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
@@ -81,12 +84,12 @@ namespace corsika::process {
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;
+ std::cout << "intersection times: " << t1 / 1_s << "; " << t2 / 1_s << " "
+ << vtn.GetModelProperties().GetName() << std::endl;
if (t1.magnitude() > 0)
intersections.emplace_back(t1, &vtn);
else if (t2.magnitude() > 0)
- throw std::runtime_error("inside other volume");
+ std::cout << "inside other volume" << std::endl;
}
};
@@ -116,11 +119,13 @@ namespace corsika::process {
min = minIter->first;
}
- std::cout << " t-intersect: " << min << " " << minIter->second->GetModelProperties().GetName() << std::endl;
- std::cout << "point of intersection: " << line.GetPosition(min).GetCoordinates() << std::endl;
+ std::cout << " t-intersect: " << min << " "
+ << minIter->second->GetModelProperties().GetName() << std::endl;
+ //~ std::cout << "point of intersection: " <<
+ //line.GetPosition(min).GetCoordinates() << std::endl;
- return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min), velocity.norm() * min,
- minIter->second);
+ return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
+ velocity.norm() * min, minIter->second);
}
};
--
GitLab