some progress...

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});
     }
   }
 

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());

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") {

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();
 

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

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.

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.

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>(

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>

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
-

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)>);

Processes/TrackingLine/testTrackingLine.cc

-
 /*
  * (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,

Processes/TrackingLine/testTrackingLineStack.h

-
 /*
  * (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