Environment/IMediumModel.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_IMediumModel_h
-#define _include_IMediumModel_h
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::environment {
-
-  class IMediumModel {
-  public:
-    virtual ~IMediumModel() = default;
-
-    virtual corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const&) const = 0;
-
-    // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
-    // approach for now, only lines are supported
-    virtual corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::LengthType) const = 0;
-
-    virtual corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::GrammageType) const = 0;
-
-    virtual NuclearComposition const& GetNuclearComposition() const = 0;
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/InhomogeneousMedium.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_environment_InhomogeneousMedium_h_
-#define _include_environment_InhomogeneousMedium_h_
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-
-/**
- * A general inhomogeneous medium. The mass density distribution TDensityFunction must be
- * a \f$C^2\f$-function.
- */
-
-namespace corsika::environment {
-
-  template <class T, class TDensityFunction>
-  class InhomogeneousMedium : public T {
-    NuclearComposition const fNuclComp;
-    TDensityFunction const fDensityFunction;
-
-  public:
-    template <typename... Args>
-    InhomogeneousMedium(NuclearComposition pNuclComp, Args&&... rhoArgs)
-        : fNuclComp(pNuclComp)
-        , fDensityFunction(rhoArgs...) {}
-
-    corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const& p) const override {
-      return fDensityFunction.EvaluateAt(p);
-    }
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
-        corsika::units::si::LengthType pTo) const override {
-      return fDensityFunction.IntegrateGrammage(pLine, pTo);
-    }
-
-    corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
-        corsika::units::si::GrammageType pGrammage) const override {
-      return fDensityFunction.ArclengthFromGrammage(pLine, pGrammage);
-    }
-  };
-
-} // namespace corsika::environment
-#endif

Environment/LinearApproximationIntegrator.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_environment_LinearApproximationIntegrator_h_
-#define _include_environment_LinearApproximationIntegrator_h_
-
-#include <limits>
-
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-
-namespace corsika::environment {
-  template <class TDerived>
-  class LinearApproximationIntegrator {
-    auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
-
-  public:
-    auto IntegrateGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-        corsika::units::si::LengthType length) const {
-      auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0));
-      auto const c1 = GetImplementation().fRho.FirstDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-      return (c0 + 0.5 * c1 * length) * length;
-    }
-
-    auto ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-        corsika::units::si::GrammageType grammage) const {
-      auto const c0 = GetImplementation().fRho(line.GetPosition(0));
-      auto const c1 = GetImplementation().fRho.FirstDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-
-      return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
-    }
-
-    auto MaximumLength(corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-                       [[maybe_unused]] double relError) const {
-      using namespace corsika::units::si;
-      [[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-
-      // todo: provide a real, working implementation
-      return 1_m * std::numeric_limits<double>::infinity();
-    }
-  };
-} // namespace corsika::environment
-
-#endif

Environment/NameModel.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_NameModel_h
-#define _include_NameModel_h
-
-#include <string>
-#include <utility>
-
-namespace corsika::environment {
-
-  template <typename T>
-  struct NameModel : public T {
-    virtual std::string const& GetName() const = 0;
-    virtual ~NameModel() = default;
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/NuclearComposition.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_NuclearComposition_h
-#define _include_NuclearComposition_h
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <cassert>
-#include <numeric>
-#include <random>
-#include <stdexcept>
-#include <vector>
-
-namespace corsika::environment {
-  class NuclearComposition {
-    std::vector<float> const fNumberFractions; //!< relative fractions of number density
-    std::vector<corsika::particles::Code> const
-        fComponents; //!< particle codes of consitutents
-
-    double const fAvgMassNumber;
-
-    template <class AConstIterator, class BConstIterator>
-    class WeightProviderIterator {
-      AConstIterator fAIter;
-      BConstIterator fBIter;
-
-    public:
-      using value_type = double;
-      using iterator_category = std::input_iterator_tag;
-      using pointer = value_type*;
-      using reference = value_type&;
-      using difference_type = ptrdiff_t;
-
-      WeightProviderIterator(AConstIterator a, BConstIterator b)
-          : fAIter(a)
-          , fBIter(b) {}
-
-      value_type operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
-
-      WeightProviderIterator& operator++() { // prefix ++
-        ++fAIter;
-        ++fBIter;
-        return *this;
-      }
-
-      auto operator==(WeightProviderIterator other) { return fAIter == other.fAIter; }
-
-      auto operator!=(WeightProviderIterator other) { return !(*this == other); }
-    };
-
-  public:
-    NuclearComposition(std::vector<corsika::particles::Code> pComponents,
-                       std::vector<float> pFractions)
-        : fNumberFractions(pFractions)
-        , fComponents(pComponents)
-        , fAvgMassNumber(std::inner_product(
-              pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
-              [](double x, double y) { return x + y; },
-              [](auto const& compID, auto const& fraction) {
-                return corsika::particles::GetNucleusA(compID) * fraction;
-              })) {
-      assert(pComponents.size() == pFractions.size());
-      auto const sumFractions =
-          std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
-
-      if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
-        throw std::runtime_error("element fractions do not add up to 1");
-      }
-    }
-
-    template <typename TFunction>
-    auto WeightedSum(TFunction func) const {
-      using ResultQuantity = decltype(func(*fComponents.cbegin()));
-
-      auto const sum = [](auto x, auto y) { return x + y; };
-      auto const prod = [&](auto const compID, auto const fraction) {
-        return func(compID) * fraction;
-      };
-
-      if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
-        return std::inner_product(
-            fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
-            ResultQuantity::zero(), // .zero() is defined for quantity types only
-            sum, prod);
-      } else {
-        return std::inner_product(
-            fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
-            ResultQuantity(0), // in other cases we have to use a bare 0
-            sum, prod);
-      }
-    }
-
-    auto size() const { return fNumberFractions.size(); }
-
-    auto const& GetFractions() const { return fNumberFractions; }
-    auto const& GetComponents() const { return fComponents; }
-    auto const GetAverageMassNumber() const { return fAvgMassNumber; }
-
-    template <class TRNG>
-    corsika::particles::Code SampleTarget(
-        std::vector<corsika::units::si::CrossSectionType> const& sigma,
-        TRNG& randomStream) const {
-      using namespace corsika::units::si;
-
-      assert(sigma.size() == fNumberFractions.size());
-
-      std::discrete_distribution channelDist(
-          WeightProviderIterator<decltype(fNumberFractions.begin()),
-                                 decltype(sigma.begin())>(fNumberFractions.begin(),
-                                                          sigma.begin()),
-          WeightProviderIterator<decltype(fNumberFractions.begin()),
-                                 decltype(sigma.end())>(fNumberFractions.end(),
-                                                        sigma.end()));
-
-      auto const iChannel = channelDist(randomStream);
-      return fComponents[iChannel];
-    }
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/SlidingPlanarExponential.h

-/*
- * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_Environment_SlidingPlanarExponential_h_
-#define _include_Environment_SlidingPlanarExponential_h_
-
-#include <corsika/environment/FlatExponential.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::environment {
-
-  // clang-format off
-  /**
-   * The SlidingPlanarExponential models mass density as
-   * \f[
-   *   \varrho(r) = \varrho_0 \exp\left( \frac{|p_0 - r|}{\lambda} \right).
-   * \f]
-   * For grammage/length conversion, the density distribution is approximated as
-   * locally flat at the starting point \f$ r_0 \f$ of the trajectory with the axis pointing
-   * from \f$ p_0 \f$ to \f$ r_0 \f$.
-   */
-   //clang-format on
-   
-  template <class T>
-  class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
-                                   public T {
-    NuclearComposition const fNuclComp;
-
-    using Base = BaseExponential<SlidingPlanarExponential<T>>;
-
-  public:
-    SlidingPlanarExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
-                             units::si::LengthType vLambda, NuclearComposition vNuclComp)
-        : Base(vP0, vRho, vLambda)
-        , fNuclComp(vNuclComp) {}
-
-    units::si::MassDensityType GetMassDensity(
-        geometry::Point const& vP) const override {
-      auto const height = (vP - Base::fP0).norm();
-      return Base::fRho0 * exp(Base::fInvLambda * height);
-    }
-
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    units::si::GrammageType IntegratedGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::LengthType vL) const override {
-      auto const axis = (vLine.GetR0() - Base::fP0).normalized();
-      return Base::IntegratedGrammage(vLine, vL, axis);
-    }
-
-    units::si::LengthType ArclengthFromGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::GrammageType vGrammage) const override {
-      auto const axis = (vLine.GetR0() - Base::fP0).normalized();
-      return Base::ArclengthFromGrammage(vLine, vGrammage, axis);
-    }
-  };
-
-} // namespace corsika::environment
-#endif

Environment/VolumeTreeNode.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_VolumeTreeNode_H
-#define _include_VolumeTreeNode_H
-
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/geometry/Volume.h>
-#include <memory>
-#include <vector>
-
-namespace corsika::environment {
-
-  class Empty {}; //<! intended for usage as default template argument
-
-  template <typename TModelProperties = Empty>
-  class VolumeTreeNode {
-  public:
-    using IModelProperties = TModelProperties;
-    using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
-    using IMPSharedPtr = std::shared_ptr<IModelProperties>;
-    using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
-
-    VolumeTreeNode(VolUPtr pVolume = nullptr)
-        : fGeoVolume(std::move(pVolume)) {}
-
-    //! convenience function equivalent to Volume::Contains
-    bool Contains(corsika::geometry::Point const& p) const {
-      return fGeoVolume->Contains(p);
-    }
-
-    VolumeTreeNode<IModelProperties> const* Excludes(
-        corsika::geometry::Point const& p) const {
-      auto exclContainsIter =
-          std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
-                       [&](auto const& s) { return bool(s->Contains(p)); });
-
-      return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
-    }
-
-    /** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
-     * \class Point \p p, or nullptr iff \p p is not contained in this volume.
-     */
-    VolumeTreeNode<IModelProperties> const* GetContainingNode(
-        corsika::geometry::Point const& p) const {
-      if (!Contains(p)) { return nullptr; }
-
-      if (auto const childContainsIter =
-              std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
-                           [&](auto const& s) { return bool(s->Contains(p)); });
-          childContainsIter == fChildNodes.cend()) // not contained in any of the children
-      {
-        if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
-        {
-          return exclContainsIter->GetContainingNode(p);
-        } else {
-          return this;
-        }
-      } else {
-        return (*childContainsIter)->GetContainingNode(p);
-      }
-    }
-
-    /**
-     * Traverses the VolumeTree pre- or post-order and calls the functor  \p func for each
-     * node. \p func takes a reference to VolumeTreeNode as argument. The return value \p
-     * func is ignored.
-     */
-    template <typename TCallable, bool preorder = true>
-    void walk(TCallable func) {
-      if constexpr (preorder) { func(*this); }
-
-      std::for_each(fChildNodes.begin(), fChildNodes.end(),
-                    [&](auto& v) { v->walk(func); });
-
-      if constexpr (!preorder) { func(*this); };
-    }
-
-    void AddChild(VTNUPtr pChild) {
-      pChild->fParentNode = this;
-      fChildNodes.push_back(std::move(pChild));
-      // It is a bad idea to return an iterator to the inserted element
-      // because it might get invalidated when the vector needs to grow
-      // later and the caller won't notice.
-    }
-
-    void ExcludeOverlapWith(VTNUPtr const& pNode) {
-      fExcludedNodes.push_back(pNode.get());
-    }
-
-    auto* GetParent() const { return fParentNode; };
-
-    auto const& GetChildNodes() const { return fChildNodes; }
-
-    auto const& GetExcludedNodes() const { return fExcludedNodes; }
-
-    auto const& GetVolume() const { return *fGeoVolume; }
-
-    auto const& GetModelProperties() const { return *fModelProperties; }
-
-    bool HasModelProperties() const { return fModelProperties.get() != nullptr; }
-
-    template <typename ModelProperties, typename... Args>
-    auto SetModelProperties(Args&&... args) {
-      static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
-                    "unusable type provided");
-
-      fModelProperties = std::make_shared<ModelProperties>(std::forward<Args>(args)...);
-      return fModelProperties;
-    }
-
-    void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
-
-    template <class MediumType, typename... Args>
-    static auto CreateMedium(Args&&... args) {
-      static_assert(std::is_base_of_v<IMediumModel, MediumType>,
-                    "unusable type provided, needs to be derived from \"IMediumModel\"");
-
-      return std::make_shared<MediumType>(std::forward<Args>(args)...);
-    }
-
-  private:
-    std::vector<VTNUPtr> fChildNodes;
-    std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
-    VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
-    VolUPtr fGeoVolume;
-    IMPSharedPtr fModelProperties;
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/testEnvironment.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-
-#include <corsika/environment/DensityFunction.h>
-#include <corsika/environment/FlatExponential.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/InhomogeneousMedium.h>
-#include <corsika/environment/LinearApproximationIntegrator.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/SlidingPlanarExponential.h>
-#include <corsika/environment/VolumeTreeNode.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <catch2/catch.hpp>
-
-using namespace corsika::geometry;
-using namespace corsika::environment;
-using namespace corsika::particles;
-using namespace corsika::units::si;
-using namespace corsika;
-
-CoordinateSystem const& gCS =
-    RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-Point const gOrigin(gCS, {0_m, 0_m, 0_m});
-
-TEST_CASE("HomogeneousMedium") {
-  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
-                                             std::vector<float>{1.f});
-  HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
-}
-
-TEST_CASE("FlatExponential") {
-  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
-                                             std::vector<float>{1.f});
-
-  Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
-  LengthType const lambda = 3_m;
-  auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
-  FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
-                                             protonComposition);
-  auto const tEnd = 5_s;
-
-  SECTION("horizontal") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {20_cm / second, 0_m / second, 0_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-
-    REQUIRE((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
-    REQUIRE((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
-  }
-
-  SECTION("vertical") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 0_m / second, 5_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-    LengthType const length = 2 * lambda;
-    GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
-
-    REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
-    REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
-  }
-
-  SECTION("escape grammage") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 0_m / second, -5_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-
-    GrammageType const escapeGrammage = rho0 * lambda;
-
-    REQUIRE(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
-    REQUIRE(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
-            std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
-  }
-
-  SECTION("inclined") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 5_m / second, 5_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-    double const cosTheta = M_SQRT1_2;
-    LengthType const length = 2 * lambda;
-    GrammageType const exact =
-        rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
-    REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
-    REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
-  }
-}
-
-TEST_CASE("SlidingPlanarExponential") {
-  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
-                                             std::vector<float>{1.f});
-
-  LengthType const lambda = 3_m;
-  auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
-  auto const tEnd = 5_s;
-
-  SlidingPlanarExponential<IMediumModel> const medium(gOrigin, rho0, lambda,
-                                                      protonComposition);
-
-  SECTION("density") {
-    CHECK(medium.GetMassDensity({gCS, {0_m, 0_m, 3_m}}) /
-              medium.GetMassDensity({gCS, {0_m, 3_m, 0_m}}) ==
-          Approx(1));
-  }
-
-  SECTION("vertical") {
-    Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
-    FlatExponential<IMediumModel> const flat(gOrigin, axis, rho0, lambda,
-                                             protonComposition);
-    Line const line({gCS, {0_m, 0_m, 1_m}},
-                    Vector<SpeedType::dimension_type>(
-                        gCS, {0_m / second, 0_m / second, 5_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-
-    CHECK(medium.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
-          flat.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
-    CHECK(medium.IntegratedGrammage(trajectory, 2_m).magnitude() ==
-          flat.IntegratedGrammage(trajectory, 2_m).magnitude());
-    CHECK(medium.ArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude() ==
-          flat.ArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude());
-  }
-}
-
-auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
-
-struct Exponential {
-  auto operator()(corsika::geometry::Point const& p) const {
-    return exp(p.GetCoordinates()[0] / 1_m) * rho0;
-  }
-
-  template <int N>
-  auto Derivative(Point const& p, Vector<dimensionless_d> const& v) const {
-    return v.GetComponents()[0] * (*this)(p) /
-           corsika::units::si::detail::static_pow<N>(1_m);
-  }
-
-  auto FirstDerivative(Point const& p, Vector<dimensionless_d> const& v) const {
-    return Derivative<1>(p, v);
-  }
-
-  auto SecondDerivative(Point const& p, Vector<dimensionless_d> const& v) const {
-    return Derivative<2>(p, v);
-  }
-};
-
-TEST_CASE("InhomogeneousMedium") {
-  Vector direction(gCS, QuantityVector<dimensionless_d>(1, 0, 0));
-
-  Line line(gOrigin, Vector<SpeedType::dimension_type>(
-                         gCS, {20_m / second, 0_m / second, 0_m / second}));
-
-  auto const tEnd = 5_s;
-  Trajectory<Line> const trajectory(line, tEnd);
-
-  Exponential const e;
-  DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
-
-  SECTION("DensityFunction") {
-    REQUIRE(e.Derivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
-            Approx(1));
-    REQUIRE(rho.EvaluateAt(gOrigin) == e(gOrigin));
-  }
-
-  auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };
-  auto const exactLength = [](auto X) { return 1_m * log(1 + X / (rho0 * 1_m)); };
-
-  auto constexpr l = 15_cm;
-
-  NuclearComposition const composition{{Code::Proton}, {1.f}};
-  InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
-
-  SECTION("Integration") {
-    REQUIRE(rho.IntegrateGrammage(trajectory, l) / exactGrammage(l) ==
-            Approx(1).epsilon(1e-2));
-    REQUIRE(rho.ArclengthFromGrammage(trajectory, exactGrammage(l)) /
-                exactLength(exactGrammage(l)) ==
-            Approx(1).epsilon(1e-2));
-    REQUIRE(rho.MaximumLength(trajectory, 1e-2) >
-            l); // todo: write reasonable test when implementation is working
-
-    REQUIRE(rho.IntegrateGrammage(trajectory, l) ==
-            inhMedium.IntegratedGrammage(trajectory, l));
-    REQUIRE(rho.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
-            inhMedium.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
-  }
-}

Framework/CMakeLists.txt

-
-add_subdirectory (Utilities)
-add_subdirectory (Units)
-add_subdirectory (Geometry)
-add_subdirectory (Particles)
-add_subdirectory (Logging)
-add_subdirectory (StackInterface)
-add_subdirectory (ProcessSequence)
-add_subdirectory (Cascade)
-add_subdirectory (Random)

Framework/Cascade/CMakeLists.txt

-
-# namespace of library -> location of header files
-set (
-  CORSIKAcascade_NAMESPACE
-  corsika/cascade
-  )
-
-# header files of this library
-set (
-  CORSIKAcascade_HEADERS
-  Cascade.h
-  testCascade.h
-  )
-
-add_library (CORSIKAcascade INTERFACE)
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
-
-# include directive for upstream code
-target_include_directories (
-  CORSIKAcascade
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  )
-
-# install library
-install (
-  FILES ${CORSIKAcascade_HEADERS}
-  DESTINATION include/${CORSIKAcascade_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
-add_executable (
-  testCascade
-  testCascade.cc
-  )
-
-target_link_libraries (
-  testCascade
-  #  CORSIKAutls
-  CORSIKArandom
-  ProcessSibyll
-  CORSIKAcascade
-  ProcessStackInspector
-  ProcessTrackingLine
-  ProcessNullModel
-  CORSIKAstackinterface
-  CORSIKAprocesses
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAprocesssequence
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-CORSIKA_ADD_TEST(testCascade)

Framework/Cascade/Cascade.dox

-/**
-
-
-  Here are have to explain the corsika::cascade::Cascade class and its
-  functionality.
-  
-
-*/
\ No newline at end of file

Framework/Cascade/Cascade.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_corsika_cascade_Cascade_h_
-#define _include_corsika_cascade_Cascade_h_
-
-#include <corsika/environment/Environment.h>
-#include <corsika/process/ProcessReturn.h>
-#include <corsika/random/ExponentialDistribution.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/random/UniformRealDistribution.h>
-#include <corsika/stack/SecondaryView.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/setup/SetupTrajectory.h>
-
-/*  see Issue 161, we need to include SetupStack only because we need
-    to globally define StackView. This is clearly not nice and should
-    be changed, when possible. It might be that StackView needs to be
-    templated in Cascade, but this would be even worse... so we don't
-    do that until it is really needed.
- */
-#include <corsika/setup/SetupStack.h>
-
-#include <cassert>
-#include <cmath>
-#include <iostream>
-#include <limits>
-#include <type_traits>
-
-#include <boost/type_index.hpp>
-using boost::typeindex::type_id_with_cvr;
-
-/**
- * The cascade namespace assembles all objects needed to simulate full particles cascades.
- */
-
-namespace corsika::cascade {
-
-  /**
-   * \class Cascade
-   *
-   * The Cascade class is constructed from template arguments making
-   * it very versatile. Via the template arguments physics models are
-   * plugged into the cascade simulation.
-   *
-   * <b>TTracking</b> must be a class according to the
-   * TrackingInterface providing the functions:
-   * <code>auto GetTrack(Particle const& p)</auto>,
-   * with the return type <code>geometry::Trajectory<corsika::geometry::Line>
-   * </code>
-   *
-   * <b>TProcessList</b> must be a ProcessSequence.   *
-   * <b>Stack</b> is the storage object for particle data, i.e. with
-   * Particle class type <code>Stack::ParticleType</code>
-   *
-   *
-   */
-
-  template <typename TTracking, typename TProcessList, typename TStack,
-            /*
-              TStackView is needed as template parameter because of issue 161 and the
-              inability of clang to understand "MakeView" so far.
-             */
-            typename TStackView = corsika::setup::StackView>
-  class Cascade {
-    using Particle = typename TStack::ParticleType;
-    using VolumeTreeNode =
-        std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
-    using MediumInterface = typename VolumeTreeNode::IModelProperties;
-
-    // we only want fully configured objects
-    Cascade() = delete;
-
-  public:
-    /**
-     * Cascade class cannot be default constructed, but needs a valid
-     * list of physics processes for configuration at construct time.
-     */
-    Cascade(corsika::environment::Environment<MediumInterface> const& env, TTracking& tr,
-            TProcessList& pl, TStack& stack)
-        : fEnvironment(env)
-        , fTracking(tr)
-        , fProcessSequence(pl)
-        , fStack(stack) {}
-
-    /**
-     * The Init function is called before the actual cascade simulations.
-     * All components of the Cascade simulation must be configured here.
-     */
-    void Init() {
-      fProcessSequence.Init();
-      fStack.Init();
-    }
-
-    /**
-     * set the nodes for all particles on the stack according to their numerical
-     * position
-     */
-    void SetNodes() {
-      std::for_each(fStack.begin(), fStack.end(), [&](auto& p) {
-        auto const* numericalNode =
-            fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
-        p.SetNode(numericalNode);
-      });
-    }
-
-    /**
-     * The Run function is the main simulation loop, which processes
-     * particles from the Stack until the Stack is empty.
-     */
-    void Run() {
-      SetNodes();
-
-      while (!fStack.IsEmpty()) {
-        while (!fStack.IsEmpty()) {
-          auto pNext = fStack.GetNextParticle();
-          Step(pNext);
-          fProcessSequence.DoStack(fStack);
-        }
-        // do cascade equations, which can put new particles on Stack,
-        // thus, the double loop
-        // DoCascadeEquations();
-      }
-    }
-
-  private:
-    /**
-     * The Step function is executed for each particle from the
-     * stack. It will calcualte geometric transport of the particles,
-     * and apply continuous and stochastic processes to it, which may
-     * lead to energy losses, scattering, absorption, decays and the
-     * production of secondary particles.
-     *
-     * New particles produced in one step are subject to further
-     * processing, e.g. thinning, etc.
-     */
-    void Step(Particle& vParticle) {
-      using namespace corsika;
-      using namespace corsika::units::si;
-
-      // determine geometric tracking
-      auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
-      [[maybe_unused]] auto const& dummy_nextVol = nextVol;
-
-      // determine combined total interaction length (inverse)
-      InverseGrammageType const total_inv_lambda =
-          fProcessSequence.GetTotalInverseInteractionLength(vParticle);
-
-      // sample random exponential step length in grammage
-      corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
-      GrammageType const next_interact = expDist(fRNG);
-
-      std::cout << "total_inv_lambda=" << total_inv_lambda
-                << ", next_interact=" << next_interact << std::endl;
-
-      auto const* currentLogicalNode = vParticle.GetNode();
-
-      // assert that particle stays outside void Universe if it has no
-      // model properties set
-      assert(currentLogicalNode != &*fEnvironment.GetUniverse() ||
-             fEnvironment.GetUniverse()->HasModelProperties());
-
-      // convert next_step from grammage to length
-      LengthType const distance_interact =
-          currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
-                                                                         next_interact);
-
-      // determine the maximum geometric step length
-      LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
-      std::cout << "distance_max=" << distance_max << std::endl;
-
-      // determine combined total inverse decay time
-      InverseTimeType const total_inv_lifetime =
-          fProcessSequence.GetTotalInverseLifetime(vParticle);
-
-      // sample random exponential decay time
-      corsika::random::ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
-      TimeType const next_decay = expDistDecay(fRNG);
-      std::cout << "total_inv_lifetime=" << total_inv_lifetime
-                << ", next_decay=" << next_decay << std::endl;
-
-      // convert next_decay from time to length [m]
-      LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
-                                        vParticle.GetEnergy() * units::constants::c;
-
-      // take minimum of geometry, interaction, decay for next step
-      auto const min_distance =
-          std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
-
-      std::cout << " move particle by : " << min_distance << std::endl;
-
-      // here the particle is actually moved along the trajectory to new position:
-      // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
-      vParticle.SetPosition(step.PositionFromArclength(min_distance));
-      // .... also update time, momentum, direction, ...
-      vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
-
-      step.LimitEndTo(min_distance);
-
-      // apply all continuous processes on particle + track
-      process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
-
-      if (status == process::EProcessReturn::eParticleAbsorbed) {
-        std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
-                  << vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl;
-        vParticle.Delete();
-        return;
-      }
-
-      std::cout << "sth. happening before geometric limit ? "
-                << ((min_distance < geomMaxLength) ? "yes" : "no") << std::endl;
-
-      /*
-        Create SecondaryView object on Stack. The data container
-        remains untouched and identical, and 'projectil' is identical
-        to 'vParticle' above this line. However,
-        projectil.AddSecondaries populate the SecondaryView, which can
-        then be used afterwards for further processing. Thus: it is
-        important to use projectle (and not vParticle) for Interaction,
-        and Decay!
-       */
-
-      TStackView secondaries(vParticle);
-      [[maybe_unused]] auto projectile = secondaries.GetProjectile();
-
-      if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
-
-        // check whether decay or interaction limits this step the
-        // outcome of decay or interaction MAY be a) new particles in
-        // secondaries, b) the projectile particle deleted (or
-        // changed)
-
-        if (min_distance == distance_interact) {
-          std::cout << "collide" << std::endl;
-
-          InverseGrammageType const current_inv_length =
-              fProcessSequence.GetTotalInverseInteractionLength(vParticle);
-
-          random::UniformRealDistribution<InverseGrammageType> uniDist(
-              current_inv_length);
-          const auto sample_process = uniDist(fRNG);
-          InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
-          fProcessSequence.SelectInteraction(vParticle, projectile, sample_process,
-                                             inv_lambda_count);
-        } else if (min_distance == distance_decay) {
-          std::cout << "decay" << std::endl;
-          InverseTimeType const actual_decay_time =
-              fProcessSequence.GetTotalInverseLifetime(vParticle);
-
-          random::UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
-          const auto sample_process = uniDist(fRNG);
-          InverseTimeType inv_decay_count = 0 / second;
-          fProcessSequence.SelectDecay(vParticle, projectile, sample_process,
-                                       inv_decay_count);
-        } else { // step-length limitation within volume
-          std::cout << "step-length limitation" << std::endl;
-        }
-
-        fProcessSequence.DoSecondaries(secondaries);
-        vParticle.Delete(); // todo: this should be reviewed. Where
-                            // exactly are particles best deleted, and
-                            // where they should NOT be
-                            // deleted... maybe Delete function should
-                            // be "protected" and not accessible to physics
-
-        auto const assertion = [&] {
-          auto const* numericalNodeAfterStep =
-              fEnvironment.GetUniverse()->GetContainingNode(vParticle.GetPosition());
-          return numericalNodeAfterStep == currentLogicalNode;
-        };
-
-        assert(assertion()); // numerical and logical nodes don't match
-      } else {               // boundary crossing, step is limited by volume boundary
-        std::cout << "boundary crossing! next node = " << nextVol << std::endl;
-        vParticle.SetNode(nextVol);
-        // DoBoundary may delete the particle (or not)
-        fProcessSequence.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
-      }
-    }
-
-  private:
-    corsika::environment::Environment<MediumInterface> const& fEnvironment;
-    TTracking& fTracking;
-    TProcessList& fProcessSequence;
-    TStack& fStack;
-    corsika::random::RNG& fRNG =
-        corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
-  }; // namespace corsika::cascade
-
-} // namespace corsika::cascade
-
-#endif

Framework/Cascade/testCascade.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/cascade/testCascade.h>
-
-#include <corsika/cascade/Cascade.h>
-
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/null_model/NullModel.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::units::si;
-using namespace corsika::geometry;
-
-#include <iostream>
-#include <limits>
-using namespace std;
-
-auto MakeDummyEnv() {
-  TestEnvironmentType env; // dummy environment
-  auto& universe = *(env.GetUniverse());
-
-  auto theMedium = TestEnvironmentType::CreateNode<Sphere>(
-      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
-      100_km * std::numeric_limits<double>::infinity());
-
-  using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
-  theMedium->SetModelProperties<MyHomogeneousModel>(
-      1_g / (1_cm * 1_cm * 1_cm),
-      environment::NuclearComposition(
-          std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
-
-  universe.AddChild(std::move(theMedium));
-
-  return env;
-}
-
-class ProcessSplit : public process::InteractionProcess<ProcessSplit> {
-
-  int fCalls = 0;
-  GrammageType fX0;
-
-public:
-  ProcessSplit(GrammageType const X0)
-      : fX0(X0) {}
-
-  template <typename Particle>
-  corsika::units::si::GrammageType GetInteractionLength(Particle const&) const {
-    return fX0;
-  }
-
-  template <typename TProjectile>
-  corsika::process::EProcessReturn DoInteraction(TProjectile& vP) {
-    fCalls++;
-    const HEPEnergyType E = vP.GetEnergy();
-    vP.AddSecondary(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            vP.GetPID(), E / 2, vP.GetMomentum(), vP.GetPosition(), vP.GetTime()});
-    vP.AddSecondary(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            vP.GetPID(), E / 2, vP.GetMomentum(), vP.GetPosition(), vP.GetTime()});
-    return EProcessReturn::eInteracted;
-  }
-
-  void Init() { fCalls = 0; }
-
-  int GetCalls() const { return fCalls; }
-};
-
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
-  int fCount = 0;
-  int fCalls = 0;
-  HEPEnergyType fEcrit;
-
-public:
-  ProcessCut(HEPEnergyType e)
-      : fEcrit(e) {}
-
-  template <typename TStack>
-  EProcessReturn DoSecondaries(TStack& vS) {
-    fCalls++;
-    auto p = vS.begin();
-    while (p != vS.end()) {
-      HEPEnergyType E = p.GetEnergy();
-      if (E < fEcrit) {
-        p.Delete();
-        fCount++;
-      } else {
-        ++p; // next particle
-      }
-    }
-    cout << "ProcessCut::DoSecondaries size=" << vS.GetSize() << " count=" << fCount
-         << endl;
-    return EProcessReturn::eOk;
-  }
-
-  void Init() {
-    fCalls = 0;
-    fCount = 0;
-  }
-
-  int GetCount() const { return fCount; }
-  int GetCalls() const { return fCalls; }
-};
-
-TEST_CASE("Cascade", "[Cascade]") {
-  random::RNGManager& rmng = random::RNGManager::GetInstance();
-  rmng.RegisterRandomStream("cascade");
-
-  auto env = MakeDummyEnv();
-  tracking_line::TrackingLine tracking;
-
-  stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true);
-  null_model::NullModel nullModel;
-
-  const GrammageType X0 = 20_g / square(1_cm);
-  const HEPEnergyType Ecrit = 85_MeV;
-  ProcessSplit split(X0);
-  ProcessCut cut(Ecrit);
-  auto sequence = nullModel << stackInspect << split << cut;
-  TestCascadeStack stack;
-
-  cascade::Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack,
-                   TestCascadeStackView>
-      EAS(env, tracking, sequence, stack);
-  CoordinateSystem const& rootCS =
-      RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-  stack.Clear();
-  HEPEnergyType E0 = 100_GeV;
-  stack.AddParticle(
-      std::tuple<particles::Code, units::si::HEPEnergyType,
-                 corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-          particles::Code::Electron, E0,
-          corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
-          Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
-  EAS.Init();
-  EAS.Run();
-
-  CHECK(cut.GetCount() == 2048);
-  CHECK(cut.GetCalls() == 2047);
-  CHECK(split.GetCalls() == 2047);
-}

Framework/Cascade/testCascade.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _Framework_Cascade_testCascade_h
-#define _Framework_Cascade_testCascade_h
-
-#include <corsika/environment/Environment.h>
-#include <corsika/setup/SetupStack.h>
-
-using TestEnvironmentType =
-    corsika::environment::Environment<corsika::environment::IMediumModel>;
-
-template <typename T>
-using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>;
-
-// 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 TestCascadeStack = corsika::stack::CombinedStack<
-    typename corsika::setup::detail::ParticleDataStack::StackImpl,
-    GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
-
-/*
-  See also Issue 161
-*/
-#if defined(__clang__)
-using TestCascadeStackView =
-    corsika::stack::SecondaryView<typename TestCascadeStack::StackImpl,
-                                  StackWithGeometryInterface>;
-#elif defined(__GNUC__) || defined(__GNUG__)
-using TestCascadeStackView = corsika::stack::MakeView<TestCascadeStack>::type;
-#endif
-
-#endif

Framework/Geometry/BaseTrajectory.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_BASETRAJECTORY_H
-#define _include_BASETRAJECTORY_H
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <string>
-
-namespace corsika::geometry {
-
-  /*!
-   * Interface / base class for trajectories.
-   */
-  class BaseTrajectory {
-
-    BaseTrajectory() = delete;
-
-  public:
-    BaseTrajectory(corsika::units::si::TimeType start, corsika::units::si::TimeType end)
-        : fTStart(start)
-        , fTEnd(end) {}
-
-    //!< for \f$ t = 0 \f$, the starting Point shall be returned.
-    virtual Point GetPosition(corsika::units::si::TimeType) const = 0;
-
-    //!< the Point is return from u=0 (start) to u=1 (end)
-    virtual Point GetPosition(double u) const = 0;
-
-    /*!
-     * returns the length between two points of the trajectory
-     * parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
-     */
-
-    virtual corsika::units::si::TimeType TimeFromArclength(
-        corsika::units::si::LengthType) const = 0;
-
-    virtual LengthType ArcLength(corsika::units::si::TimeType t1,
-                                 corsika::units::si::TimeType t2) const = 0;
-
-    virtual corsika::units::si::TimeType GetDuration(
-        corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
-      return t2 - t1;
-    }
-
-    virtual Point GetEndpoint() const { return GetPosition(fTEnd); }
-    virtual Point GetStartpoint() const { return GetPosition(fTStart); }
-
-  protected:
-    corsika::units::si::TimeType const fTStart, fTEnd;
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/BaseVector.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_BASEVECTOR_H_
-#define _include_BASEVECTOR_H_
-
-#include <corsika/geometry/CoordinateSystem.h>
-#include <corsika/geometry/QuantityVector.h>
-
-namespace corsika::geometry {
-
-  /*!
-   * Common base class for Vector and Point. Currently it does basically nothing.
-   */
-
-  template <typename dim>
-  class BaseVector {
-  protected:
-    QuantityVector<dim> qVector;
-    CoordinateSystem const* cs;
-
-  public:
-    BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
-        : qVector(pQVector)
-        , cs(&pCS) {}
-
-    auto const& GetCoordinateSystem() const { return *cs; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/CMakeLists.txt

-
-set (
-  GEOMETRY_SOURCES
-  CoordinateSystem.cc
-  )
-
-set (
-  GEOMETRY_HEADERS
-  Vector.h
-  Point.h
-  Line.h
-  Sphere.h
-  Plane.h
-  Volume.h
-  CoordinateSystem.h
-  RootCoordinateSystem.h
-  Helix.h
-  BaseVector.h
-  QuantityVector.h
-  Trajectory.h
-  FourVector.h
-  )
-
-set (
-  GEOMETRY_NAMESPACE
-  corsika/geometry
-  )
-
-add_library (CORSIKAgeometry STATIC ${GEOMETRY_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAgeometry ${GEOMETRY_NAMESPACE} ${GEOMETRY_HEADERS})
-
-set_target_properties (
-  CORSIKAgeometry
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  PUBLIC_HEADER "${GEOMETRY_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAutilities
-  )
-
-target_include_directories (
-  CORSIKAgeometry
-  SYSTEM
-  PUBLIC    ${EIGEN3_INCLUDE_DIR}
-  INTERFACE ${EIGEN3_INCLUDE_DIR}
-  )
-
-target_include_directories (
-  CORSIKAgeometry
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS CORSIKAgeometry
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${GEOMETRY_NAMESPACE}
-  )
-
-# --------------------
-# code unit testing
-add_executable (testGeometry testGeometry.cc)
-target_link_libraries (
-  testGeometry
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-CORSIKA_ADD_TEST(testGeometry)
-
-add_executable (testFourVector testFourVector.cc) 
-target_link_libraries (
-  testFourVector
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-CORSIKA_ADD_TEST(testFourVector)

Framework/Geometry/CoordinateSystem.cc

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/geometry/CoordinateSystem.h>
-#include <stdexcept>
-
-using namespace corsika::geometry;
-
-/**
- * returns the transformation matrix necessary to transform primitives with coordinates
- * in \a pFrom to \a pTo, e.g.
- * \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
- * (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
- * the indicated CoordinateSystem).
- */
-EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom,
-                                                   CoordinateSystem const& pTo) {
-  CoordinateSystem const* a{&pFrom};
-  CoordinateSystem const* b{&pTo};
-  CoordinateSystem const* commonBase{nullptr};
-
-  while (a != b && b != nullptr) {
-    a = &pFrom;
-
-    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::runtime_error("no connection between coordinate systems found!");
-  }
-
-  EigenTransform t = EigenTransform::Identity();
-  auto* p = &pFrom;
-
-  while (p != commonBase) {
-    t = p->GetTransform() * t;
-    p = p->GetReference();
-  }
-
-  p = &pTo;
-
-  while (p != commonBase) {
-    t = t * p->GetTransform().inverse(Eigen::TransformTraits::Isometry);
-    p = p->GetReference();
-  }
-
-  return t;
-}

Framework/Geometry/CoordinateSystem.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_COORDINATESYSTEM_H_
-#define _include_COORDINATESYSTEM_H_
-
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/sgn.h>
-#include <Eigen/Dense>
-#include <stdexcept>
-
-typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
-typedef Eigen::Translation<double, 3> EigenTranslation;
-
-namespace corsika::geometry {
-
-  class RootCoordinateSystem;
-  template <typename T>
-  class Vector;
-
-  using corsika::units::si::length_d;
-
-  class CoordinateSystem {
-    CoordinateSystem const* reference = nullptr;
-    EigenTransform transf;
-
-    CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
-        : reference(&reference)
-        , transf(transf) {}
-
-    CoordinateSystem()
-        : // for creating the root CS
-        transf(EigenTransform::Identity()) {}
-
-  protected:
-    static auto CreateCS() { return CoordinateSystem(); }
-    friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
-                                                    /// create ONE unique root CS
-
-  public:
-    static EigenTransform GetTransformation(CoordinateSystem const& c1,
-                                            CoordinateSystem const& c2);
-
-    auto& operator=(const CoordinateSystem& pCS) {
-      reference = pCS.reference;
-      transf = pCS.transf;
-      return *this;
-    }
-
-    auto translate(QuantityVector<length_d> vector) const {
-      EigenTransform const translation{EigenTranslation(vector.eVector)};
-
-      return CoordinateSystem(*this, translation);
-    }
-
-    template <typename TDim>
-    auto RotateToZ(Vector<TDim> vVec) const {
-      auto const a = vVec.normalized().GetComponents(*this).eVector;
-      auto const a1 = a(0), a2 = a(1);
-
-      auto const s = utl::sgn(a(2));
-      auto const c = 1 / (1 + s * a(2));
-
-      Eigen::Matrix3d A, B;
-
-      if (s > 0) {
-        A << 1, 0, -a1,                     // comment to prevent clang-format
-            0, 1, -a2,                      // .
-            a1, a2, 1;                      // .
-        B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-            -a1 * a2 * c, -a2 * a2 * c, 0,  // .
-            0, 0, -(a1 * a1 + a2 * a2) * c; // .
-
-      } else {
-        A << 1, 0, a1,                      // .
-            0, -1, -a2,                     // .
-            a1, a2, -1;                     // .
-        B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-            +a1 * a2 * c, +a2 * a2 * c, 0,  // .
-            0, 0, (a1 * a1 + a2 * a2) * c;  // .
-      }
-
-      return CoordinateSystem(*this, EigenTransform(A + B));
-    }
-
-    template <typename TDim>
-    auto rotate(QuantityVector<TDim> axis, double angle) const {
-      if (axis.eVector.isZero()) {
-        throw std::runtime_error("null-vector given as axis parameter");
-      }
-
-      EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
-
-      return CoordinateSystem(*this, rotation);
-    }
-
-    template <typename TDim>
-    auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
-                            QuantityVector<TDim> axis, double angle) {
-      if (axis.eVector.isZero()) {
-        throw std::runtime_error("null-vector given as axis parameter");
-      }
-
-      EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
-                                  EigenTranslation(translation.eVector)};
-
-      return CoordinateSystem(*this, transf);
-    }
-
-    auto const* GetReference() const { return reference; }
-
-    auto const& GetTransform() const { return transf; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/FourVector.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_corsika_framework_geometry_fourvector_h_
-#define _include_corsika_framework_geometry_fourvector_h_
-
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <iostream>
-#include <type_traits>
-
-namespace corsika::geometry {
-
-  /**
-     FourVector supports "full" units, e.g. E in [GeV/c] and p in [GeV],
-     or also t in [s] and r in [m], etc.
-
-     However, for HEP applications it is also possible to use E and p
-     both in [GeV].
-
-     The FourVector can return NormSqr and Norm, whereas Norm is
-     sqrt(abs(NormSqr)). The physical units are always calculated and
-     returned properly.
-
-     FourVector can also return if it is TimeLike, SpaceLike or PhotonLike.
-
-     When a FourVector is initialized with a lvalue reference, this is
-     also used for the internal storage, which should lead to complete
-     disappearance of the FourVector class during optimization.
-   */
-
-  template <typename TimeType, typename SpaceVecType>
-  class FourVector {
-
-  public:
-    using SpaceType = typename std::decay<SpaceVecType>::type::Quantity;
-
-    //! check the types and the physical units here:
-    static_assert(
-        std::is_same<typename std::decay<TimeType>::type, SpaceType>::value ||
-            std::is_same<typename std::decay<TimeType>::type,
-                         decltype(std::declval<SpaceType>() / corsika::units::si::meter *
-                                  corsika::units::si::second)>::value,
-        "Units of time-like and space-like coordinates must either be idential "
-        "(e.g. GeV) or [E/c]=[p]");
-
-  public:
-    FourVector(const TimeType& eT, const SpaceVecType& eS)
-        : fTimeLike(eT)
-        , fSpaceLike(eS) {}
-
-    TimeType GetTimeLikeComponent() const { return fTimeLike; }
-    SpaceVecType& GetSpaceLikeComponents() { return fSpaceLike; }
-    const SpaceVecType& GetSpaceLikeComponents() const { return fSpaceLike; }
-
-    auto GetNormSqr() const { return GetTimeSquared() - fSpaceLike.squaredNorm(); }
-
-    SpaceType GetNorm() const { return sqrt(abs(GetNormSqr())); }
-
-    bool IsTimelike() const {
-      return GetTimeSquared() < fSpaceLike.squaredNorm();
-    } //! Norm2 < 0
-
-    bool IsSpacelike() const {
-      return GetTimeSquared() > fSpaceLike.squaredNorm();
-    } //! Norm2 > 0
-
-    /* this is not numerically stable
-    bool IsPhotonlike() const {
-      return GetTimeSquared() == fSpaceLike.squaredNorm();
-    } //! Norm2 == 0
-    */
-
-    FourVector& operator+=(const FourVector& b) {
-      fTimeLike += b.fTimeLike;
-      fSpaceLike += b.fSpaceLike;
-      return *this;
-    }
-
-    FourVector& operator-=(const FourVector& b) {
-      fTimeLike -= b.fTimeLike;
-      fSpaceLike -= b.fSpaceLike;
-      return *this;
-    }
-
-    FourVector& operator*=(const double b) {
-      fTimeLike *= b;
-      fSpaceLike *= b;
-      return *this;
-    }
-
-    FourVector& operator/=(const double b) {
-      fTimeLike /= b;
-      fSpaceLike.GetComponents() /= b; // TODO: WHY IS THIS??????
-      return *this;
-    }
-
-    FourVector& operator/(const double b) {
-      *this /= b;
-      return *this;
-    }
-
-    /**
-       Note that the product between two 4-vectors assumes that you use
-       the same "c" convention for both. Only the LHS vector is checked
-       for this. You cannot mix different conventions due to
-       unit-checking.
-     */
-    SpaceType operator*(const FourVector& b) {
-      if constexpr (std::is_same<typename std::decay<TimeType>::type,
-                                 decltype(std::declval<SpaceType>() /
-                                          corsika::units::si::meter *
-                                          corsika::units::si::second)>::value)
-        return fTimeLike * b.fTimeLike *
-                   (corsika::units::constants::c * corsika::units::constants::c) -
-               fSpaceLike.norm();
-      else
-        return fTimeLike * fTimeLike - fSpaceLike.norm();
-    }
-
-  private:
-    /**
-       This function is automatically compiled to use of ignore the
-       extra factor of "c" for the time-like quantity
-     */
-    auto GetTimeSquared() const {
-      if constexpr (std::is_same<typename std::decay<TimeType>::type,
-                                 decltype(std::declval<SpaceType>() /
-                                          corsika::units::si::meter *
-                                          corsika::units::si::second)>::value)
-        return fTimeLike * fTimeLike *
-               (corsika::units::constants::c * corsika::units::constants::c);
-      else
-        return fTimeLike * fTimeLike;
-    }
-
-  protected:
-    //! the data members
-    TimeType fTimeLike;
-    SpaceVecType fSpaceLike;
-
-    //! the friends: math operators
-    template <typename T, typename U>
-    friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
-    operator+(const FourVector<T, U>&, const FourVector<T, U>&);
-
-    template <typename T, typename U>
-    friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
-    operator-(const FourVector<T, U>&, const FourVector<T, U>&);
-
-    template <typename T, typename U>
-    friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
-    operator*(const FourVector<T, U>&, const double);
-
-    template <typename T, typename U>
-    friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
-    operator/(const FourVector<T, U>&, const double);
-  };
-
-  /**
-      The math operator+
-   */
-  template <typename TimeType, typename SpaceVecType>
-  inline FourVector<typename std::decay<TimeType>::type,
-                    typename std::decay<SpaceVecType>::type>
-  operator+(const FourVector<TimeType, SpaceVecType>& a,
-            const FourVector<TimeType, SpaceVecType>& b) {
-    return FourVector<typename std::decay<TimeType>::type,
-                      typename std::decay<SpaceVecType>::type>(
-        a.fTimeLike + b.fTimeLike, a.fSpaceLike + b.fSpaceLike);
-  }
-
-  /**
-     The math operator-
-  */
-  template <typename TimeType, typename SpaceVecType>
-  inline FourVector<typename std::decay<TimeType>::type,
-                    typename std::decay<SpaceVecType>::type>
-  operator-(const FourVector<TimeType, SpaceVecType>& a,
-            const FourVector<TimeType, SpaceVecType>& b) {
-    return FourVector<typename std::decay<TimeType>::type,
-                      typename std::decay<SpaceVecType>::type>(
-        a.fTimeLike - b.fTimeLike, a.fSpaceLike - b.fSpaceLike);
-  }
-
-  /**
-     The math operator*
-  */
-  template <typename TimeType, typename SpaceVecType>
-  inline FourVector<typename std::decay<TimeType>::type,
-                    typename std::decay<SpaceVecType>::type>
-  operator*(const FourVector<TimeType, SpaceVecType>& a, const double b) {
-    return FourVector<typename std::decay<TimeType>::type,
-                      typename std::decay<SpaceVecType>::type>(a.fTimeLike * b,
-                                                               a.fSpaceLike * b);
-  }
-
-  /**
-      The math operator/
-   */
-  template <typename TimeType, typename SpaceVecType>
-  inline FourVector<typename std::decay<TimeType>::type,
-                    typename std::decay<SpaceVecType>::type>
-  operator/(const FourVector<TimeType, SpaceVecType>& a, const double b) {
-    return FourVector<typename std::decay<TimeType>::type,
-                      typename std::decay<SpaceVecType>::type>(a.fTimeLike / b,
-                                                               a.fSpaceLike / b);
-  }
-
-} // namespace corsika::geometry
-
-#endif