From cce55ccc70f4422d7d6d0a002b325a1d0169364f Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Fri, 22 Mar 2019 00:10:21 -0300
Subject: [PATCH] some progress...
---
Documentation/Examples/cascade_example.cc | 83 ++++++++++---------
Framework/Cascade/testCascade.cc | 1 -
Framework/Geometry/testGeometry.cc | 6 +-
.../HadronicElasticModel.cc | 12 ++-
.../HadronicElasticModel.h | 4 -
Processes/Sibyll/Interaction.cc | 4 +-
Processes/Sibyll/NuclearInteraction.cc | 3 +-
Processes/Sibyll/testSibyll.cc | 7 +-
Processes/StackInspector/StackInspector.cc | 2 +-
Processes/TrackingLine/TrackingLine.cc | 5 +-
Processes/TrackingLine/TrackingLine.h | 2 +-
Processes/TrackingLine/testTrackingLine.cc | 29 ++++---
.../TrackingLine/testTrackingLineStack.h | 36 +++-----
13 files changed, 91 insertions(+), 103 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 9ff2e73cd..8b1f90be6 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -15,13 +15,13 @@
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
-#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/environment/NameModel.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NameModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h>
@@ -216,25 +216,29 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
-struct MyBoundaryCrossingProcess : public BoundaryCrossingProcess<MyBoundaryCrossingProcess> {
- //~ environment::BaseNodeType const& fA, fB;
-
- MyBoundaryCrossingProcess() {}
-
- //~ MyBoundaryCrossingProcess(environment::BaseNodeType const& a, environment::BaseNodeType const& b) : fA(a), fB(b) {}
-
- template <typename Particle>
- EProcessReturn DoBoundaryCrossing(Particle& p, environment::BaseNodeType const& from, environment::BaseNodeType const& to) {
- std::cout << "boundary crossing! from:" << &from << "; to: " << &to << std::endl;
-
- //~ if ((&fA == &from && &fB == &to) || (&fA == &to && &fB == &from)) {
- p.Delete();
- //~ }
-
- return EProcessReturn::eOk;
- }
-
- void Init() {}
+struct MyBoundaryCrossingProcess
+ : public BoundaryCrossingProcess<MyBoundaryCrossingProcess> {
+ //~ environment::BaseNodeType const& fA, fB;
+
+ MyBoundaryCrossingProcess() {}
+
+ //~ MyBoundaryCrossingProcess(environment::BaseNodeType const& a,
+ //environment::BaseNodeType const& b) : fA(a), fB(b) {}
+
+ 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;
+
+ //~ if ((&fA == &from && &fB == &to) || (&fA == &to && &fB == &from)) {
+ p.Delete();
+ //~ }
+
+ return EProcessReturn::eOk;
+ }
+
+ void Init() {}
};
//
@@ -246,34 +250,34 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- environment::Environment env;
+ using EnvType = environment::Environment<setup::IEnvironmentModel>;
+ EnvType env;
auto& universe = *(env.GetUniverse());
- auto outerMedium = environment::Environment::CreateNode<Sphere>(
+ auto outerMedium = environment::Environment<EnvType>::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
// fraction of oxygen
const float fox = 0.20946;
using MyHomogeneousModel = environment::HomogeneousMedium<setup::IEnvironmentModel>;
- outerMedium->SetModelProperties<setup::IEnvironmentModel>("outer",
- 1_kg / (1_m * 1_m * 1_m),
+ outerMedium->SetModelProperties<setup::IEnvironmentModel>(
+ "outer", 1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{(float)1. - fox, fox}));
-
- auto innerMedium = environment::Environment::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 2000_m);
- innerMedium->SetModelProperties<setup::IEnvironmentModel>("inner",
- 1_kg / (1_m * 1_m * 1_m),
+ auto innerMedium = environment::Environment<EnvType>::CreateNode<Sphere>(
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 2000_m);
+
+ innerMedium->SetModelProperties<setup::IEnvironmentModel>(
+ "inner", 1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{(float)1. - fox, fox}));
-
+
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
@@ -291,14 +295,16 @@ int main() {
ProcessCut cut(20_GeV);
random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
- process::HadronicElasticModel::HadronicElasticInteraction
- hadronicElastic(env);
+ process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic(env);
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
- auto sequence = p0 /*<< sibyll << sibyllNuc << decay << cut*/ << hadronicElastic << MyBoundaryCrossingProcess() << trackWriter;
+ auto sequence =
+ p0 /*<< sibyll << sibyllNuc << decay << cut*/ << hadronicElastic
+ << MyBoundaryCrossingProcess()
+ << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -329,10 +335,9 @@ int main() {
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
for (int l = 0; l < 100; ++l) {
- stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{
- beamCode, E0, plab, pos, 0_ns});
+ stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point,
+ units::si::TimeType>{beamCode, E0, plab, pos, 0_ns});
}
}
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 7d473eecc..5ec522a67 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -45,7 +45,6 @@ using namespace corsika::geometry;
#include <iostream>
using namespace std;
-
auto MakeDummyEnv() {
TestEnvironmentType env; // dummy environment
auto& universe = *(env.GetUniverse());
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index d73ae99b2..dae272be1 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -177,8 +177,10 @@ TEST_CASE("Trajectories") {
CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3));
-
- CHECK((base.NormalizedDirection().GetComponents(rootCS) - QuantityVector<dimensionless_d>{1, 0, 0}).eVector.norm() == Approx(0).margin(absMargin));
+
+ CHECK((base.NormalizedDirection().GetComponents(rootCS) -
+ QuantityVector<dimensionless_d>{1, 0, 0})
+ .eVector.norm() == Approx(0).margin(absMargin));
}
SECTION("Helix") {
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.cc b/Processes/HadronicElasticModel/HadronicElasticModel.cc
index 681d8fc0d..42a42e0ab 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.cc
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.cc
@@ -24,16 +24,15 @@
#include <iostream>
using namespace corsika::setup;
-using Particle = Stack::ParticleType;
-using Track = Trajectory;
+using Particle = corsika::setup::Stack::ParticleType;
+using Track = corsika::setup::Trajectory;
namespace corsika::process::HadronicElasticModel {
void HadronicElasticInteraction::Init() {}
- HadronicElasticInteraction::HadronicElasticInteraction(
- environment::Environment const& env, units::si::CrossSectionType x,
- units::si::CrossSectionType y)
+ HadronicElasticInteraction::HadronicElasticInteraction(units::si::CrossSectionType x,
+ units::si::CrossSectionType y)
: fX(x)
, fY(y) {}
@@ -87,8 +86,7 @@ namespace corsika::process::HadronicElasticModel {
using namespace units::si;
using namespace units::constants;
- const auto* currentNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto* currentNode = p.GetNode();
const auto& composition = currentNode->GetModelProperties().GetNuclearComposition();
const auto& components = composition.GetComponents();
diff --git a/Processes/HadronicElasticModel/HadronicElasticModel.h b/Processes/HadronicElasticModel/HadronicElasticModel.h
index 9760569ee..d0826a304 100644
--- a/Processes/HadronicElasticModel/HadronicElasticModel.h
+++ b/Processes/HadronicElasticModel/HadronicElasticModel.h
@@ -19,10 +19,6 @@
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
-namespace corsika::environment {
- class Environment;
-}
-
namespace corsika::process::HadronicElasticModel {
/**
* A simple model for elastic hadronic interactions based on the formulas
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 683730554..290c72b2a 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -271,8 +271,8 @@ namespace corsika::process::sibyll {
cross_section_of_components[i] = sigProd;
}
- const auto targetCode = mediumComposition.SampleTarget(
- cross_section_of_components, fRNG);
+ const auto targetCode =
+ mediumComposition.SampleTarget(cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 56a0c67f7..23ec7510d 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -353,7 +353,8 @@ namespace corsika::process::sibyll {
[[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS
}
- const auto targetCode = mediumComposition.SampleTarget(cross_section_of_components, fRNG);
+ const auto targetCode =
+ mediumComposition.SampleTarget(cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 6d286624e..318bb5c23 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -88,9 +88,10 @@ TEST_CASE("SibyllInterface", "[processes]") {
environment::Environment<environment::IMediumModel> env;
auto& universe = *(env.GetUniverse());
- auto theMedium = environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
- geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ auto theMedium =
+ environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 8552da191..a22742165 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -16,8 +16,8 @@
#include <corsika/logging/Logger.h>
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/cascade/testCascade.h>
+#include <corsika/setup/SetupTrajectory.h>
#include <iostream>
#include <limits>
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index 4806e7110..6b119c890 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -9,12 +9,12 @@
* the license.
*/
-#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/environment/Environment.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
#include <algorithm>
#include <iostream>
@@ -27,7 +27,7 @@ namespace corsika::process::tracking_line {
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
TimeOfIntersection(corsika::geometry::Line const& line,
- geometry::Sphere const& sphere) {
+ geometry::Sphere const& sphere) {
auto const delta = line.GetR0() - sphere.GetCenter();
auto const v = line.GetV0();
auto const vSqNorm = v.squaredNorm();
@@ -51,4 +51,3 @@ namespace corsika::process::tracking_line {
}
}
} // namespace corsika::process::tracking_line
-
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 55c865295..ef482fd3d 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -74,7 +74,7 @@ namespace corsika::process {
auto const& children = currentLogicalVolumeNode->GetChildNodes();
auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
- std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
+ std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
auto addIfIntersects = [&](auto const& vtn, auto const& nextNode) {
static_assert(std::is_same_v<decltype(vtn), decltype(nextNode)>);
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 19efb14f3..26dce27d1 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -1,4 +1,3 @@
-
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -36,10 +35,10 @@ using namespace std;
using namespace corsika::units::si;
TEST_CASE("TrackingLine") {
- corsika::environment::Environment env; // dummy environment
+ environment::Environment<environment::Empty> env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
- tracking_line::TrackingLine<DummyStack, setup::Trajectory> tracking(env);
+ tracking_line::TrackingLine tracking;
SECTION("intersection with sphere") {
Point const origin(cs, {0_m, 0_m, 0_m});
@@ -52,7 +51,7 @@ TEST_CASE("TrackingLine") {
setup::Trajectory traj(line, 12345_s);
auto const opt =
- tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
+ tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
REQUIRE(opt.has_value());
auto [t1, t2] = opt.value();
@@ -60,25 +59,29 @@ TEST_CASE("TrackingLine") {
REQUIRE(t2 / 11_s == Approx(1));
auto const optNoIntersection =
- tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
+ tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
REQUIRE_FALSE(optNoIntersection.has_value());
}
SECTION("maximally possible propagation") {
auto& universe = *(env.GetUniverse());
- //~ std::cout << env.GetUniverse().get() << std::endl;
-
auto const radius = 20_m;
- auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+ auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
universe.AddChild(std::move(theMedium));
-
- Point p0(cs, 0_m, 0_m, 0_m);
- auto* const node = universe.GetContainingNode(p0);
- DummyParticle p(1_GeV, Vector<MOMENTUM>(cs, 0_GeV, 0_GeV, 1_GeV),
- p0, node);
+
+ TestTrackingLineStack stack;
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::MuPlus,
+ 1_GeV,
+ {cs, {0_GeV, 0_GeV, 1_GeV}},
+ {cs, {0_m, 0_m, 0_km}},
+ 0_ns});
+ auto p = stack.GetNextParticle();
Point const origin(cs, {0_m, 0_m, 0_m});
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
diff --git a/Processes/TrackingLine/testTrackingLineStack.h b/Processes/TrackingLine/testTrackingLineStack.h
index 45975d032..5c3de074e 100644
--- a/Processes/TrackingLine/testTrackingLineStack.h
+++ b/Processes/TrackingLine/testTrackingLineStack.h
@@ -1,4 +1,3 @@
-
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -17,34 +16,19 @@
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/setup/SetupStack.h>
-typedef corsika::units::si::hepmomentum_d MOMENTUM;
-
-struct DummyParticle {
- corsika::units::si::HEPEnergyType fEnergy;
- corsika::geometry::Vector<MOMENTUM> fMomentum;
- corsika::geometry::Point fPosition;
- corsika::environment::VolumeTreeNode<> const* fNodePtr;
+using TestEnvironmentType = corsika::environment::Environment<corsika::environment::Empty>;
- DummyParticle(corsika::units::si::HEPEnergyType pEnergy,
- corsika::geometry::Vector<MOMENTUM> pMomentum,
- corsika::geometry::Point pPosition,
- corsika::environment::BaseNodeType const* pNodePtr)
- : fEnergy(pEnergy)
- , fMomentum(pMomentum)
- , fPosition(pPosition)
- , fNodePtr(pNodePtr) {}
+template <typename T>
+using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>;
- auto GetEnergy() const { return fEnergy; }
- auto GetMomentum() const { return fMomentum; }
- auto GetPosition() const { return fPosition; }
- auto GetPID() const { return corsika::particles::Code::Unknown; }
- auto* GetNode() const { return fNodePtr; }
-};
+// combine particle data stack with geometry information for tracking
+template <typename StackIter>
+using StackWithGeometryInterface =
+ corsika::stack::CombinedParticleInterface<corsika::setup::detail::ParticleDataStack::PIType,
+ SetupGeometryDataInterface, StackIter>;
+using TestTrackingLineStack = corsika::stack::CombinedStack<typename corsika::setup::detail::ParticleDataStack::StackImpl, GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
-struct DummyStack {
- using ParticleType = DummyParticle;
- using StackIterator = DummyParticle;
-};
#endif
--
GitLab