diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc
index 73063442169cf2b67f342edebadfe502cae06ae3..528c0d4436a2ccf00f226642fc5eda20ce8afd55 100644
--- a/Documentation/Examples/geometry_example.cc
+++ b/Documentation/Examples/geometry_example.cc
@@ -2,7 +2,7 @@
 #include <fwk/Sphere.h>
 #include <fwk/Point.h>
 #include <fwk/CoordinateSystem.h>
-#include <Units/PhysicalUnits.h>
+#include <fwk/PhysicalUnits.h>
 
 #include <iostream>
 #include <typeinfo>
diff --git a/Documentation/Examples/helix_example.cc b/Documentation/Examples/helix_example.cc
index 46b9f7161d864c00f13855f9b73ca8eed9a70ecf..c365c83b0fd300b3fac3fa23040031fc8883a7b0 100644
--- a/Documentation/Examples/helix_example.cc
+++ b/Documentation/Examples/helix_example.cc
@@ -1,46 +1,46 @@
-#include <Units/PhysicalUnits.h>
-#include <fwk/Vector.h>
 #include <fwk/CoordinateSystem.h>
-#include <fwk/Point.h>
 #include <fwk/Helix.h>
+#include <fwk/PhysicalUnits.h>
+#include <fwk/Point.h>
+#include <fwk/Vector.h>
+
+#include <array>
 #include <cstdlib>
 #include <iostream>
-#include <array>
 
-using namespace phys::units; 
-using namespace phys::units::io; // support stream << unit
+using namespace phys::units;
+using namespace phys::units::io;       // support stream << unit
 using namespace phys::units::literals; // support unit literals like 5_m;
 
-int main()
-{
-    CoordinateSystem root;
-    
-    fwk::Point const r0(root, {0_m, 0_m, 0_m});
-    auto const omegaC = 2 * M_PI * 1_Hz;
-    fwk::Vector<speed_d> vPar(root, {0_m / second, 0_m / second, 10_cm / second});
-    fwk::Vector<speed_d> vPerp(root, {1_m / second, 0_m / second, 0_m / second});
-    
-    fwk::Helix h(r0, omegaC, vPar, vPerp);
-    
-    auto constexpr t0 = 0_s;
-    auto constexpr t1 = 1_s;
-    auto constexpr dt = 1_us;
-    auto constexpr n = long((t1 - t0) / dt) + 1;
-    
-    auto arr = std::make_unique<std::array<std::array<double, 4>, n>>();
-    auto& positions = *arr;
-    
-    for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i)
-    {
-        auto const r = h.GetPosition(t).GetCoordinates();
-        
-        positions[i][0] = t / 1_s;
-        positions[i][1] = r[0] / 1_m;
-        positions[i][2] = r[1] / 1_m;
-        positions[i][3] = r[2] / 1_m;
-    }
-    
-    std::cout << positions[n-2][0] << " " << positions[n-2][1] << " " << positions[n-2][2] << " " << positions[n-2][3] << std::endl;
-    
-    return EXIT_SUCCESS;
+int main() {
+  fwk::CoordinateSystem root;
+
+  fwk::Point const r0(root, {0_m, 0_m, 0_m});
+  auto const omegaC = 2 * M_PI * 1_Hz;
+  fwk::Vector<speed_d> vPar(root, {0_m / second, 0_m / second, 10_cm / second});
+  fwk::Vector<speed_d> vPerp(root, {1_m / second, 0_m / second, 0_m / second});
+
+  fwk::Helix h(r0, omegaC, vPar, vPerp);
+
+  auto constexpr t0 = 0_s;
+  auto constexpr t1 = 1_s;
+  auto constexpr dt = 1_us;
+  auto constexpr n = long((t1 - t0) / dt) + 1;
+
+  auto arr = std::make_unique<std::array<std::array<double, 4>, n>>();
+  auto& positions = *arr;
+
+  for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i) {
+    auto const r = h.GetPosition(t).GetCoordinates();
+
+    positions[i][0] = t / 1_s;
+    positions[i][1] = r[0] / 1_m;
+    positions[i][2] = r[1] / 1_m;
+    positions[i][3] = r[2] / 1_m;
+  }
+
+  std::cout << positions[n - 2][0] << " " << positions[n - 2][1] << " "
+            << positions[n - 2][2] << " " << positions[n - 2][3] << std::endl;
+
+  return EXIT_SUCCESS;
 }
diff --git a/Documentation/Examples/stack_example.cc b/Documentation/Examples/stack_example.cc
index c526b5887a886c6061b517abad2eda69e8deb839..0529656635f311d6ad216f5a6854802940da3d7f 100644
--- a/Documentation/Examples/stack_example.cc
+++ b/Documentation/Examples/stack_example.cc
@@ -1,34 +1,29 @@
-#include <ParticleStack/StackOne.h>
+#include <fwk/StackOne.h>
 //#include <corsika/StackOne.h>
 
-#include <iostream>
 #include <iomanip>
+#include <iostream>
 
 using namespace std;
 
-void fill(stack::StackOne& s)
-{
-  for (int i=0; i<11; ++i) {
+//using namespace fwk;
+
+void fill(stack::StackOne& s) {
+  for (int i = 0; i < 11; ++i) {
     auto p = s.NewParticle();
     p.SetId(i);
-    p.SetEnergy(1.5*i);
+    p.SetEnergy(1.5 * i);
   }
 }
 
-void
-read(stack::StackOne& s)
-{
+void read(stack::StackOne& s) {
   cout << "found Stack with " << s.GetSize() << " particles. " << endl;
   double Etot = 0;
-  for (auto p : s) {
-    Etot += p.GetEnergy();
-  }
+  for (auto p : s) { Etot += p.GetEnergy(); }
   cout << "Etot=" << Etot << endl;
 }
 
-int
-main()
-{
+int main() {
   stack::StackOne s;
   fill(s);
   read(s);
diff --git a/Framework/Geometry/CoordinateSystem.cc b/Framework/Geometry/CoordinateSystem.cc
index 9c9f97e39f4088a5b2adcc2c1ac1ac24e69f33d1..2763ae2046d518c6c3e3a6ae0b4a2a6d595408ac 100644
--- a/Framework/Geometry/CoordinateSystem.cc
+++ b/Framework/Geometry/CoordinateSystem.cc
@@ -1,54 +1,45 @@
-#include <Geometry/CoordinateSystem.h>
-
-EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2)
-{
-    CoordinateSystem const* a{&c1};
-    CoordinateSystem const* b{&c2};
-    CoordinateSystem const* commonBase{nullptr};
-    
-    while (a != b && b != nullptr)
-    {
-        a = &c1;
-        
-        while (a != b && a != nullptr)
-        {
-            a = a->GetReference();
-        }
-        
-        if (a == b)
-            break;
-        
-        b = b->GetReference();
-    }
-    
-    if (a == b && a != nullptr)
-    {
-        commonBase = a;
-        
-    }
-    else
-    {
-        throw std::string("no connection between coordinate systems found!");
-    }
-    
-    
-    EigenTransform t = EigenTransform::Identity();
-    
-    auto* p = &c1;
-    
-    while (p != commonBase)
-    {
-        t = p->GetTransform() * t;
-        p = p->GetReference();
-    }
-    
-    p = &c2;
-    
-    while (p != commonBase)
-    {
-        t = p->GetTransform().inverse(Eigen::TransformTraits::Isometry) * t;
-        p = p->GetReference();
-    }
-    
-    return t;
+#include <fwk/CoordinateSystem.h>
+
+using namespace fwk;
+
+EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& c1,
+                                                   CoordinateSystem const& c2) {
+  CoordinateSystem const* a{&c1};
+  CoordinateSystem const* b{&c2};
+  CoordinateSystem const* commonBase{nullptr};
+
+  while (a != b && b != nullptr) {
+    a = &c1;
+
+    while (a != b && a != nullptr) { a = a->GetReference(); }
+
+    if (a == b) break;
+
+    b = b->GetReference();
+  }
+
+  if (a == b && a != nullptr) {
+    commonBase = a;
+
+  } else {
+    throw std::string("no connection between coordinate systems found!");
+  }
+
+  EigenTransform t = EigenTransform::Identity();
+
+  auto* p = &c1;
+
+  while (p != commonBase) {
+    t = p->GetTransform() * t;
+    p = p->GetReference();
+  }
+
+  p = &c2;
+
+  while (p != commonBase) {
+    t = p->GetTransform().inverse(Eigen::TransformTraits::Isometry) * t;
+    p = p->GetReference();
+  }
+
+  return t;
 }
diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index 238054ca7db7698f302efcd09cc90abc14a688b5..20348ffe73cbc320555ddf02c50b2e9099b30a31 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/Framework/Geometry/CoordinateSystem.h
@@ -1,72 +1,63 @@
 #ifndef _include_COORDINATESYSTEM_H_
 #define _include_COORDINATESYSTEM_H_
 
-#include <Geometry/QuantityVector.h>
-#include <Units/PhysicalUnits.h>
+#include <fwk/PhysicalUnits.h>
+#include <fwk/QuantityVector.h>
 
 #include <Eigen/Dense>
 
 typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
 typedef Eigen::Translation<double, 3> EigenTranslation;
 
-class CoordinateSystem
-{
+namespace fwk {
+
+  class CoordinateSystem {
     CoordinateSystem const* reference = nullptr;
     EigenTransform transf;
-    
-    CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf) :
-        reference(&reference),
-        transf(transf)
-    {
-        
-    }
-    
-public:
-    static EigenTransform GetTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2);
-
-    CoordinateSystem() : // for creating the root CS
-        transf(EigenTransform::Identity())
-    {
-        
-    }
-    
-    auto& operator=(const CoordinateSystem& pCS)
-    {
-        reference = pCS.reference;
-        transf = pCS.transf;
-        return *this;
-    }    
-    
-    auto translate(QuantityVector<phys::units::length_d> vector) const
-    {
-        EigenTransform const translation{EigenTranslation(vector.eVector)};
-        
-        return CoordinateSystem(*this, translation);
+
+    CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
+        : reference(&reference)
+        , transf(transf) {}
+
+  public:
+    static EigenTransform GetTransformation(CoordinateSystem const& c1,
+                                            CoordinateSystem const& c2);
+
+    CoordinateSystem()
+        : // for creating the root CS
+        transf(EigenTransform::Identity()) {}
+
+    auto& operator=(const CoordinateSystem& pCS) {
+      reference = pCS.reference;
+      transf = pCS.transf;
+      return *this;
     }
-    
-    auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const
-    {
-        EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
-        
-        return CoordinateSystem(*this, rotation);
+
+    auto translate(QuantityVector<phys::units::length_d> vector) const {
+      EigenTransform const translation{EigenTranslation(vector.eVector)};
+
+      return CoordinateSystem(*this, translation);
     }
-    
-    auto translateAndRotate(QuantityVector<phys::units::length_d> translation, QuantityVector<phys::units::length_d> axis, double angle)
-    {
-        EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) * EigenTranslation(translation.eVector)};
-        
-        return CoordinateSystem(*this, transf);
+
+    auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
+      EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
+
+      return CoordinateSystem(*this, rotation);
     }
-    
-    auto const* GetReference() const
-    {
-        return reference;
+
+    auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
+                            QuantityVector<phys::units::length_d> axis, double angle) {
+      EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
+                                  EigenTranslation(translation.eVector)};
+
+      return CoordinateSystem(*this, transf);
     }
-    
-    auto const& GetTransform() const
-    {
-        return transf;
-    }    
-};
+
+    auto const* GetReference() const { return reference; }
+
+    auto const& GetTransform() const { return transf; }
+  };
+
+} // namespace fwk
 
 #endif
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 1234940b3f16aa1e4345338ab9e36835f2255d89..38ba50c5cb4aa2beac1f59b1daafd9d19d2317f1 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -1,45 +1,44 @@
 #ifndef _include_HELIX_H_
 #define _include_HELIX_H_
 
-#include <fwk/Vector.h>
 #include <fwk/Point.h>
-#include <Units/PhysicalUnits.h>
+#include <fwk/Vector.h>
+
+#include <fwk/PhysicalUnits.h>
 
 #include <cmath>
 
 namespace fwk {
 
-class Helix // TODO: inherit from to-be-implemented "Trajectory"
-{
+  class Helix // TODO: inherit from to-be-implemented "Trajectory"
+  {
     using SpeedVec = Vector<Speed::dimension_type>;
-    
+
     Point const r0;
     Frequency const omegaC;
     SpeedVec const vPar;
     SpeedVec vPerp, uPerp;
-    
+
     Length const radius;
-    
-public:
+
+  public:
     Helix(Point const& pR0, phys::units::quantity<phys::units::frequency_d> pOmegaC,
-        SpeedVec const& pvPar, SpeedVec const& pvPerp) :
-        r0(pR0), omegaC(pOmegaC), vPar(pvPar), vPerp(pvPerp), uPerp(vPerp.cross(vPar.normalized())),
-        radius(pvPar.norm() / abs(pOmegaC))
-    {
-        
-    }
-    
-    auto GetPosition(Time t) const
-    {
-        return r0 + vPar * t + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
-    }
-    
-    auto GetRadius() const
-    {
-        return radius;
+          SpeedVec const& pvPar, SpeedVec const& pvPerp)
+        : r0(pR0)
+        , omegaC(pOmegaC)
+        , vPar(pvPar)
+        , vPerp(pvPerp)
+        , uPerp(vPerp.cross(vPar.normalized()))
+        , radius(pvPar.norm() / abs(pOmegaC)) {}
+
+    auto GetPosition(Time t) const {
+      return r0 + vPar * t +
+             (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
     }
-};
 
-} // end namesapce
- 
+    auto GetRadius() const { return radius; }
+  };
+
+} // namespace fwk
+
 #endif
diff --git a/Framework/Geometry/Point.h b/Framework/Geometry/Point.h
index 6607f1bdc71c9a58d79e8f435e0bc7c0b87fa034..47f5e954d0f50938148b1048063e4300ae594c29 100644
--- a/Framework/Geometry/Point.h
+++ b/Framework/Geometry/Point.h
@@ -4,8 +4,8 @@
 #include <fwk/BaseVector.h>
 #include <fwk/QuantityVector.h>
 #include <fwk/Vector.h>
-#include <Units/PhysicalUnits.h>
 
+#include <fwk/PhysicalUnits.h>
 
 namespace fwk {
 
@@ -13,61 +13,51 @@ namespace fwk {
    * A Point represents a point in position space. It is defined by its
    * coordinates with respect to some CoordinateSystem.
    */
-  class Point : public BaseVector<phys::units::length_d>
-  {
+  class Point : public BaseVector<phys::units::length_d> {
   public:
-  Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector) :
-    BaseVector<phys::units::length_d>(pCS, pQVector)
-      {
-      }
-    
-  Point(CoordinateSystem const& cs, Length x, Length y, Length z) :
-        BaseVector<phys::units::length_d>(cs, {x, y, z})
-    {
-    }
-    
-    auto GetCoordinates() const
-    {
+    Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector)
+        : BaseVector<phys::units::length_d>(pCS, pQVector) {}
+
+    Point(CoordinateSystem const& cs, Length x, Length y, Length z)
+        : BaseVector<phys::units::length_d>(cs, {x, y, z}) {}
+
+    auto GetCoordinates() const { return BaseVector<phys::units::length_d>::qVector; }
+
+    auto GetCoordinates(CoordinateSystem const& pCS) const {
+      if (&pCS == BaseVector<phys::units::length_d>::cs) {
         return BaseVector<phys::units::length_d>::qVector;
+      } else {
+        return QuantityVector<phys::units::length_d>(
+            CoordinateSystem::GetTransformation(*BaseVector<phys::units::length_d>::cs,
+                                                pCS) *
+            BaseVector<phys::units::length_d>::qVector.eVector);
+      }
     }
-    
-    auto GetCoordinates(CoordinateSystem const& pCS) const
-    {
-        if (&pCS == BaseVector<phys::units::length_d>::cs)
-        {
-            return BaseVector<phys::units::length_d>::qVector;
-        }
-        else
-        {
-            return QuantityVector<phys::units::length_d>(CoordinateSystem::GetTransformation(*BaseVector<phys::units::length_d>::cs, pCS) * BaseVector<phys::units::length_d>::qVector.eVector);
-        }
-    }
-    
+
     /*!
      * transforms the Point into another CoordinateSystem by changing its
      * coordinates interally
      */
-    void rebase(CoordinateSystem const& pCS)
-    {
-        BaseVector<phys::units::length_d>::qVector = GetCoordinates(pCS);
-        BaseVector<phys::units::length_d>::cs = &pCS;
+    void rebase(CoordinateSystem const& pCS) {
+      BaseVector<phys::units::length_d>::qVector = GetCoordinates(pCS);
+      BaseVector<phys::units::length_d>::cs = &pCS;
     }
-    
-    Point operator+(Vector<phys::units::length_d> const& pVec) const
-    {        
-        return Point(*BaseVector<phys::units::length_d>::cs, GetCoordinates() + pVec.GetComponents(*BaseVector<phys::units::length_d>::cs));
+
+    Point operator+(Vector<phys::units::length_d> const& pVec) const {
+      return Point(
+          *BaseVector<phys::units::length_d>::cs,
+          GetCoordinates() + pVec.GetComponents(*BaseVector<phys::units::length_d>::cs));
     }
-    
+
     /*!
      * returns the distance Vector between two points
      */
-    Vector<phys::units::length_d> operator-(Point const& pB) const
-    {
-        auto& cs = *BaseVector<phys::units::length_d>::cs;
-        return Vector<phys::units::length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
+    Vector<phys::units::length_d> operator-(Point const& pB) const {
+      auto& cs = *BaseVector<phys::units::length_d>::cs;
+      return Vector<phys::units::length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
     }
-};
+  };
+
+} // namespace fwk
 
-} // end namespace
-  
 #endif
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index 4f9209f8e2ee76bb6549f892712e7cf2cea563fc..0694f0d196a2f47237b5a2de10532978178c3d11 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -1,139 +1,118 @@
 #ifndef _include_QUANTITYVECTOR_H_
 #define _include_QUANTITYVECTOR_H_
 
-#include <Units/PhysicalUnits.h>
+#include <fwk/PhysicalUnits.h>
 
 #include <Eigen/Dense>
+
 #include <iostream>
 #include <utility>
 
-/*!
- * A QuantityVector is a three-component container based on Eigen::Vector3d
- * with a phys::units::dimension. Arithmethic operators are defined that
- * propagate the dimensions by dimensional analysis.
- */
+namespace fwk {
 
-template <typename dim>
-class QuantityVector
-{
-protected:
-    using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity corresponding to the dimension
-    
-public:
+  /*!
+   * A QuantityVector is a three-component container based on Eigen::Vector3d
+   * with a phys::units::dimension. Arithmethic operators are defined that
+   * propagate the dimensions by dimensional analysis.
+   */
+
+  template <typename dim>
+  class QuantityVector {
+  protected:
+    using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity
+                                                         // corresponding to the dimension
+
+  public:
     Eigen::Vector3d eVector; //!< the actual container where the raw numbers are stored
-    
+
     typedef dim dimension; //!< should be a phys::units::dimension
 
-    QuantityVector(Quantity a, Quantity b, Quantity c) :
-        eVector{a.magnitude(), b.magnitude(), c.magnitude()}
-    {
-    }
-    
-    QuantityVector(Eigen::Vector3d pBareVector) :
-        eVector(pBareVector)
-    {
-    }
-    
-    auto operator[](size_t index) const
-    {
-        return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
+    QuantityVector(Quantity a, Quantity b, Quantity c)
+        : eVector{a.magnitude(), b.magnitude(), c.magnitude()} {}
+
+    QuantityVector(Eigen::Vector3d pBareVector)
+        : eVector(pBareVector) {}
+
+    auto operator[](size_t index) const {
+      return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
     }
-    
-    auto norm() const
-    {
-        return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
+
+    auto norm() const {
+      return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
     }
-    
-    auto squaredNorm() const
-    {
-        using QuantitySquared = decltype(std::declval<Quantity>() * std::declval<Quantity>());
-        return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
+
+    auto squaredNorm() const {
+      using QuantitySquared =
+          decltype(std::declval<Quantity>() * std::declval<Quantity>());
+      return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
     }
-    
-    auto operator+(QuantityVector<dim> const& pQVec) const
-    {
-        return QuantityVector<dim>(eVector + pQVec.eVector);
+
+    auto operator+(QuantityVector<dim> const& pQVec) const {
+      return QuantityVector<dim>(eVector + pQVec.eVector);
     }
-    
-    auto operator-(QuantityVector<dim> const& pQVec) const
-    {
-        return QuantityVector<dim>(eVector - pQVec.eVector);
+
+    auto operator-(QuantityVector<dim> const& pQVec) const {
+      return QuantityVector<dim>(eVector - pQVec.eVector);
     }
-    
+
     template <typename ScalarDim>
-    auto operator*(phys::units::quantity<ScalarDim, double> const p) const
-    {
-        using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
-        
-        if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
-        {
-            return QuantityVector<phys::units::dimensionless_d>(eVector * p.magnitude());
-        }
-        else
-        {
-            return QuantityVector<typename ResQuantity::dimension_type>(eVector * p.magnitude());
-        }
+    auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
+      using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
+
+      if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
+                                                              // a "Quantity" anymore
+      {
+        return QuantityVector<phys::units::dimensionless_d>(eVector * p.magnitude());
+      } else {
+        return QuantityVector<typename ResQuantity::dimension_type>(eVector *
+                                                                    p.magnitude());
+      }
     }
-    
+
     template <typename ScalarDim>
-    auto operator/(phys::units::quantity<ScalarDim, double> const p) const
-    {
-        return (*this) * (1 / p);
-    }
-    
-    auto operator*(double const p) const
-    {
-        return QuantityVector<dim>(eVector * p);
-    }
-    
-    auto operator/(double const p) const
-    {
-        return QuantityVector<dim>(eVector / p);
-    }
-    
-    auto& operator/=(double const p)
-    {
-        eVector /= p;
-        return *this;
+    auto operator/(phys::units::quantity<ScalarDim, double> const p) const {
+      return (*this) * (1 / p);
     }
-    
-    auto& operator*=(double const p)
-    {
-        eVector *= p;
-        return *this;
-    }
-    
-    auto& operator+=(QuantityVector<dim> const& pQVec)
-    {
-        eVector += pQVec.eVector;
-        return *this;
+
+    auto operator*(double const p) const { return QuantityVector<dim>(eVector * p); }
+
+    auto operator/(double const p) const { return QuantityVector<dim>(eVector / p); }
+
+    auto& operator/=(double const p) {
+      eVector /= p;
+      return *this;
     }
-    
-    auto& operator-=(QuantityVector<dim> const& pQVec)
-    {
-        eVector -= pQVec.eVector;
-        return *this;
+
+    auto& operator*=(double const p) {
+      eVector *= p;
+      return *this;
     }
-    
-    auto& operator-() const
-    {
-        return QuantityVector<dim>(-eVector);
+
+    auto& operator+=(QuantityVector<dim> const& pQVec) {
+      eVector += pQVec.eVector;
+      return *this;
     }
-    
-    auto normalized() const
-    {
-        return (*this) * (1 / norm());
+
+    auto& operator-=(QuantityVector<dim> const& pQVec) {
+      eVector -= pQVec.eVector;
+      return *this;
     }
-};
+
+    auto& operator-() const { return QuantityVector<dim>(-eVector); }
+
+    auto normalized() const { return (*this) * (1 / norm()); }
+  };
+
+} // end namespace fwk
 
 template <typename dim>
-auto& operator<<(std::ostream& os, QuantityVector<dim> qv)
-{
-    using Quantity = phys::units::quantity<dim, double>;
-    
-    os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2)
-       << ") " << phys::units::to_unit_symbol<dim, double>(Quantity(phys::units::detail::magnitude_tag, 1));
-    return os;
+auto& operator<<(std::ostream& os, fwk::QuantityVector<dim> qv) {
+  using Quantity = phys::units::quantity<dim, double>;
+
+  os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2) << ") "
+     << phys::units::to_unit_symbol<dim, double>(
+            Quantity(phys::units::detail::magnitude_tag, 1));
+  return os;
 }
 
 #endif
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 634ffe31ec1599565110d297cc1ad2bc6910e3e1..39d3dcb426f3b28ce95f377a6cb99022c6314bb7 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -1,30 +1,25 @@
 #ifndef _include_SPHERE_H_
 #define _include_SPHERE_H_
 
+#include <fwk/PhysicalUnits.h>
 #include <fwk/Point.h>
-#include <Units/PhysicalUnits.h>
 
 namespace fwk {
 
-class Sphere
-{
+  class Sphere {
     Point center;
     Length const radius;
-    
-public:
-    Sphere(Point const& pCenter, Length const pRadius) :
-        center(pCenter), radius(pRadius)
-    {
-        
-    }
-        
-    auto isInside(Point const& p) const
-    {
-        return radius*radius > (center - p).squaredNorm();
-    }
 
-};
+  public:
+    Sphere(Point const& pCenter, Length const pRadius)
+        : center(pCenter)
+        , radius(pRadius) {}
+
+    auto isInside(Point const& p) const {
+      return radius * radius > (center - p).squaredNorm();
+    }
+  };
 
-}// end namespace
+} // namespace fwk
 
 #endif
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index 113d8701813d46e9023a7b0cd127501c9e282c73..f7191230ca844350c1ed46e142a1e94cbd2c7595 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -3,14 +3,15 @@
 
 #include <fwk/BaseVector.h>
 #include <fwk/QuantityVector.h>
-#include <Units/PhysicalUnits.h>
+
+#include <fwk/PhysicalUnits.h>
 
 /*!
  * A Vector represents a 3-vector in Euclidean space. It is defined by components
  * given in a specific CoordinateSystem. It has a physical dimension ("unit")
  * as part of its type, so you cannot mix up e.g. electric with magnetic fields
  * (but you could calculate their cross-product to get an energy flux vector).
- * 
+ *
  * When transforming coordinate systems, a Vector is subject to the rotational
  * part only and invariant under translations.
  */
@@ -18,74 +19,57 @@
 namespace fwk {
 
   template <typename dim>
-    class Vector : public BaseVector<dim>
-{
+  class Vector : public BaseVector<dim> {
     using Quantity = phys::units::quantity<dim, double>;
-    
-public:
-    Vector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector) :
-        BaseVector<dim>(pCS, pQVector)
-    {
-    }
 
-    Vector(CoordinateSystem const& cs, Quantity x, Quantity y, Quantity z) :
-        BaseVector<dim>(cs, QuantityVector<dim>(x, y, z))
-    {
-    }
-    
+  public:
+    Vector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
+        : BaseVector<dim>(pCS, pQVector) {}
+
+    Vector(CoordinateSystem const& cs, Quantity x, Quantity y, Quantity z)
+        : BaseVector<dim>(cs, QuantityVector<dim>(x, y, z)) {}
+
     /*!
      * returns a QuantityVector with the components given in the "home"
      * CoordinateSystem of the Vector
      */
-    auto GetComponents() const
-    {
-        return BaseVector<dim>::qVector;
-    }
-    
+    auto GetComponents() const { return BaseVector<dim>::qVector; }
+
     /*!
      * returns a QuantityVector with the components given in an arbitrary
      * CoordinateSystem
      */
-    auto GetComponents(CoordinateSystem const& pCS) const
-    {
-        if (&pCS == BaseVector<dim>::cs)
-        {
-            return BaseVector<dim>::qVector;
-        }
-        else
-        {
-            return QuantityVector<dim>(CoordinateSystem::GetTransformation(*BaseVector<dim>::cs, pCS).linear() * BaseVector<dim>::qVector.eVector);
-        }
+    auto GetComponents(CoordinateSystem const& pCS) const {
+      if (&pCS == BaseVector<dim>::cs) {
+        return BaseVector<dim>::qVector;
+      } else {
+        return QuantityVector<dim>(
+            CoordinateSystem::GetTransformation(*BaseVector<dim>::cs, pCS).linear() *
+            BaseVector<dim>::qVector.eVector);
+      }
     }
-    
+
     /*!
      * transforms the Vector into another CoordinateSystem by changing
      * its components internally
      */
-    void rebase(CoordinateSystem const& pCS)
-    {
-        BaseVector<dim>::qVector = GetComponents(pCS);        
-        BaseVector<dim>::cs = &pCS;
+    void rebase(CoordinateSystem const& pCS) {
+      BaseVector<dim>::qVector = GetComponents(pCS);
+      BaseVector<dim>::cs = &pCS;
     }
-    
+
     /*!
      * returns the norm/length of the Vector. Before using this method,
      * think about whether squaredNorm() might be cheaper for your computation.
      */
-    auto norm() const
-    {
-        return BaseVector<dim>::qVector.norm();
-    }
-    
+    auto norm() const { return BaseVector<dim>::qVector.norm(); }
+
     /*!
      * returns the squared norm of the Vector. Before using this method,
      * think about whether norm() might be cheaper for your computation.
      */
-    auto squaredNorm() const
-    {
-        return BaseVector<dim>::qVector.squaredNorm();
-    }
-    
+    auto squaredNorm() const { return BaseVector<dim>::qVector.squaredNorm(); }
+
     /*!
      * returns a Vector \f$ \vec{v}_{\parallel} \f$ which is the parallel projection
      * of this vector \f$ \vec{v}_1 \f$ along another Vector \f$ \vec{v}_2 \f$ given by
@@ -94,114 +78,100 @@ public:
      *   \f]
      */
     template <typename dim2>
-    auto parallelProjectionOnto(Vector<dim2> const& pVec, CoordinateSystem const& pCS) const
-    {
-        auto const ourCompVec = GetComponents(pCS);
-        auto const otherCompVec = pVec.GetComponents(pCS);
-        auto const& a = ourCompVec.eVector;
-        auto const& b = otherCompVec.eVector;
-        
-        return Vector<dim>(pCS, QuantityVector<dim>(b * ((a.dot(b)) / b.squaredNorm())));
+    auto parallelProjectionOnto(Vector<dim2> const& pVec,
+                                CoordinateSystem const& pCS) const {
+      auto const ourCompVec = GetComponents(pCS);
+      auto const otherCompVec = pVec.GetComponents(pCS);
+      auto const& a = ourCompVec.eVector;
+      auto const& b = otherCompVec.eVector;
+
+      return Vector<dim>(pCS, QuantityVector<dim>(b * ((a.dot(b)) / b.squaredNorm())));
     }
-    
+
     template <typename dim2>
-    auto parallelProjectionOnto(Vector<dim2> const& pVec) const
-    {
-        return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
+    auto parallelProjectionOnto(Vector<dim2> const& pVec) const {
+      return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
     }
-    
-    auto operator+(Vector<dim> const& pVec) const
-    {
-        auto const components = GetComponents(*BaseVector<dim>::cs) + pVec.GetComponents(*BaseVector<dim>::cs);
-        return Vector<dim>(*BaseVector<dim>::cs, components);
+
+    auto operator+(Vector<dim> const& pVec) const {
+      auto const components =
+          GetComponents(*BaseVector<dim>::cs) + pVec.GetComponents(*BaseVector<dim>::cs);
+      return Vector<dim>(*BaseVector<dim>::cs, components);
     }
-    
-    auto operator-(Vector<dim> const& pVec) const
-    {
-        auto const components = GetComponents() - pVec.GetComponents(*BaseVector<dim>::cs);
-        return Vector<dim>(*BaseVector<dim>::cs, components);
+
+    auto operator-(Vector<dim> const& pVec) const {
+      auto const components = GetComponents() - pVec.GetComponents(*BaseVector<dim>::cs);
+      return Vector<dim>(*BaseVector<dim>::cs, components);
     }
-    
-    auto& operator*=(double const p)
-    {
-        BaseVector<dim>::qVector *= p;
-        return *this;
+
+    auto& operator*=(double const p) {
+      BaseVector<dim>::qVector *= p;
+      return *this;
     }
-    
+
     template <typename ScalarDim>
-    auto operator*(phys::units::quantity<ScalarDim, double> const p) const
-    {
-        using ProdQuantity = phys::units::detail::Product<dim, ScalarDim, double, double>;
-        
-        if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
-        {
-            return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
-        }
-        else
-        {
-            return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);        
-        }
+    auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
+      using ProdQuantity = phys::units::detail::Product<dim, ScalarDim, double, double>;
+
+      if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless,
+                                                               // not a "Quantity" anymore
+      {
+        return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs,
+                                                    BaseVector<dim>::qVector * p);
+      } else {
+        return Vector<typename ProdQuantity::dimension_type>(
+            *BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
+      }
     }
-    
+
     template <typename ScalarDim>
-    auto operator/(phys::units::quantity<ScalarDim, double> const p) const
-    {
-        return (*this) * (1 / p);
+    auto operator/(phys::units::quantity<ScalarDim, double> const p) const {
+      return (*this) * (1 / p);
     }
-    
-    auto operator*(double const p) const
-    {        
-        return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
-    }
-    
-    auto operator/(double const p) const
-    {        
-        return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector / p);
+
+    auto operator*(double const p) const {
+      return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
     }
-    
-    auto& operator+=(Vector<dim> const& pVec)
-    {
-        BaseVector<dim>::qVector += pVec.GetComponents(*BaseVector<dim>::cs);
-        return *this;
+
+    auto operator/(double const p) const {
+      return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector / p);
     }
-    
-    auto& operator-=(Vector<dim> const& pVec)
-    {
-        BaseVector<dim>::qVector -= pVec.GetComponents(*BaseVector<dim>::cs);
-        return *this;
+
+    auto& operator+=(Vector<dim> const& pVec) {
+      BaseVector<dim>::qVector += pVec.GetComponents(*BaseVector<dim>::cs);
+      return *this;
     }
-    
-    auto& operator-() const
-    {
-        return Vector<dim>(*BaseVector<dim>::cs, - BaseVector<dim>::qVector);
+
+    auto& operator-=(Vector<dim> const& pVec) {
+      BaseVector<dim>::qVector -= pVec.GetComponents(*BaseVector<dim>::cs);
+      return *this;
     }
-    
-    auto normalized() const
-    {
-        return (*this) * (1 / norm());
+
+    auto& operator-() const {
+      return Vector<dim>(*BaseVector<dim>::cs, -BaseVector<dim>::qVector);
     }
-    
+
+    auto normalized() const { return (*this) * (1 / norm()); }
+
     template <typename dim2>
-    auto cross(Vector<dim2> pV) const
-    {
-        auto const c1 = GetComponents().eVector;
-        auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
-        auto const bareResult = c1.cross(c2);
-        
-        using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
-        
-        if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
-        {
-            return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, bareResult);
-        }
-        else
-        {
-            return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs, bareResult);        
-        }
-    }
-    
-};
+    auto cross(Vector<dim2> pV) const {
+      auto const c1 = GetComponents().eVector;
+      auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
+      auto const bareResult = c1.cross(c2);
+
+      using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
+
+      if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless,
+                                                               // not a "Quantity" anymore
+      {
+        return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, bareResult);
+      } else {
+        return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs,
+                                                             bareResult);
+      }
+    }
+  };
+
+} // namespace fwk
 
-} // end namespace
-  
 #endif
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index ec0ab2c8e03c3cd3076301d0c4e97a74fd0557df..2af6e2c50d21a5008f13c65d4ec490ea5eb363d8 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -1,7 +1,7 @@
 #define CATCH_CONFIG_MAIN  // This tells Catch to provide a main() - only do this in one cpp file
 #include <catch2/catch.hpp>
 
-#include <Units/PhysicalUnits.h>
+#include <fwk/PhysicalUnits.h>
 
 using namespace phys::units;
 using namespace phys::units::literals;
diff --git a/Framework/ParticleStack/CMakeLists.txt b/Framework/ParticleStack/CMakeLists.txt
index 66534852893a0aa85e45bc4438170a30444ec06f..6b9547ec2a79e66e5b124a7b2cf5aa72b8b70e68 100644
--- a/Framework/ParticleStack/CMakeLists.txt
+++ b/Framework/ParticleStack/CMakeLists.txt
@@ -1,14 +1,22 @@
 
-set (PUBLIC_HEADER_FILES StackOne.h)
+set (CORSIKAstack_HEADERS StackOne.h)
+set (CORSIKAstack_NAMESPACE fwk)
 
 add_library (CORSIKAstack INTERFACE)
 
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAstack ${CORSIKAstack_NAMESPACE} ${CORSIKAstack_HEADERS})
+
 #set_target_properties (CORSIKAstack PROPERTIES VERSION ${PROJECT_VERSION})
 #set_target_properties (CORSIKAstack PROPERTIES SOVERSION 1)
 
 #set_target_properties (CORSIKAstack PROPERTIES PUBLIC_HEADER "${STACK_HEADERS}")
 
-#target_link_libraries (CORSIKAstackinterface CORSIKAunits)
+target_link_libraries (
+  CORSIKAstack
+  INTERFACE
+  CORSIKAstackint
+  CORSIKAunits
+  )
 
 #target_include_directories (CORSIKAstack PRIVATE   ${EIGEN3_INCLUDE_DIR})
 target_include_directories (
@@ -21,12 +29,12 @@ target_include_directories (
 #add_custom_command (
 #  CORSIKAstack
 #  PRE_BUILD
-#  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${PUBLIC_HEADER_FILES} ${PROJECT_BINARY_DIR}/include/corsika})
+#  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CORSIKAstack_HEADERS} ${PROJECT_BINARY_DIR}/include/corsika})
 #  )
 
 install (
   FILES
   StackOne.h
   DESTINATION
-  include/Stack
+  include/${CORSIKAstack_NAMESPACE}
   )
diff --git a/Framework/ParticleStack/StackOne.h b/Framework/ParticleStack/StackOne.h
index e41007a6a97aea7ba0d91402df5cd022372778b0..e7d4285f738274515c8a538e61994dc3beb89349 100644
--- a/Framework/ParticleStack/StackOne.h
+++ b/Framework/ParticleStack/StackOne.h
@@ -1,62 +1,53 @@
 #ifndef _include_stackone_h_
 #define _include_stackone_h_
 
-#include <vector>
 #include <string>
+#include <vector>
 
-#include <StackInterface/Stack.h>
-//#include <corsika/Stack.h>
-
+#include <fwk/Stack.h>
 
 namespace stack {
-  
-  
+
   /**
      Example of a particle object on the stack.
    */
-  
-  template<typename _Stack>
-  class ParticleReadOne : public StackIteratorInfo<_Stack, ParticleReadOne<_Stack> >
-    {
-      using StackIteratorInfo<_Stack, ParticleReadOne>::GetIndex;
-      using StackIteratorInfo<_Stack, ParticleReadOne>::GetStack;
-      
-    public:
-      void SetId(const int id) { GetStack().SetId(GetIndex(), id); }
-      void SetEnergy(const double e) { GetStack().SetEnergy(GetIndex(), e); }
-      
-      int GetId() const { return  GetStack().GetId(GetIndex()); }
-      double GetEnergy() const { return GetStack().GetEnergy(GetIndex()); }
-      
-      double GetPDG() const { return 0; } // ConvertToPDG(GetId()); }  
-      void SetPDG(double v) { GetStack().SetId(0, 0); } //fIndex, ConvertFromPDG(v)); }
-    };
-  
-
- 
-  
+
+  template <typename _Stack>
+  class ParticleReadOne : public StackIteratorInfo<_Stack, ParticleReadOne<_Stack> > {
+    using StackIteratorInfo<_Stack, ParticleReadOne>::GetIndex;
+    using StackIteratorInfo<_Stack, ParticleReadOne>::GetStack;
+
+  public:
+    void SetId(const int id) { GetStack().SetId(GetIndex(), id); }
+    void SetEnergy(const double e) { GetStack().SetEnergy(GetIndex(), e); }
+
+    int GetId() const { return GetStack().GetId(GetIndex()); }
+    double GetEnergy() const { return GetStack().GetEnergy(GetIndex()); }
+
+    double GetPDG() const { return 0; }               // ConvertToPDG(GetId()); }
+    void SetPDG(double v) { GetStack().SetId(0, 0); } // fIndex, ConvertFromPDG(v)); }
+  };
 
   /**
      Memory implementation of the most simple particle stack object.
    */
-  
-  class StackOneImpl
-  {    
+
+  class StackOneImpl {
   private:
     /// the actual memory to store particle data
 
     std::vector<int> fId;
     std::vector<double> fData;
-    
+
   public:
     void Clear() { fData.clear(); }
-    
-    int GetSize() const { return fData.size();  }
+
+    int GetSize() const { return fData.size(); }
     int GetCapacity() const { return fData.size(); }
-        
+
     void SetId(const int i, const int id) { fId[i] = id; }
     void SetEnergy(const int i, const double e) { fData[i] = e; }
-    
+
     const int GetId(const int i) const { return fId[i]; }
     const double GetEnergy(const int i) const { return fData[i]; }
 
@@ -67,15 +58,23 @@ namespace stack {
       fData[i2] = fData[i1];
       fId[i2] = fId[i1];
     }
-    
+
   protected:
-    void IncrementSize() { fData.push_back(0.); fId.push_back(0.); }
-    void DecrementSize() { if (fData.size()>0) { fData.pop_back(); fId.pop_back(); } }
+    void IncrementSize() {
+      fData.push_back(0.);
+      fId.push_back(0.);
+    }
+    void DecrementSize() {
+      if (fData.size() > 0) {
+        fData.pop_back();
+        fId.pop_back();
+      }
+    }
   };
-  
+
   typedef StackIterator<StackOneImpl, ParticleReadOne<StackOneImpl> > ParticleOne;
   typedef Stack<StackOneImpl, ParticleOne> StackOne;
-  
-} // end namespace
-  
+
+} // namespace stack
+
 #endif
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index 2d65ec9277b812b13560010207399303a7fd1014..c31aa2612e49d189e13f50218b8198f915df0fcd 100644
--- a/Framework/Particles/testParticles.cc
+++ b/Framework/Particles/testParticles.cc
@@ -1,7 +1,8 @@
-#define CATCH_CONFIG_MAIN  // This tells Catch to provide a main() - only do this in one cpp file
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+                          // cpp file
 #include <catch2/catch.hpp>
 
-#include <Units/PhysicalUnits.h>
+#include <fwk/PhysicalUnits.h>
 
 #include <fwk/Particles.h>
 
@@ -10,21 +11,16 @@ using namespace phys::units::literals;
 
 using namespace fwk::particle;
 
-TEST_CASE( "Particles", "[Particles]" )
-{  
-  SECTION( "Types" )
-    {
-      REQUIRE( Electron::GetType()==InternalParticleCode::Electron );
-    }
-
-    SECTION( "Data" )
-    {
-      REQUIRE( Electron::GetMass()/0.511_MeV==Approx(1) );
-      REQUIRE( Electron::GetMass()/GetMass(InternalParticleCode::Electron)==Approx(1) );
-      REQUIRE( Electron::GetCharge()/phys::units::e==Approx(-1) );
-      REQUIRE( Positron::GetCharge()/phys::units::e==Approx(+1) );
-      REQUIRE( GetElectricCharge(Positron::GetAntiParticle())/phys::units::e==Approx(-1) );
-      REQUIRE( Electron::GetName() == "e-" );
-    }
+TEST_CASE("Particles", "[Particles]") {
+  SECTION("Types") { REQUIRE(Electron::GetType() == InternalParticleCode::Electron); }
 
+  SECTION("Data") {
+    REQUIRE(Electron::GetMass() / 0.511_MeV == Approx(1));
+    REQUIRE(Electron::GetMass() / GetMass(InternalParticleCode::Electron) == Approx(1));
+    REQUIRE(Electron::GetCharge() / phys::units::e == Approx(-1));
+    REQUIRE(Positron::GetCharge() / phys::units::e == Approx(+1));
+    REQUIRE(GetElectricCharge(Positron::GetAntiParticle()) / phys::units::e ==
+            Approx(-1));
+    REQUIRE(Electron::GetName() == "e-");
+  }
 }
diff --git a/Framework/StackInterface/CMakeLists.txt b/Framework/StackInterface/CMakeLists.txt
index 21e156044e29f420d7e0442c83c60b60b37aeb14..ba9ef212faec93809bcbe5d0a80c3754a9009168 100644
--- a/Framework/StackInterface/CMakeLists.txt
+++ b/Framework/StackInterface/CMakeLists.txt
@@ -1,18 +1,25 @@
 set (
-  PUBLIC_HEADER_FILES
+  CORSIKAstackint_HEADERS
   Stack.h
   StackIterator.h
   )
 
+set (
+  CORSIKAstackint_NAMESPACE
+  fwk
+  )
+
 add_library (
-  CORSIKAstackinterface
+  CORSIKAstackint
   INTERFACE
   )
 
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAstackint ${CORSIKAstackint_NAMESPACE} ${CORSIKAstackint_HEADERS})
+
 target_include_directories (
-  CORSIKAstackinterface
+  CORSIKAstackint
   INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include>
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
   $<INSTALL_INTERFACE:include>
   )
 
@@ -20,12 +27,12 @@ target_include_directories (
 #  TARGET
 #  CORSIKAstackinterface
 #  PRE_BUILD
-#  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${PUBLIC_HEADER_FILES} ${PROJECT_BINARY_DIR}/include/corsika
+#  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CORSIKAstackinterface_HEADERS} ${PROJECT_BINARY_DIR}/include/corsika
 #  )
 
 install (
   FILES
-  ${PUBLIC_HEADER_FILES}
+  ${CORSIKAstackint_HEADERS}
   DESTINATION
-  include/corsika
+  include/${CORSIKAstackint_NAMESPACE}
   )
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index 9412286bf85fc9558308e01d117124c181bbba7b..14627ba2f4a7e4ca7d929dd1be9c471d220c0be7 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -1,8 +1,7 @@
 #ifndef _include_Stack_h__
 #define _include_Stack_h__
 
-#include <StackInterface/StackIterator.h> // to help application programmres
-//#include <corsika/StackIterator.h> // to help application programmres
+#include <fwk/StackIterator.h> // to help application programmres
 
 /**
    All classes around management of particles on a stack.
@@ -15,41 +14,43 @@ namespace stack {
      std-type begin/end function to allow integration in normal for
      loops etc.
    */
-  
-  template<typename DataImpl, typename Particle> 
+
+  template <typename DataImpl, typename Particle>
   class Stack : public DataImpl {
 
   public:
     using DataImpl::GetCapacity;
     using DataImpl::GetSize;
-    
+
     using DataImpl::Clear;
     using DataImpl::Copy;
-    
-    using DataImpl::IncrementSize;
+
     using DataImpl::DecrementSize;
-    
-  public:  
+    using DataImpl::IncrementSize;
+
+  public:
     typedef Particle iterator;
     typedef const Particle const_iterator;
 
     /// these are functions required by std containers and std loops
-    iterator begin() { return iterator(*this, 0); } 
-    iterator end() { return iterator(*this, GetSize()); } 
-    iterator last() { return iterator(*this, GetSize()-1); } 
-    
+    iterator begin() { return iterator(*this, 0); }
+    iterator end() { return iterator(*this, GetSize()); }
+    iterator last() { return iterator(*this, GetSize() - 1); }
+
     /// these are functions required by std containers and std loops
-    const_iterator cbegin() const { return const_iterator(*this, 0); } 
-    const_iterator cend() const { return const_iterator(*this, GetSize()); } 
-    const_iterator clast() const { return const_iterator(*this, GetSize()-1); } 
+    const_iterator cbegin() const { return const_iterator(*this, 0); }
+    const_iterator cend() const { return const_iterator(*this, GetSize()); }
+    const_iterator clast() const { return const_iterator(*this, GetSize() - 1); }
 
     /// increase stack size, create new particle at end of stack
-    iterator NewParticle() { IncrementSize(); return iterator(*this, GetSize()-1); }
+    iterator NewParticle() {
+      IncrementSize();
+      return iterator(*this, GetSize() - 1);
+    }
     /// delete last particle on stack by decrementing stack size
     void DeleteLast() { DecrementSize(); }
   };
 
-} // end namespace
-  
+} // namespace stack
+
 #endif
-  
diff --git a/Framework/StackInterface/StackIterator.h b/Framework/StackInterface/StackIterator.h
index b2e9a8b6c523b8330731e729e5a7321a130fe3c5..ae8b924c6396429c5bb8d9d3484bf7eae265a9fb 100644
--- a/Framework/StackInterface/StackIterator.h
+++ b/Framework/StackInterface/StackIterator.h
@@ -1,17 +1,18 @@
 #ifndef _include_StackIterator_h__
 #define _include_StackIterator_h__
 
-#include <iostream>
 #include <iomanip>
+#include <iostream>
 
 namespace stack {
 
   // forward decl.
-  template<class Stack, class Particle> class StackIteratorInfo;
+  template <class Stack, class Particle>
+  class StackIteratorInfo;
 
   /**
      @class StackIterator
-     
+
      The StackIterator is the main interface to iterator over
      particles on a stack. At the same time StackIterator is a
      Particle object by itself, thus there is no difference between
@@ -19,78 +20,97 @@ namespace stack {
 
      This allows to write code like
      \verbatim
-     for (auto& p : theStack) { p.SetEnergy(newEnergy); }  
+     for (auto& p : theStack) { p.SetEnergy(newEnergy); }
      \endverbatim
 
      The template argument Stack determines the type of Stack object
      the data is stored in. A pointer to the Stack object is part of
      the StackIterator. In addition to Stack the iterator only knows
-     the index fIndex in the Stack data. 
+     the index fIndex in the Stack data.
 
      The template argument Particles acts as a policy to provide
      readout function of Particle data from the stack. The Particle
      class must know how to retrieve information from the Stack data
      for a particle entry at any index fIndex.
   */
-  
-  template<typename Stack, typename Particle>
-  class StackIterator : public Particle
-  {
+
+  template <typename Stack, typename Particle>
+  class StackIterator : public Particle {
     friend Stack;
     friend Particle;
-    friend StackIteratorInfo<Stack,Particle>;
-    
+    friend StackIteratorInfo<Stack, Particle>;
+
   private:
     int fIndex;
-    
+
     //#warning stacks should not be copied because of this:
     Stack* fData;
-    
+
   public:
-    StackIterator() : fData(0), fIndex(0) { }
-    StackIterator(Stack& data, const int index) : fData(&data), fIndex(index) { }
-    StackIterator(const StackIterator& mit) : fData(mit.fData), fIndex(mit.fIndex) { }
-    
-    StackIterator& operator++() { ++fIndex; return *this; }
-    StackIterator operator++(int) { StackIterator tmp(*this); ++fIndex; return tmp; }
+    StackIterator()
+        : fData(0)
+        , fIndex(0) {}
+    StackIterator(Stack& data, const int index)
+        : fData(&data)
+        , fIndex(index) {}
+    StackIterator(const StackIterator& mit)
+        : fData(mit.fData)
+        , fIndex(mit.fIndex) {}
+
+    StackIterator& operator++() {
+      ++fIndex;
+      return *this;
+    }
+    StackIterator operator++(int) {
+      StackIterator tmp(*this);
+      ++fIndex;
+      return tmp;
+    }
     bool operator==(const StackIterator& rhs) { return fIndex == rhs.fIndex; }
     bool operator!=(const StackIterator& rhs) { return fIndex != rhs.fIndex; }
-    
+
     StackIterator& operator*() { return *this; }
     const StackIterator& operator*() const { return *this; }
-    
+
   protected:
     int GetIndex() const { return fIndex; }
     Stack& GetStack() { return *fData; }
     const Stack& GetStack() const { return *fData; }
 
     // this is probably not needed rigth now:
-    //inline StackIterator<Stack,Particle>& BaseRef() { return static_cast<StackIterator<Stack, Particle>&>(*this); }
-    //inline const StackIterator<Stack,Particle>& BaseRef() const { return static_cast<const StackIterator<Stack, Particle>&>(*this); }
+    // inline StackIterator<Stack,Particle>& BaseRef() { return
+    // static_cast<StackIterator<Stack, Particle>&>(*this); } inline const
+    // StackIterator<Stack,Particle>& BaseRef() const { return static_cast<const
+    // StackIterator<Stack, Particle>&>(*this); }
   };
 
-  
-
   /**
      @class StackIteratorInfo
-     
-     This is the class where custom ... 
+
+     This is the class where custom ...
      Internal helper class for StackIterator. Document better...
    */
-  
-  template<typename _Stack, typename Particle>
+
+  template <typename _Stack, typename Particle>
   class StackIteratorInfo {
-    
-    friend Particle;  
+
+    friend Particle;
+
   private:
     StackIteratorInfo() {}
-    
+
   protected:
-    inline _Stack& GetStack() { return static_cast<StackIterator<_Stack, Particle>*>(this)->GetStack(); }
-    inline int GetIndex() const { return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetIndex(); }
-    inline const _Stack& GetStack() const { return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetStack(); }
+    inline _Stack& GetStack() {
+      return static_cast<StackIterator<_Stack, Particle>*>(this)->GetStack();
+    }
+    inline int GetIndex() const {
+      return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetIndex();
+    }
+    inline const _Stack& GetStack() const {
+      return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetStack();
+    }
   };
 
 } // end namespace stack
-  
+
 #endif
diff --git a/Framework/Units/CMakeLists.txt b/Framework/Units/CMakeLists.txt
index 7881afc4f9264391915ca73b5a71dd21b06f9054..d0818953175fa1edc576fccf1475a01c795c0ba9 100644
--- a/Framework/Units/CMakeLists.txt
+++ b/Framework/Units/CMakeLists.txt
@@ -1,14 +1,23 @@
 
 add_library (CORSIKAunits INTERFACE)
 
+set (CORSIKAunits_NAMESPACE fwk)
+set (
+  CORSIKAunits_HEADERS
+  PhysicalUnits.h
+  PhysicalConstants.h
+  )
+
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAunits ${CORSIKAunits_NAMESPACE} ${CORSIKAunits_HEADERS})
+
 target_include_directories (CORSIKAunits
   INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework>
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
   $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/ThirdParty>
   $<INSTALL_INTERFACE:include>
   )
 
-install (FILES PhysicalUnits.h DESTINATION include/Units)
+install (FILES ${CORSIKAunits_HEADERS} DESTINATION include/${CORSIKAunits_NAMESPACE})
 
 # code testing
 add_executable (testUnits testUnits.cc)
diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc
index 5cf6e93c87bfb9304dd48e9ac4282699553b1ad2..e963c6796ee5d746abc0c29bd006c5d07c394d0e 100644
--- a/Framework/Units/testUnits.cc
+++ b/Framework/Units/testUnits.cc
@@ -1,7 +1,7 @@
 #define CATCH_CONFIG_MAIN  // This tells Catch to provide a main() - only do this in one cpp file
 #include <catch2/catch.hpp>
 
-#include <Units/PhysicalUnits.h>
+#include <fwk/PhysicalUnits.h>
 
 #include <array>