diff --git a/corsika/detail/media/Environment.inl b/corsika/detail/media/Environment.inl
index cc5a945519bed11dffccfacf9084a1008ad2e074..434d9c77cfa443a2c42969501f96964a7b2b69a7 100644
--- a/corsika/detail/media/Environment.inl
+++ b/corsika/detail/media/Environment.inl
@@ -14,26 +14,37 @@
namespace corsika {
- template <typename IEnvironmentModel>
- auto& Environment<IEnvironmentModel>::GetUniverse() {
- return fUniverse; }
-
- template <typename IEnvironmentModel>
- auto const& Environment<IEnvironmentModel>::GetUniverse() const { return fUniverse; }
-
- template <typename IEnvironmentModel>
- auto const& Environment<IEnvironmentModel>::GetCoordinateSystem() const { return fCoordinateSystem; }
-
- // factory method for creation of VolumeTreeNodes
- template <typename IEnvironmentModel>
- template <class TVolumeType, typename... TVolumeArgs>
- auto Environment<IEnvironmentModel>::CreateNode(TVolumeArgs&&... args) {
- static_assert(std::is_base_of_v<corsika::Volume, TVolumeType>,
- "unusable type provided, needs to be derived from "
- "\"corsika::Volume\"");
-
- return std::make_unique<BaseNodeType>(
- std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
- }
-
-}
\ No newline at end of file
+ template <typename IEnvironmentModel>
+ Environment<IEnvironmentModel>::Environment()
+ : coordinateSystem_{RootCoordinateSystem::getInstance().GetRootCoordinateSystem()}
+ , universe_(std::make_unique<BaseNodeType>(
+ std::make_unique<Universe>(coordinateSystem_))) {}
+
+ template <typename IEnvironmentModel>
+ typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr& Environment<IEnvironmentModel>::getUniverse() {
+ return universe_;
+ }
+
+ template <typename IEnvironmentModel>
+ typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr const& Environment<IEnvironmentModel>::getUniverse() const {
+ return universe_;
+ }
+
+ template <typename IEnvironmentModel>
+ CoordinateSystem const& Environment<IEnvironmentModel>::getCoordinateSystem() const {
+ return coordinateSystem_;
+ }
+
+ // factory method for creation of VolumeTreeNodes
+ template <typename IEnvironmentModel>
+ template <class TVolumeType, typename... TVolumeArgs>
+ std::unique_ptr< VolumeTreeNode<IEnvironmentModel> > Environment<IEnvironmentModel>::CreateNode(TVolumeArgs&&... args) {
+ static_assert(std::is_base_of_v<corsika::Volume, TVolumeType>,
+ "unusable type provided, needs to be derived from "
+ "\"corsika::Volume\"");
+
+ return std::make_unique<BaseNodeType>(
+ std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
+ }
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
index c3b4c7fe64ec803a295a22b379b21716b2c5a3c5..6dc5519d4927917a1ff6ddf9921fc78fe619d2be 100644
--- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
+++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
@@ -70,7 +70,7 @@ namespace corsika {
}
void LayeredSphericalAtmosphereBuilder::assemble(Environment<IMediumModel>& env) {
- auto& universe = env.GetUniverse();
+ auto& universe = env.getUniverse();
auto* outmost = universe.get();
while (!layers_.empty()) {
diff --git a/corsika/detail/media/NuclearComposition.inl b/corsika/detail/media/NuclearComposition.inl
new file mode 100644
index 0000000000000000000000000000000000000000..ba56e6ea3569375f3ac75a23a38363d36f82e341
--- /dev/null
+++ b/corsika/detail/media/NuclearComposition.inl
@@ -0,0 +1,144 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <cassert>
+#include <functional>
+#include <numeric>
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+namespace corsika {
+
+ template <class AConstIterator, class BConstIterator>
+ NuclearComposition::WeightProviderIterator<
+ AConstIterator, BConstIterator>::WeightProviderIterator(AConstIterator a,
+ BConstIterator b)
+ : aIter_(a)
+ , bIter_(b) {}
+
+ template <class AConstIterator, class BConstIterator>
+ typename NuclearComposition::WeightProviderIterator<AConstIterator,
+ BConstIterator>::value_type
+ NuclearComposition::WeightProviderIterator<AConstIterator, BConstIterator>::operator*()
+ const {
+ return ((*aIter_) * (*bIter_)).magnitude();
+ }
+
+ template <class AConstIterator, class BConstIterator>
+ NuclearComposition::WeightProviderIterator<AConstIterator, BConstIterator>&
+ NuclearComposition::WeightProviderIterator<AConstIterator,
+ BConstIterator>::operator++() { // prefix ++
+ ++aIter_;
+ ++bIter_;
+ return *this;
+ }
+
+ template <class AConstIterator, class BConstIterator>
+ auto
+ NuclearComposition::WeightProviderIterator<AConstIterator, BConstIterator>::operator==(
+ WeightProviderIterator other) {
+ return aIter_ == other.aIter_;
+ }
+
+ template <class AConstIterator, class BConstIterator>
+ auto
+ NuclearComposition::WeightProviderIterator<AConstIterator, BConstIterator>::operator!=(
+ WeightProviderIterator other) {
+ return !(*this == other);
+ }
+
+ NuclearComposition::NuclearComposition(std::vector<corsika::Code> pComponents,
+ std::vector<float> pFractions)
+ : numberFractions_(pFractions)
+ , components_(pComponents)
+ , avgMassNumber_(std::inner_product(
+ pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
+ std::plus<double>(), [](auto const compID, auto const fraction) -> double {
+ if (IsNucleus(compID)) {
+ return GetNucleusA(compID) * fraction;
+ } else {
+ return GetMass(compID) / units::si::ConvertSIToHEP(constants::u) * 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");
+ }
+ updateHash();
+ }
+
+ template <typename TFunction>
+ auto NuclearComposition::WeightedSum(TFunction func) const {
+ using ResultQuantity = decltype(func(*components_.cbegin()));
+
+ 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(
+ components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+ ResultQuantity::zero(), // .zero() is defined for quantity types only
+ std::plus<ResultQuantity>(), prod);
+ } else {
+ return std::inner_product(
+ components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+ ResultQuantity(0), // in other cases we have to use a bare 0
+ std::plus<ResultQuantity>(), prod);
+ }
+ }
+
+ auto NuclearComposition::size() const { return numberFractions_.size(); }
+
+ auto const& NuclearComposition::GetFractions() const { return numberFractions_; }
+ auto const& NuclearComposition::GetComponents() const { return components_; }
+ auto const NuclearComposition::GetAverageMassNumber() const { return avgMassNumber_; }
+
+ template <class TRNG>
+ corsika::Code NuclearComposition::SampleTarget(
+ std::vector<units::si::CrossSectionType> const& sigma,
+ TRNG& randomStream) const {
+ using namespace units::si;
+
+ assert(sigma.size() == numberFractions_.size());
+
+ std::discrete_distribution channelDist(
+ WeightProviderIterator<decltype(numberFractions_.begin()),
+ decltype(sigma.begin())>(numberFractions_.begin(),
+ sigma.begin()),
+ WeightProviderIterator<decltype(numberFractions_.begin()), decltype(sigma.end())>(
+ numberFractions_.end(), sigma.end()));
+
+ auto const iChannel = channelDist(randomStream);
+ return components_[iChannel];
+ }
+
+ // Note: when this class ever modifies its internal data, the hash
+ // must be updated, too!
+ size_t NuclearComposition::hash() const { return hash_; }
+
+ void NuclearComposition::updateHash() {
+ std::vector<std::size_t> hashes;
+ for (float ifrac : GetFractions()) hashes.push_back(std::hash<float>{}(ifrac));
+ for (corsika::Code icode : GetComponents())
+ hashes.push_back(std::hash<int>{}(static_cast<int>(icode)));
+ std::size_t h = std::hash<double>{}(GetAverageMassNumber());
+ for (std::size_t ih : hashes) h = h ^ (ih << 1);
+ hash_ = h;
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/media/ShowerAxis.inl b/corsika/detail/media/ShowerAxis.inl
new file mode 100644
index 0000000000000000000000000000000000000000..2f5691794c6496878ababeeefdfee5eec4a8a68a
--- /dev/null
+++ b/corsika/detail/media/ShowerAxis.inl
@@ -0,0 +1,107 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/logging/Logging.hpp>
+
+#include <string>
+
+namespace corsika {
+
+ template <typename TEnvModel>
+ ShowerAxis::ShowerAxis(Point const& pStart,
+ corsika::Vector<units::si::length_d> length,
+ Environment<TEnvModel> const& env,
+ int steps)
+ : pointStart_(pStart)
+ , length_(length)
+ , max_length_(length_.norm())
+ , steplength_(max_length_ / steps)
+ , axis_normalized_(length / max_length_)
+ , X_(steps + 1) {
+ auto const* const universe = env.getUniverse().get();
+
+ auto rho = [pStart, length, universe](double x) {
+ auto const p = pStart + length * x;
+ auto const* node = universe->GetContainingNode(p);
+ return node->GetModelProperties().getMassDensity(p).magnitude();
+ };
+
+ double error;
+ int k = 0;
+ X_[0] = units::si::GrammageType::zero();
+ auto sum = units::si::GrammageType::zero();
+
+ for (int i = 1; i <= steps; ++i) {
+ auto const x_prev = (i - 1.) / steps;
+ auto const d_prev = max_length_ * x_prev;
+ auto const x = double(i) / steps;
+ auto const r = boost::math::quadrature::gauss_kronrod<double, 15>::integrate(
+ rho, x_prev, x, 15, 1e-9, &error);
+ auto const result =
+ units::si::MassDensityType(phys::units::detail::magnitude_tag, r) * max_length_;
+
+ sum += result;
+ X_[i] = sum;
+
+ for (; sum > k * X_binning_; ++k) {
+ d_.emplace_back(d_prev + k * X_binning_ * steplength_ / result);
+ }
+ }
+
+ assert(std::is_sorted(X_.cbegin(), X_.cend()));
+ assert(std::is_sorted(d_.cbegin(), d_.cend()));
+ }
+
+ GrammageType ShowerAxis::X(LengthType l) const {
+ auto const fractionalBin = l / steplength_;
+ int const lower = fractionalBin; // indices of nearest X support points
+ auto const lambda = fractionalBin - lower;
+ decltype(X_.size()) const upper = lower + 1;
+
+ if (lower < 0) {
+ CORSIKA_LOG_ERROR("cannot extrapolate to points behind point of injection l={} m",
+ l / 1_m);
+ throw std::runtime_error("cannot extrapolate to points behind point of injection");
+ }
+
+ if (upper >= X_.size()) {
+ const std::string err =
+ fmt::format("shower axis too short, cannot extrapolate (l / max_length_ = {} )",
+ l / max_length_);
+ CORSIKA_LOG_ERROR(err);
+ throw std::runtime_error(err.c_str());
+ }
+
+ assert(0 <= lambda && lambda <= 1.);
+
+ CORSIKA_LOG_TRACE("ShowerAxis::X l={} m, lower={}, lambda={}, upper={}", l / 1_m, lower,
+ lambda, upper);
+
+ // linear interpolation between X[lower] and X[upper]
+ return X_[upper] * lambda + X_[lower] * (1 - lambda);
+ }
+
+ LengthType ShowerAxis::steplength() const { return steplength_; }
+
+ GrammageType ShowerAxis::maximumX() const { return *X_.rbegin(); }
+
+ GrammageType ShowerAxis::minimumX() const { return GrammageType::zero(); }
+
+ GrammageType ShowerAxis::projectedX(Point const& p) const {
+ auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
+ return X(projectedLength);
+ }
+
+ corsika::Vector<units::si::dimensionless_d> const& ShowerAxis::getDirection() const {
+ return axis_normalized_;
+ }
+
+ Point const& ShowerAxis::getStart() const { return pointStart_; }
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/media/Environment.hpp b/corsika/media/Environment.hpp
index 51cae1847d7ff5fc48d09f17b459bbd740626019..fd2ccc74ccaec39f5f128642f08d82ee2e6e99d5 100644
--- a/corsika/media/Environment.hpp
+++ b/corsika/media/Environment.hpp
@@ -8,12 +8,14 @@
#pragma once
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/VolumeTreeNode.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/media/IMediumModel.hpp>
+
#include <corsika/media/Universe.hpp>
-#include <corsika/media/VolumeTreeNode.hpp>
+
#include <limits>
namespace corsika {
@@ -23,32 +25,22 @@ namespace corsika {
public:
using BaseNodeType = VolumeTreeNode<IEnvironmentModel>;
- Environment()
- : fCoordinateSystem(
- corsika::RootCoordinateSystem::getInstance().GetRootCoordinateSystem())
- , fUniverse(std::make_unique<BaseNodeType>(
- std::make_unique<Universe>(fCoordinateSystem))) {}
-
- // using IEnvironmentModel = corsika::IEnvironmentModel;
+ Environment();
- inline auto& GetUniverse();
+ typename BaseNodeType::VTNUPtr& getUniverse();
+ typename BaseNodeType::VTNUPtr const& getUniverse() const;
- inline auto const& GetUniverse() const;
-
- inline auto const& GetCoordinateSystem() const;
+ CoordinateSystem const& getCoordinateSystem() const;
// factory method for creation of VolumeTreeNodes
template <class TVolumeType, typename... TVolumeArgs>
- static auto CreateNode(TVolumeArgs&&... args);
+ static std::unique_ptr<BaseNodeType> CreateNode(TVolumeArgs&&... args);
private:
- corsika::CoordinateSystem const& fCoordinateSystem;
- typename BaseNodeType::VTNUPtr fUniverse;
+ CoordinateSystem const& coordinateSystem_;
+ typename BaseNodeType::VTNUPtr universe_;
};
- // using SetupBaseNodeType = VolumeTreeNode<corsika::IEnvironmentModel>;
- // using SetupEnvironment = Environment<corsika::IEnvironmentModel>;
-
} // namespace corsika
#include <corsika/detail/media/Environment.inl>
diff --git a/corsika/media/IMagneticFieldModel.hpp b/corsika/media/IMagneticFieldModel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c7c62fab8e3c18c145450a4f0402977eb10bedf1
--- /dev/null
+++ b/corsika/media/IMagneticFieldModel.hpp
@@ -0,0 +1,47 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+ /**
+ * An interface for magnetic field models.
+ *
+ * This is the base interface for magnetic field model mixins.
+ *
+ */
+ template <typename Model>
+ class IMagneticFieldModel : public Model {
+
+ // a type-alias for a magnetic field vector
+ using MagneticFieldVector =
+ corsika::geometry::Vector<corsika::phys::units::si::magnetic_flux_density_d>;
+
+ public:
+ /**
+ * Evaluate the magnetic field at a given location.
+ *
+ * @param point The location to evaluate the field at.
+ * @returns The magnetic field vector at that point.
+ */
+ virtual auto GetMagneticField(corsika::geometry::Point const&) const
+ -> MagneticFieldVector = 0;
+
+ /**
+ * A virtual default destructor.
+ */
+ virtual ~IMagneticFieldModel() = default; // LCOV_EXCL_LINE
+
+ }; // END: class MagneticField
+
+} // namespace corsika::environment
diff --git a/corsika/media/MediumProperties.hpp b/corsika/media/MediumProperties.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f01c128e0bc3f91c7cbf6d91c5b001a528591e3a
--- /dev/null
+++ b/corsika/media/MediumProperties.hpp
@@ -0,0 +1,92 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+namespace corsika {
+
+ /**
+ * Medium types are useful most importantly for effective models
+ * like energy losses. a particular medium (mixture of components)
+ * may have specif properties not reflected by its mixture of
+ * components.
+ **/
+
+ enum class MediumType {
+ Unknown,
+ Element,
+ RadioactiveElement,
+ InorganicCompound,
+ OrganicCompound,
+ Polymer,
+ Mixture,
+ BiologicalDosimetry
+ };
+
+ enum class State { Unknown, Solid, Liquid, Gas, DiatomicGas };
+
+ enum class Medium : int16_t;
+ using MediumIntType = std::underlying_type<Medium>::type;
+
+ /**
+ * \struct MediumData
+ *
+ * Simple object to group together a number of properties
+ *
+ **/
+ struct MediumData {
+ std::string name_;
+ std::string pretty_name_;
+ double weight_;
+ int weight_significant_figure_;
+ int weight_error_last_digit_;
+ double Z_over_A_;
+ double sternheimer_density_;
+ double corrected_density_;
+ State state_;
+ MediumType type_;
+ std::string symbol_;
+ double Ieff_;
+ double Cbar_;
+ double x0_;
+ double x1_;
+ double aa_;
+ double sk_;
+ double dlt0_;
+
+ std::string name() const { return name_; }
+ std::string pretty_name() const { return pretty_name_; }
+ double weight() const { return weight_; }
+ int weight_significant_figure() const { return weight_significant_figure_; }
+ int weight_error_last_digit() const { return weight_error_last_digit_; }
+ double Z_over_A() const { return Z_over_A_; }
+ double sternheimer_density() const { return sternheimer_density_; }
+ double corrected_density() const { return corrected_density_; }
+ State state() const { return state_; }
+ MediumType type() const { return type_; }
+ std::string symbol() const { return symbol_; }
+ double Ieff() const { return Ieff_; }
+ double Cbar() const { return Cbar_; }
+ double x0() const { return x0_; }
+ double x1() const { return x1_; }
+ double aa() const { return aa_; }
+ double sk() const { return sk_; }
+ double dlt0() const { return dlt0_; }
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/media/GeneratedMediaProperties.inc>
+
+namespace corsika::environment {
+
+ constexpr MediumData const& mediumData(Medium const m) {
+ return detail::medium_data[static_cast<MediumIntType>(m)];
+ }
+
+} // namespace corsika::environment
diff --git a/corsika/media/NuclearComposition.hpp b/corsika/media/NuclearComposition.hpp
index f7f4bd5a25391901ff59846afdfbe67798ce6763..feb30c3da97c4932738147ad4ea69f747370588e 100644
--- a/corsika/media/NuclearComposition.hpp
+++ b/corsika/media/NuclearComposition.hpp
@@ -1,5 +1,5 @@
-n/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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
@@ -19,13 +19,19 @@ n/*
#include <vector>
namespace corsika {
+ class NuclearComposition {
+ std::vector<float> const numberFractions_; //!< relative fractions of number density
+ std::vector<corsika::Code> const
+ components_; //!< particle codes of consitutents
+
+ double const avgMassNumber_;
- namespace nuc_comp::detail {
+ std::size_t hash_;
template <class AConstIterator, class BConstIterator>
class WeightProviderIterator {
- AConstIterator fAIter;
- BConstIterator fBIter;
+ AConstIterator aIter_;
+ BConstIterator bIter_;
public:
using value_type = double;
@@ -34,103 +40,43 @@ namespace corsika {
using reference = value_type&;
using difference_type = ptrdiff_t;
- WeightProviderIterator(AConstIterator a, BConstIterator b)
- : fAIter(a)
- , fBIter(b) {}
+ WeightProviderIterator(AConstIterator a, BConstIterator b);
- value_type operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
+ value_type operator*() const;
- WeightProviderIterator& operator++() { // prefix ++
- ++fAIter;
- ++fBIter;
- return *this;
- }
+ WeightProviderIterator& operator++();
- auto operator==(WeightProviderIterator other) { return fAIter == other.fAIter; }
+ auto operator==(WeightProviderIterator other);
- auto operator!=(WeightProviderIterator other) { return !(*this == other); }
+ auto operator!=(WeightProviderIterator other);
};
- } // namespace nuc_comp::detail
-
- class NuclearComposition {
- std::vector<float> const fNumberFractions; //!< relative fractions of number density
- std::vector<corsika::Code> const fComponents; //!< particle codes of consitutents
-
- double const fAvgMassNumber;
-
- /*
- * FIXME: smelling here... nested class why? ...done
- * FIXME: promote it to namespace detail ... done
- */
-
public:
NuclearComposition(std::vector<corsika::Code> pComponents,
- std::vector<float> pFractions)
- : fNumberFractions(pFractions)
- , fComponents(pComponents)
- , fAvgMassNumber(std::inner_product(
- pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
- std::plus<double>(), [](auto const compID, auto const fraction) -> double {
- if (is_nucleus(compID)) {
- return nucleus_A(compID) * fraction;
- } else {
- return mass(compID) / ConvertSIToHEP(constants::u) * 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");
- }
- }
+ std::vector<float> pFractions);
template <typename TFunction>
- auto WeightedSum(TFunction func) const {
- using ResultQuantity = decltype(func(*fComponents.cbegin()));
-
- 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
- std::plus<ResultQuantity>(), prod);
- } else {
- return std::inner_product(
- fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
- ResultQuantity(0), // in other cases we have to use a bare 0
- std::plus<ResultQuantity>(), 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; }
+ auto WeightedSum(TFunction func) const ;
+
+ auto size() const;
+
+ auto const& GetFractions() const;
+ auto const& GetComponents() const ;
+ auto const GetAverageMassNumber() const;
template <class TRNG>
- corsika::Code SampleTarget(std::vector<CrossSectionType> const& sigma,
- TRNG& randomStream) const {
-
- assert(sigma.size() == fNumberFractions.size());
-
- std::discrete_distribution channelDist(
- nuc_comp::detail::WeightProviderIterator<decltype(fNumberFractions.begin()),
- decltype(sigma.begin())>(
- fNumberFractions.begin(), sigma.begin()),
- nuc_comp::detail::WeightProviderIterator<decltype(fNumberFractions.begin()),
- decltype(sigma.end())>(
- fNumberFractions.end(), sigma.end()));
-
- auto const iChannel = channelDist(randomStream);
- return fComponents[iChannel];
- }
+ Code SampleTarget(
+ std::vector<units::si::CrossSectionType> const& sigma,
+ TRNG& randomStream) const;
+
+ // Note: when this class ever modifies its internal data, the hash
+ // must be updated, too!
+ size_t hash() const;
+
+ private:
+ void updateHash();
};
} // namespace corsika
+
+#include <corsika/detail/media/NuclearComposition.inl>
diff --git a/corsika/media/NuclearComposition_old.hpp b/corsika/media/NuclearComposition_old.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..82577f24545c7d4598e5ca09a84989a3e358bd33
--- /dev/null
+++ b/corsika/media/NuclearComposition_old.hpp
@@ -0,0 +1,136 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <cassert>
+#include <functional>
+#include <numeric>
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+namespace corsika {
+
+ namespace nuc_comp::detail {
+
+ 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); }
+ };
+
+ } // namespace nuc_comp::detail
+
+ class NuclearComposition {
+ std::vector<float> const fNumberFractions; //!< relative fractions of number density
+ std::vector<corsika::Code> const fComponents; //!< particle codes of consitutents
+
+ double const fAvgMassNumber;
+
+ /*
+ * FIXME: smelling here... nested class why? ...done
+ * FIXME: promote it to namespace detail ... done
+ */
+
+ public:
+ NuclearComposition(std::vector<corsika::Code> pComponents,
+ std::vector<float> pFractions)
+ : fNumberFractions(pFractions)
+ , fComponents(pComponents)
+ , fAvgMassNumber(std::inner_product(
+ pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
+ std::plus<double>(), [](auto const compID, auto const fraction) -> double {
+ if (IsNucleus(compID)) {
+ return GetNucleusA(compID) * fraction;
+ } else {
+ return GetMass(compID) / ConvertSIToHEP(constants::u) * 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 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
+ std::plus<ResultQuantity>(), prod);
+ } else {
+ return std::inner_product(
+ fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
+ ResultQuantity(0), // in other cases we have to use a bare 0
+ std::plus<ResultQuantity>(), 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::Code SampleTarget(std::vector<CrossSectionType> const& sigma,
+ TRNG& randomStream) const {
+
+ assert(sigma.size() == fNumberFractions.size());
+
+ std::discrete_distribution channelDist(
+ nuc_comp::detail::WeightProviderIterator<decltype(fNumberFractions.begin()),
+ decltype(sigma.begin())>(
+ fNumberFractions.begin(), sigma.begin()),
+ nuc_comp::detail::WeightProviderIterator<decltype(fNumberFractions.begin()),
+ decltype(sigma.end())>(
+ fNumberFractions.end(), sigma.end()));
+
+ auto const iChannel = channelDist(randomStream);
+ return fComponents[iChannel];
+ }
+ };
+
+} // namespace corsika
diff --git a/corsika/media/ShowerAxis.hpp b/corsika/media/ShowerAxis.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..49855030520c40bc0505bbd1b7ffd8359a152822
--- /dev/null
+++ b/corsika/media/ShowerAxis.hpp
@@ -0,0 +1,80 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <cassert>
+#include <cstdlib>
+#include <fstream>
+#include <functional>
+#include <iterator>
+#include <memory>
+#include <stdexcept>
+#include <vector>
+
+#include <iostream>
+
+#include <boost/math/quadrature/gauss_kronrod.hpp>
+
+namespace corsika {
+
+ /**
+ * \class ShowerAxis
+ *
+ * The environment::ShowerAxis is created from a geometry::Point and
+ * a geometry::Vector and inside an Environment. It internally uses
+ * a table with steps=10000 (default) rows for interpolation.
+ *
+ * The shower axis can convert location in the shower into a
+ * projected grammage along the shower axis.
+ *
+ **/
+
+ class ShowerAxis {
+ public:
+ template <typename TEnvModel>
+ ShowerAxis(corsika::Point const& pStart,
+ corsika::Vector<units::si::length_d> length,
+ Environment<TEnvModel> const& env, int steps = 10'000);
+
+ units::si::LengthType steplength() const;
+
+ units::si::GrammageType maximumX() const;
+
+ units::si::GrammageType minimumX() const;
+
+ units::si::GrammageType projectedX(Point const& p) const;
+
+ units::si::GrammageType X(units::si::LengthType) const;
+
+ Vector<units::si::dimensionless_d> const& getDirection() const;
+
+ Point const& getStart() const;
+
+ private:
+ Point const pointStart_;
+ Vector<units::si::length_d> const length_;
+ units::si::LengthType const max_length_, steplength_;
+ Vector<units::si::dimensionless_d> const axis_normalized_;
+ std::vector<units::si::GrammageType> X_;
+
+ // for storing the lengths corresponding to equidistant X values
+ units::si::GrammageType const X_binning_ = std::invoke([]() {
+ using namespace units::si;
+ return 1_g / 1_cm / 1_cm;
+ });
+ std::vector<units::si::LengthType> d_;
+ };
+} // namespace corsika
+
+#include <corsika/detail/media/ShowerAxis.inl>
\ No newline at end of file
diff --git a/corsika/media/UniformMagneticField.hpp b/corsika/media/UniformMagneticField.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cdde91eabf653c23acdda9b433adbcf158150a44
--- /dev/null
+++ b/corsika/media/UniformMagneticField.hpp
@@ -0,0 +1,68 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+ /**
+ * A uniform (constant) magnetic field.
+ *
+ * This class returns the same magnetic field vector
+ * for all evaluated locations.
+ *
+ */
+ template <typename T>
+ class UniformMagneticField : public T {
+
+ // a type-alias for a magnetic field vector
+ using MagneticFieldVector = corsika::Vector<phys::units::magnetic_flux_density_d>;
+
+ MagneticFieldVector B_; ///< The constant magnetic field we use.
+
+ public:
+ /**
+ * Construct a UniformMagneticField.
+ *
+ * This is initialized with a fixed magnetic field
+ * and returns this magnetic field at all locations.
+ *
+ * @param field The fixed magnetic field to return.
+ */
+ template <typename... Args>
+ UniformMagneticField(MagneticFieldVector const& B, Args&&... args)
+ : T(std::forward<Args>(args)...)
+ , B_(B) {}
+
+ /**
+ * Evaluate the magnetic field at a given location.
+ *
+ * @param point The location to evaluate the field at.
+ * @returns The magnetic field vector.
+ */
+ MagneticFieldVector GetMagneticField(
+ corsika::geometry::Point const&) const final override {
+ return B_;
+ }
+
+ /**
+ * Set the magnetic field returned by this instance.
+ *
+ * @param point The location to evaluate the field at.
+ * @returns The magnetic field vector.
+ */
+ auto SetMagneticField(MagneticFieldVector const& Bfield) -> void { B_ = Bfield; }
+
+ }; // END: class MagneticField
+
+} // namespace corsika
+
+#include <corsika/detail/media/UniformMagneticField.inl>
\ No newline at end of file
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index 005c6235074eff3c84f04b93a4dafda6ca629518..989958a45df8ce7ce97ffe03b07e335a847d3a2d 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -312,7 +312,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
CHECK(builder.size() == 2);
auto const builtEnv = builder.assemble();
- auto const& univ = builtEnv.GetUniverse();
+ auto const& univ = builtEnv.getUniverse();
CHECK(builder.size() == 0);
CHECK(univ->GetChildNodes().size() == 1);
diff --git a/tests/media/testShowerAxis.cpp b/tests/media/testShowerAxis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f3df2264a017cb3187fbf63819818c220fb1d496
--- /dev/null
+++ b/tests/media/testShowerAxis.cpp
@@ -0,0 +1,80 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/media/DensityFunction.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/VolumeTreeNode.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <catch2/catch.hpp>
+
+//using namespace
+using namespace corsika;
+
+const auto density = 1_kg / (1_m * 1_m * 1_m);
+
+auto setupEnvironment(Code vTargetCode) {
+ // setup environment, geometry
+ auto env = std::make_unique<Environment<IMediumModel>>();
+ auto& universe = *(env->getUniverse());
+ const CoordinateSystem& cs = env->getCoordinateSystem();
+
+ auto theMedium =
+ Environment<IMediumModel>::CreateNode<Sphere>(
+ Point{cs, 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ density, NuclearComposition(std::vector<Code>{vTargetCode},
+ std::vector<float>{1.}));
+
+ auto const* nodePtr = theMedium.get();
+ universe.AddChild(std::move(theMedium));
+
+ return std::make_tuple(std::move(env), &cs, nodePtr);
+}
+
+TEST_CASE("Homogeneous Density") {
+ auto [env, csPtr, nodePtr] = setupEnvironment(Code::Nitrogen);
+ auto const& cs = *csPtr;
+ [[maybe_unused]] auto const& env_dummy = env;
+ [[maybe_unused]] auto const& node_dummy = nodePtr;
+
+ auto const observationHeight = 0_km;
+ auto const injectionHeight = 10_km;
+ auto const t = -observationHeight + injectionHeight;
+ Point const showerCore{cs, 0_m, 0_m, observationHeight};
+ Point const injectionPos = showerCore + Vector<dimensionless_d>{cs, {0, 0, 1}} * t;
+
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos),
+ *env, 20};
+
+ CHECK(showerAxis.steplength() == 500_m);
+
+ CHECK(showerAxis.maximumX() / (10_km * density) == Approx(1).epsilon(1e-8));
+
+ CHECK(showerAxis.minimumX() == 0_g / square(1_cm));
+
+ const Point p{cs, 10_km, 20_km, 8.3_km};
+ CHECK(showerAxis.projectedX(p) / (1.7_km * density) == Approx(1).epsilon(1e-8));
+
+ const units::si::LengthType d = 6.789_km;
+ CHECK(showerAxis.X(d) / (d * density) == Approx(1).epsilon(1e-8));
+
+ const Vector<dimensionless_d> dir{cs, {0, 0, -1}};
+ CHECK(showerAxis.getDirection().GetComponents(cs) == dir.GetComponents(cs));
+ CHECK(showerAxis.getStart().GetCoordinates() == injectionPos.GetCoordinates());
+}