From c2a13032295b8dac49d0ca6162a8d7f38b7d1f3f Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 11 Nov 2020 19:21:03 +0100
Subject: [PATCH] migrated setup

---
 corsika/detail/setup/SetupEnvironment.inl |  40 ++++
 corsika/detail/setup/SetupStack.hpp       |  48 +++++
 corsika/detail/setup/SetupStack.inl       |  41 ++++
 corsika/setup/SetupEnvironment.hpp        |  43 ++++-
 corsika/setup/SetupStack.hpp              | 217 ++++++----------------
 corsika/setup/SetupTrajectory.hpp         |  49 -----
 6 files changed, 225 insertions(+), 213 deletions(-)
 create mode 100644 corsika/detail/setup/SetupEnvironment.inl
 create mode 100644 corsika/detail/setup/SetupStack.hpp
 create mode 100644 corsika/detail/setup/SetupStack.inl

diff --git a/corsika/detail/setup/SetupEnvironment.inl b/corsika/detail/setup/SetupEnvironment.inl
new file mode 100644
index 000000000..7b16b1c67
--- /dev/null
+++ b/corsika/detail/setup/SetupEnvironment.inl
@@ -0,0 +1,40 @@
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+
+#include <limits>
+
+namespace corsika::setup::testing {
+
+  inline std::tuple<std::unique_ptr<setup::Environment>, CoordinateSystem const*,
+                    setup::Environment::BaseNodeType const*>
+  setup_environment(Code vTargetCode) {
+
+    auto env = std::make_unique<setup::Environment>();
+    auto& universe = *(env->GetUniverse());
+    CoordinateSystem const& cs = env->GetCoordinateSystem();
+
+    /**
+     * our world is a sphere at 0,0,0 with R=infty
+     */
+    auto world = setup::Environment::CreateNode<Sphere>(
+        Point{cs, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+    /**
+     * construct suited environment medium model:
+     */
+    using MyHomogeneousModel = MediumPropertyModel<
+        UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+    world->SetModelProperties<MyHomogeneousModel>(
+        Medium::AirDry1Atm, Vector(cs, 0_T, 0_T, 1_T), 1_kg / (1_m * 1_m * 1_m),
+        NuclearComposition(std::vector<Code>{vTargetCode}, std::vector<float>{1.}));
+
+    setup::Environment::BaseNodeType const* nodePtr = world.get();
+    universe.AddChild(std::move(world));
+
+    return std::make_tuple(std::move(env), &cs, nodePtr);
+  }
+
+} // namespace corsika::setup::testing
diff --git a/corsika/detail/setup/SetupStack.hpp b/corsika/detail/setup/SetupStack.hpp
new file mode 100644
index 000000000..712556461
--- /dev/null
+++ b/corsika/detail/setup/SetupStack.hpp
@@ -0,0 +1,48 @@
+#pragma once
+
+#include <corsika/framework/stack/CombinedStack.hpp>
+#include <corsika/framework/stack/node/GeometryNodeStackExtension.hpp>
+#include <corsika/framework/stack/nuclear_extension/NuclearStackExtension.hpp>
+#include <corsika/framework/stack/history/HistorySecondaryProducer.hpp>
+#include <corsika/framework/stack/history/HistoryStackExtension.hpp>
+
+#include <corsika/setup/SetupEnvironment.hpp>
+
+namespace corsika {
+
+  namespace setup::detail {
+
+    // ------------------------------------------
+    // add geometry node tracking data to stack:
+
+    // the GeometryNode stack needs to know the type of geometry-nodes from the
+    // environment:
+    template <typename TStackIter>
+    using SetupGeometryDataInterface =
+        typename MakeGeometryDataInterface<TStackIter, setup::Environment>::type;
+
+    // combine particle data stack with geometry information for tracking
+    template <typename TStackIter>
+    using StackWithGeometryInterface =
+        CombinedParticleInterface<nuclear_extension::ParticleDataStack::MPIType,
+                                  SetupGeometryDataInterface, TStackIter>;
+
+    using StackWithGeometry =
+        CombinedStack<typename nuclear_extension::ParticleDataStack::StackImpl,
+                      GeometryData<setup::Environment>, StackWithGeometryInterface>;
+
+    // ------------------------------------------
+    // Add [optional] history data to stack, too:
+
+    // combine dummy stack with geometry information for tracking
+    template <typename TStackIter>
+    using StackWithHistoryInterface =
+        CombinedParticleInterface<StackWithGeometry::MPIType, HistoryEventDataInterface,
+                                  TStackIter>;
+
+    using StackWithHistory = CombinedStack<typename StackWithGeometry::StackImpl,
+                                           HistoryEventData, StackWithHistoryInterface>;
+
+  } // namespace setup::detail
+
+} // namespace corsika
diff --git a/corsika/detail/setup/SetupStack.inl b/corsika/detail/setup/SetupStack.inl
new file mode 100644
index 000000000..3c1133555
--- /dev/null
+++ b/corsika/detail/setup/SetupStack.inl
@@ -0,0 +1,41 @@
+#pragma once
+
+/**
+ * \function setup_stack
+ *
+ * standard stack setup for unit tests.
+ * \todo This can be moved to "test" directory, when available.
+ */
+
+namespace corsika::setup::testing {
+
+  inline std::tuple<std::unique_ptr<setup::Stack>, std::unique_ptr<setup::StackView>>
+  setup_stack(Code vProjectileType, int vA, int vZ, HEPEnergyType vMomentum,
+              setup::Environment::BaseNodeType* const vNodePtr,
+              CoordinateSystem const& cs) {
+
+    auto stack = std::make_unique<setup::Stack>();
+
+    Point const origin(cs, {0_m, 0_m, 0_m});
+    MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+
+    if (vProjectileType == Code::Nucleus) {
+      auto constexpr mN = constants::nucleonMass;
+      HEPEnergyType const E0 = sqrt(static_pow<2>(mN * vA) + pLab.squaredNorm());
+      auto particle = stack->AddParticle(
+          std::make_tuple(Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ));
+      particle.SetNode(vNodePtr);
+      return std::make_tuple(std::move(stack),
+                             std::make_unique<setup::StackView>(particle));
+    } else { // not a nucleus
+      HEPEnergyType const E0 =
+          sqrt(static_pow<2>(GetMass(vProjectileType)) + pLab.squaredNorm());
+      auto particle =
+          stack->AddParticle(std::make_tuple(vProjectileType, E0, pLab, origin, 0_ns));
+      particle.SetNode(vNodePtr);
+      return std::make_tuple(std::move(stack),
+                             std::make_unique<setup::StackView>(particle));
+    }
+  }
+
+} // namespace corsika::setup::testing
diff --git a/corsika/setup/SetupEnvironment.hpp b/corsika/setup/SetupEnvironment.hpp
index be970d2b1..3d88b8c18 100644
--- a/corsika/setup/SetupEnvironment.hpp
+++ b/corsika/setup/SetupEnvironment.hpp
@@ -9,10 +9,45 @@
 #pragma once
 
 #include <corsika/media/Environment.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
 #include <corsika/media/IMediumModel.hpp>
-#include <corsika/media/NameModel.hpp>
+#include <corsika/media/IMediumPropertyModel.hpp>
+#include <corsika/media/IRefractiveIndexModel.hpp>
 
 namespace corsika::setup {
-  using IEnvironmentModel = corsika::IMediumModel;
-  using SetupEnvironment = corsika::Environment<IEnvironmentModel>;
-} // namespace corsika::setup
+
+  /**
+     Definition of the default environemnt model interface. Each model
+     interface provides properties of the environment in a position
+     bdependent way.
+   */
+
+  using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+  using Environment = Environment<EnvironmentInterface>;
+
+} // end namespace corsika::setup
+
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/InhomogeneousMedium.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+
+#include <tuple>
+#include <unique_ptr>
+
+/**
+ * \function setup_environment
+ *
+ * standard environment for unit testing.
+ *
+ * \todo This can be moved to "test" directory, when available.
+ */
+namespace corsika::setup::testing {
+
+  inline std::tuple<std::unique_ptr<setup::Environment>, CoordinateSystem const*,
+                    setup::Environment::BaseNodeType const*>
+  setup_environment(Code vTargetCode);
+
+} // namespace corsika::setup::testing
+
+#include <corsika/detail/setup/SetupEnvironment.inl>
diff --git a/corsika/setup/SetupStack.hpp b/corsika/setup/SetupStack.hpp
index 8bb07507a..580e47fe2 100644
--- a/corsika/setup/SetupStack.hpp
+++ b/corsika/setup/SetupStack.hpp
@@ -6,131 +6,15 @@
  * the license.
  */
 
-#pragma once
+#include <corsika/detail/setup/SetupStack.hpp>
 
-// the basic particle data stack:
-#include <corsika/stack/SuperStupidStack.hpp>
-
-// extension with nuclear data for Code::Nucleus
-#include <corsika/stack/NuclearStackExtension.hpp>
-
-// extension with geometry information for tracking
-#include <corsika/framework/stack/CombinedStack.hpp>
-#include <corsika/media/Environment.hpp>
-
-#include <corsika/setup/SetupEnvironment.hpp>
-
-#include <tuple>
-#include <utility>
-#include <vector>
-
-  namespace detail {
-
-/**
- * @class GeometryData
- *
- * definition of stack-data object to store geometry information
- */
-class GeometryData {
-
-public:
-  using BaseNodeType = typename TEnvType::BaseNodeType;
-
-  // these functions are needed for the Stack interface
-  void Init() {}
-  void Clear() { fNode.clear(); }
-  unsigned int GetSize() const { return fNode.size(); }
-  unsigned int GetCapacity() const { return fNode.size(); }
-  void Copy(const int i1, const int i2) { fNode[i2] = fNode[i1]; }
-  void Swap(const int i1, const int i2) { std::swap(fNode[i1], fNode[i2]); }
-
-  // custom data access function
-  void SetNode(const int i, BaseNodeType const* v) { fNode[i] = v; }
-  auto const* GetNode(const int i) const { return fNode[i]; }
-
-  // these functions are also needed by the Stack interface
-  void IncrementSize() { fNode.push_back(nullptr); }
-  void DecrementSize() {
-    if (fNode.size() > 0) { fNode.pop_back(); }
-  }
-
-  // custom private data section
-private:
-  std::vector<const BaseNodeType*> fNode;
-};
-
-/**
- * @class GeometryDataInterface
- *
- * corresponding defintion of a stack-readout object, the iteractor
- * dereference operator will deliver access to these function
-// defintion of a stack-readout object, the iteractor dereference
-// operator will deliver access to these function
- */
-template <typename T, typename TEnvType>
-class GeometryDataInterface : public T {
-
-public:
-  using T::GetIndex;
-  using T::GetStackData;
-  using T::SetParticleData;
-  using BaseNodeType = typename TEnvType::BaseNodeType;
-
-  // default version for particle-creation from input data
-  void SetParticleData(const std::tuple<BaseNodeType const*> v) {
-    SetNode(std::get<0>(v));
-  }
-  void SetParticleData(GeometryDataInterface& parent,
-                       const std::tuple<BaseNodeType const*>) {
-    SetNode(parent.GetNode()); // copy Node from parent particle!
-  }
-  void SetParticleData() { SetNode(nullptr); }
-  void SetParticleData(GeometryDataInterface& parent) {
-    SetNode(parent.GetNode()); // copy Node from parent particle!
-  }
-  void SetNode(BaseNodeType const* v) { GetStackData().SetNode(GetIndex(), v); }
-  auto const* GetNode() const { return GetStackData().GetNode(GetIndex()); }
-};
+#include <array>
+#include <unique_ptr>
 
 namespace corsika::setup {
 
-    // the GeometryNode stack needs to know the type of geometry-nodes from the
-    // environment:
-    template <typename TStackIter>
-    using SetupGeometryDataInterface = typename stack::node::MakeGeometryDataInterface<
-        TStackIter, corsika::setup::Environment>::type;
-
-    //
-    // this is an auxiliary help typedef, which I don't know how to put
-    // into NuclearStackExtension.h where it belongs...
-    template <typename StackIter>
-    using ExtendedParticleInterfaceType =
-        corsika::nuclear_extension::NuclearParticleInterface<
-            corsika::super_stupid::SuperStupidStack::PIType, StackIter>;
-    //
-
-    // the particle data stack with extra nuclear information:
-    using ParticleDataStack = corsika::nuclear_extension::NuclearStackExtension<
-        corsika::super_stupid::SuperStupidStack, ExtendedParticleInterfaceType>;
-
-    // ------------------------------------------
-    // Add [optional] history data to stack, too:
-
-    // combine particle data stack with geometry information for tracking
-    template <typename StackIter>
-    using StackWithGeometryInterface =
-        corsika::CombinedParticleInterface<ParticleDataStack::PIType,
-                                           SetupGeometryDataInterface, StackIter>;
-
-    using StackWithGeometry =
-        corsika::CombinedStack<typename ParticleDataStack::StackImpl,
-                               GeometryData<setup::SetupEnvironment>,
-                               StackWithGeometryInterface>;
-
-  } // namespace detail
-
   // ---------------------------------------
-  // this is the FINAL stack we use in C8:
+  // this is the stack we use in C8 executables:
 
 #ifdef WITH_HISTORY
 
@@ -139,30 +23,7 @@ namespace corsika::setup {
    */
   using Stack = detail::StackWithHistory;
   template <typename T1, template <typename> typename M2>
-  using StackViewProducer = corsika::history::HistorySecondaryProducer<T1, M2>;
-
-  namespace detail {
-    /*
-      See Issue 161
-
-      unfortunately clang does not support this in the same way (yet) as
-      gcc, so we have to distinguish here. If clang cataches up, we
-      could remove the clang branch here and also in
-      corsika::Cascade. The gcc code is much more generic and
-      universal. If we could do the gcc version, we won't had to define
-      StackView globally, we could do it with MakeView whereever it is
-      actually needed. Keep an eye on this!
-    */
-#if defined(__clang__)
-    using TheStackView = corsika::stack::SecondaryView<
-        typename corsika::setup::Stack::StackImpl,
-        // CHECK with CLANG: corsika::setup::Stack::MPIType>;
-        corsika::setup::detail::StackWithHistoryInterface, StackViewProducer>;
-#elif defined(__GNUC__) || defined(__GNUG__)
-    using TheStackView =
-        corsika::stack::MakeView<corsika::setup::Stack, StackViewProducer>::type;
-#endif
-  } // namespace detail
+  using StackViewProducer = HistorySecondaryProducer<T1, M2>;
 
 #else // WITH_HISTORY
 
@@ -173,24 +34,42 @@ namespace corsika::setup {
   template <typename T1, template <typename> typename M2>
   using StackViewProducer = corsika::stack::DefaultSecondaryProducer<T1, M2>;
 
-  namespace detail {
-    /*
-      See Issue 161
-
-      unfortunately clang does not support this in the same way (yet) as
-      gcc, so we have to distinguish here. If clang cataches up, we
-      could remove the clang branch here and also in
-      corsika::Cascade. The gcc code is much more generic and
-      universal. If we could do the gcc version, we won't had to define
-      StackView globally, we could do it with MakeView whereever it is
-      actually needed. Keep an eye on this!
-    */
+#endif
+
+  // ---------------------------------------
+  // this is the stackitertor (particle type) we use in C8 executables:
+
+  /*
+    See Issue 161
+
+    unfortunately clang does not support this in the same way (yet) as
+    gcc, so we have to distinguish here. If clang cataches up, we
+    could remove the clang branch here and also in
+    corsika::Cascade. The gcc code is much more generic and
+    universal. If we could do the gcc version, we won't had to define
+    StackView globally, we could do it with MakeView whereever it is
+    actually needed. Keep an eye on this!
+  */
+
+#ifdef WITH_HISTORY
+
 #if defined(__clang__)
-  using StackView =
-      corsika::SecondaryView<typename corsika::setup::Stack::StackImpl,
-                             corsika::setup::detail::StackWithGeometryInterface>;
+  using StackView = SecondaryView<typename Stack::StackImpl,
+                                  // CHECK with CLANG: setup::Stack::MPIType>;
+                                  detail::StackWithHistoryInterface, StackViewProducer>;
 #elif defined(__GNUC__) || defined(__GNUG__)
-  using StackView = corsika::MakeView<corsika::setup::Stack>::type;
+  using StackView = make_view<setup::Stack, StackViewProducer>::type;
+#endif
+
+#else // WITH_HISTORY
+
+#if defined(__clang__)
+  using StackView = SecondaryView<typename setup::Stack::StackImpl,
+                                  // CHECK with CLANG:
+                                  // setup::Stack::MPIType>;
+                                  setup::detail::StackWithGeometryInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+  using StackView = make_view<setup::Stack>::type;
 #endif
   } // namespace detail
 
@@ -201,4 +80,22 @@ namespace corsika::setup {
 
   using StackView = detail::TheStackView;
 
+#endif // WITH_HISTORY
+
 } // namespace corsika::setup
+
+/**
+ * standard stack setup for unit tests. This can be moved to "test"
+ * directory, when available.
+ */
+
+namespace corsika::setup::testing {
+
+  inline std::tuple<std::unique_ptr<setup::Stack>, std::unique_ptr<setup::StackView>>
+  setup_stack(Code vProjectileType, int vA, int vZ, HEPEnergyType vMomentum,
+	      const setup::Environment::BaseNodeType* vNodePtr,
+	      CoordinateSystem const& cs);
+
+} // namespace corsika::setup::testing
+
+#include <corsika/detail/setup/SetupStack.inl>
diff --git a/corsika/setup/SetupTrajectory.hpp b/corsika/setup/SetupTrajectory.hpp
index 38b055a1b..8bed53ce8 100644
--- a/corsika/setup/SetupTrajectory.hpp
+++ b/corsika/setup/SetupTrajectory.hpp
@@ -8,61 +8,12 @@ n/*
 
 #pragma once
 
-#include <corsika/framework/geometry/Helix.hpp>
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/geometry/Trajectory.hpp>
 
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-#include <corsika/units/PhysicalUnits.h>
-
 namespace corsika::setup {
 
   /// definition of Trajectory base class, to be used in tracking and cascades
   typedef corsika::Trajectory<corsika::Line> Trajectory;
 
-  /*
-  typedef std::variant<std::monostate, corsika::Trajectory<Line>,
-                       corsika::Trajectory<Helix>>
-                       Trajectory;
-
-    tracking_leapfrog_straight::Tracking is a more simple and direct
-    leap-frog implementation. The two halve steps are coded explicitly
-    as two straight segments. Intersections with other volumes are
-    calculate only on the straight segments. This algorithm is based
-    on LineTrajectory.
-
-    tracking_line::Tracking is a pure straight tracker. It is based on
-    LineTrajectory.
-   */  
-  typedef corsika::process::tracking_leapfrog_curved::Tracking Tracking;
-  //typedef corsika::process::tracking_leapfrog_straight::Tracking Tracking;
-  //typedef corsika::process::tracking_line::Tracking Tracking;
-
-  /// definition of Trajectory base class, to be used in tracking and cascades
-  //typedef corsika::geometry::LineTrajectory Trajectory;
-  typedef corsika::geometry::LeapFrogTrajectory Trajectory;
-
-  /**
-
-     The following section is for unit testing only. Eventually it should
-     be moved to "tests".
-    
-    
-   */
-  
-  namespace testing {
-
-  /// helper visitor to modify Particle by moving along Trajectory
-  class GetDuration {
-  public:
-    TimeType operator()(std::monostate const&) {
-      return 0 * second;
-    }
-    template <typename T>
-    TimeType operator()(T const& trajectory) {
-      return trajectory.GetDuration();
-    }
-  };
-  */
 } // namespace corsika::setup
-- 
GitLab