From 911947804cda4a2977b617b61696e219b1531aca Mon Sep 17 00:00:00 2001 From: Dominik Baack <dominik.baack@tu-dortmund.de> Date: Thu, 26 Nov 2020 21:45:30 +0100 Subject: [PATCH] Several modifications to make the code more standard conform and document it --- corsika/detail/media/BaseExponential.inl | 4 +- corsika/detail/media/Environment.inl | 2 +- corsika/detail/media/FlatExponential.inl | 4 +- corsika/detail/media/HomogeneousMedium.inl | 5 +- corsika/detail/media/InhomogeneousMedium.inl | 8 +- .../LayeredSphericalAtmosphereBuilder.inl | 4 + .../media/LinearApproximationIntegrator.inl | 4 +- corsika/detail/media/MediumPropertyModel.inl | 2 +- corsika/detail/media/NuclearComposition.inl | 40 +++--- corsika/detail/media/ShowerAxis.inl | 10 +- corsika/media/Environment.hpp | 24 +++- corsika/media/IMagneticFieldModel.hpp | 4 +- corsika/media/MediumProperties.hpp | 46 +++--- corsika/media/MediumPropertyModel.hpp | 2 +- corsika/media/NuclearComposition.hpp | 66 ++++++--- corsika/media/NuclearComposition_old.hpp | 136 ------------------ corsika/media/ShowerAxis.hpp | 7 +- corsika/media/UniformMagneticField.hpp | 9 +- src/media/CMakeLists.txt | 2 +- src/media/readProperties.py | 38 ++--- tests/media/testMedium.cpp | 33 ++--- tests/media/testRefractiveIndex.cpp | 2 +- tests/media/testShowerAxis.cpp | 2 +- 23 files changed, 193 insertions(+), 261 deletions(-) delete mode 100644 corsika/media/NuclearComposition_old.hpp diff --git a/corsika/detail/media/BaseExponential.inl b/corsika/detail/media/BaseExponential.inl index 7cfac790e..96ce39943 100644 --- a/corsika/detail/media/BaseExponential.inl +++ b/corsika/detail/media/BaseExponential.inl @@ -24,7 +24,7 @@ namespace corsika { } template <typename TDerived> - units::si::GrammageType BaseExponential<TDerived>::integratedGrammage( + GrammageType BaseExponential<TDerived>::integratedGrammage( Trajectory<Line> const& line, units::si::LengthType vL, Vector<units::si::dimensionless_d> const& axis) const { if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); } @@ -40,7 +40,7 @@ namespace corsika { } template <typename TDerived> - units::si::LengthType BaseExponential<TDerived>::arclengthFromGrammage( + LengthType BaseExponential<TDerived>::getArclengthFromGrammage( Trajectory<Line> const& line, units::si::GrammageType grammage, Vector<units::si::dimensionless_d> const& axis) const { auto const uDotA = line.NormalizedDirection().dot(axis).magnitude(); diff --git a/corsika/detail/media/Environment.inl b/corsika/detail/media/Environment.inl index 434d9c77c..b1fc4bc62 100644 --- a/corsika/detail/media/Environment.inl +++ b/corsika/detail/media/Environment.inl @@ -38,7 +38,7 @@ namespace corsika { // factory method for creation of VolumeTreeNodes template <typename IEnvironmentModel> template <class TVolumeType, typename... TVolumeArgs> - std::unique_ptr< VolumeTreeNode<IEnvironmentModel> > Environment<IEnvironmentModel>::CreateNode(TVolumeArgs&&... args) { + 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\""); diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl index 34b94e0a2..e98d2822b 100644 --- a/corsika/detail/media/FlatExponential.inl +++ b/corsika/detail/media/FlatExponential.inl @@ -43,13 +43,13 @@ namespace corsika { } template <typename T> - units::si::GrammageType FlatExponential<T>::integratedGrammage( + GrammageType FlatExponential<T>::integratedGrammage( Trajectory<Line> const& line, units::si::LengthType to) const { return BaseExponential<FlatExponential<T>>::integratedGrammage(line, to, axis_); } template <typename T> - units::si::LengthType FlatExponential<T>::arclengthFromGrammage( + LengthType FlatExponential<T>::getArclengthFromGrammage( Trajectory<Line> const& line, units::si::GrammageType grammage) const { return BaseExponential<FlatExponential<T>>::arclengthFromGrammage(line, grammage, axis_); diff --git a/corsika/detail/media/HomogeneousMedium.inl b/corsika/detail/media/HomogeneousMedium.inl index adbe03c4f..090125388 100644 --- a/corsika/detail/media/HomogeneousMedium.inl +++ b/corsika/detail/media/HomogeneousMedium.inl @@ -32,14 +32,13 @@ namespace corsika { } template <typename T> - units::si::GrammageType HomogeneousMedium<T>::integratedGrammage( + GrammageType HomogeneousMedium<T>::integratedGrammage( Trajectory<Line> const&, units::si::LengthType to) const { - using namespace units::si; return to * density_; } template <typename T> - units::si::LengthType HomogeneousMedium<T>::arclengthFromGrammage( + LengthType HomogeneousMedium<T>::getArclengthFromGrammage( Trajectory<Line> const&, units::si::GrammageType grammage) const { return grammage / density_; } diff --git a/corsika/detail/media/InhomogeneousMedium.inl b/corsika/detail/media/InhomogeneousMedium.inl index ff64e9cd1..08e5f7af5 100644 --- a/corsika/detail/media/InhomogeneousMedium.inl +++ b/corsika/detail/media/InhomogeneousMedium.inl @@ -19,7 +19,7 @@ namespace corsika { template <typename T, typename TDensityFunction> template <typename... TArgs> InhomogeneousMedium<T, TDensityFunction>::InhomogeneousMedium( - NuclearComposition nuclComp, TArgs&&... rhoTArgs) + NuclearComposition const& nuclComp, TArgs&&... rhoTArgs) : nuclComp_(nuclComp) , densityFunction_(rhoTArgs...){} @@ -36,15 +36,15 @@ namespace corsika { } template <typename T, typename TDensityFunction> - units::si::GrammageType InhomogeneousMedium<T, TDensityFunction>::integratedGrammage( + GrammageType InhomogeneousMedium<T, TDensityFunction>::integratedGrammage( Trajectory<Line> const& line, units::si::LengthType to) const { return densityFunction_.integrateGrammage(line, to); } template <typename T, typename TDensityFunction> - units::si::LengthType InhomogeneousMedium<T, TDensityFunction>::arclengthFromGrammage( + LengthType InhomogeneousMedium<T, TDensityFunction>::getArclengthFromGrammage( Trajectory<Line> const& line, units::si::GrammageType grammage) const { - return densityFunction_.arclengthFromGrammage(line, grammage); + return densityFunction_.getArclengthFromGrammage(line, grammage); } } // namespace corsika diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl index 340a25af4..fa6cd992d 100644 --- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl +++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl @@ -8,6 +8,10 @@ * the license. */ +#pragma once + +#include <corsika/framework/logging/Logging.hpp> + #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp> #include <corsika/media/FlatExponential.hpp> #include <corsika/media/HomogeneousMedium.hpp> diff --git a/corsika/detail/media/LinearApproximationIntegrator.inl b/corsika/detail/media/LinearApproximationIntegrator.inl index 98d52db1c..2826251b2 100644 --- a/corsika/detail/media/LinearApproximationIntegrator.inl +++ b/corsika/detail/media/LinearApproximationIntegrator.inl @@ -29,7 +29,7 @@ namespace corsika { } template <typename TDerived> - auto LinearApproximationIntegrator<TDerived>::arclengthFromGrammage( + auto LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage( Trajectory<Line> const& line, units::si::GrammageType grammage) const { auto const c0 = getImplementation().rho_(line.GetPosition(0)); auto const c1 = getImplementation().rho_.FirstDerivative(line.GetPosition(0), @@ -39,7 +39,7 @@ namespace corsika { } template <typename TDerived> - auto LinearApproximationIntegrator<TDerived>::maximumLength( + auto LinearApproximationIntegrator<TDerived>::getMaximumLength( Trajectory<Line> const& line, [[maybe_unused]] double relError) const { using namespace units::si; [[maybe_unused]] auto const c1 = getImplementation().rho_.SecondDerivative( diff --git a/corsika/detail/media/MediumPropertyModel.inl b/corsika/detail/media/MediumPropertyModel.inl index 3e584252b..90cdd9754 100644 --- a/corsika/detail/media/MediumPropertyModel.inl +++ b/corsika/detail/media/MediumPropertyModel.inl @@ -19,7 +19,7 @@ namespace corsika { , medium_(medium) {} template <typename T> - Medium MediumPropertyModel<T>::getMedium(Point const&) const final override { + Medium MediumPropertyModel<T>::getMedium(Point const&) const { return medium_; } diff --git a/corsika/detail/media/NuclearComposition.inl b/corsika/detail/media/NuclearComposition.inl index ba56e6ea3..f1107f7b1 100644 --- a/corsika/detail/media/NuclearComposition.inl +++ b/corsika/detail/media/NuclearComposition.inl @@ -1,5 +1,5 @@ /* - * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * (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 @@ -58,17 +58,18 @@ namespace corsika { return !(*this == other); } - NuclearComposition::NuclearComposition(std::vector<corsika::Code> pComponents, - std::vector<float> pFractions) + NuclearComposition::NuclearComposition(std::vector<corsika::Code> const& pComponents, + std::vector<float> const& 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; + if (is_nucleus(compID)) { + return get_nucleus_A(compID) * fraction; } else { - return GetMass(compID) / units::si::ConvertSIToHEP(constants::u) * fraction; + return get_mass(compID) / units::si::ConvertSIToHEP(constants::u) * + fraction; } })) { assert(pComponents.size() == pFractions.size()); @@ -78,11 +79,11 @@ namespace corsika { if (!(0.999f < sumFractions && sumFractions < 1.001f)) { throw std::runtime_error("element fractions do not add up to 1"); } - updateHash(); + this->updateHash(); } template <typename TFunction> - auto NuclearComposition::WeightedSum(TFunction func) const { + auto NuclearComposition::weightedSum(TFunction const& func) const { using ResultQuantity = decltype(func(*components_.cbegin())); auto const prod = [&](auto const compID, auto const fraction) { @@ -104,14 +105,19 @@ namespace corsika { 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_; } + std::vector<float> const& NuclearComposition::getFractions() const { + return numberFractions_; + } + + std::vector<corsika::Code> 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 { + corsika::Code NuclearComposition::sampleTarget( + std::vector<units::si::CrossSectionType> const& sigma, TRNG& randomStream) const { using namespace units::si; assert(sigma.size() == numberFractions_.size()); @@ -133,10 +139,10 @@ namespace corsika { 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()) + for (float ifrac : this->getFractions()) hashes.push_back(std::hash<float>{}(ifrac)); + for (corsika::Code icode : this->getComponents()) hashes.push_back(std::hash<int>{}(static_cast<int>(icode))); - std::size_t h = std::hash<double>{}(GetAverageMassNumber()); + std::size_t h = std::hash<double>{}(this->getAverageMassNumber()); for (std::size_t ih : hashes) h = h ^ (ih << 1); hash_ = h; } diff --git a/corsika/detail/media/ShowerAxis.inl b/corsika/detail/media/ShowerAxis.inl index 2f5691794..38e58ef1a 100644 --- a/corsika/detail/media/ShowerAxis.inl +++ b/corsika/detail/media/ShowerAxis.inl @@ -15,7 +15,7 @@ namespace corsika { template <typename TEnvModel> ShowerAxis::ShowerAxis(Point const& pStart, - corsika::Vector<units::si::length_d> length, + corsika::Vector<units::si::length_d> const& length, Environment<TEnvModel> const& env, int steps) : pointStart_(pStart) @@ -58,7 +58,7 @@ namespace corsika { assert(std::is_sorted(d_.cbegin(), d_.cend())); } - GrammageType ShowerAxis::X(LengthType l) const { + GrammageType ShowerAxis::getX(LengthType l) const { auto const fractionalBin = l / steplength_; int const lower = fractionalBin; // indices of nearest X support points auto const lambda = fractionalBin - lower; @@ -87,11 +87,11 @@ namespace corsika { return X_[upper] * lambda + X_[lower] * (1 - lambda); } - LengthType ShowerAxis::steplength() const { return steplength_; } + LengthType ShowerAxis::getSteplength() const { return steplength_; } - GrammageType ShowerAxis::maximumX() const { return *X_.rbegin(); } + GrammageType ShowerAxis::getMaximumX() const { return *X_.rbegin(); } - GrammageType ShowerAxis::minimumX() const { return GrammageType::zero(); } + GrammageType ShowerAxis::getMinimumX() const { return GrammageType::zero(); } GrammageType ShowerAxis::projectedX(Point const& p) const { auto const projectedLength = (p - pointStart_).dot(axis_normalized_); diff --git a/corsika/media/Environment.hpp b/corsika/media/Environment.hpp index fd2ccc74c..f9cf69d82 100644 --- a/corsika/media/Environment.hpp +++ b/corsika/media/Environment.hpp @@ -20,6 +20,9 @@ namespace corsika { + /** Base Evnironment class + * Describes the Environment in which the shower is propagated + **/ template <typename IEnvironmentModel> class Environment { public: @@ -27,14 +30,31 @@ namespace corsika { Environment(); + /** Getters for the universe stored in the Environment + * + * @retval Retuns reference to a Universe object with infinite size + **/ + ///@{ + //* Get non const universe */ typename BaseNodeType::VTNUPtr& getUniverse(); + //* Get const universe */ typename BaseNodeType::VTNUPtr const& getUniverse() const; + ///@} + /** Getter for the CoordinateSystem used in the Environment + * + * @retval Retuns a const reference to the CoordinateSystem used + **/ CoordinateSystem const& getCoordinateSystem() const; - // factory method for creation of VolumeTreeNodes + /** Factory method for creation of VolumeTreeNodes + * @tparam TVolumeType Type of volume to be created + * @tparam TVolumeArgs Types to forward to the constructor + * @param args Parameter forwarded to the constructor of TVolumeType + * @retval Retuns unique pointer to a VolumeTreeNode with the same EnvitonmentModel as this class + **/ template <class TVolumeType, typename... TVolumeArgs> - static std::unique_ptr<BaseNodeType> CreateNode(TVolumeArgs&&... args); + static std::unique_ptr<BaseNodeType> createNode(TVolumeArgs&&... args); private: CoordinateSystem const& coordinateSystem_; diff --git a/corsika/media/IMagneticFieldModel.hpp b/corsika/media/IMagneticFieldModel.hpp index c7c62fab8..f110feb97 100644 --- a/corsika/media/IMagneticFieldModel.hpp +++ b/corsika/media/IMagneticFieldModel.hpp @@ -34,7 +34,7 @@ namespace corsika { * @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 + virtual auto getMagneticField(corsika::geometry::Point const&) const -> MagneticFieldVector = 0; /** @@ -44,4 +44,4 @@ namespace corsika { }; // END: class MagneticField -} // namespace corsika::environment +} // namespace corsika diff --git a/corsika/media/MediumProperties.hpp b/corsika/media/MediumProperties.hpp index f01c128e0..5fb0efeed 100644 --- a/corsika/media/MediumProperties.hpp +++ b/corsika/media/MediumProperties.hpp @@ -33,7 +33,7 @@ namespace corsika { enum class Medium : int16_t; using MediumIntType = std::underlying_type<Medium>::type; - /** + /** \todo documentation needs update ... * \struct MediumData * * Simple object to group together a number of properties @@ -59,34 +59,34 @@ namespace corsika { 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_; } + std::string getName() const { return name_; } + std::string getPrettyName() const { return pretty_name_; } + double getWeight() const { return weight_; } + const int& weight_significant_figure() const { return weight_significant_figure_; } + const int& weight_error_last_digit() const { return weight_error_last_digit_; } + const double& Z_over_A() const { return Z_over_A_; } + double getSternheimerDensity() const { return sternheimer_density_; } + double getCorrectedDensity() const { return corrected_density_; } + State getState() const { return state_; } + MediumType getType() const { return type_; } + std::string getSymbol() const { return symbol_; } + double getIeff() const { return Ieff_; } + double getCbar() const { return Cbar_; } + double getX0() const { return x0_; } + double getX1() const { return x1_; } + double getAA() const { return aa_; } + double getSK() const { return sk_; } + double getDlt0() const { return dlt0_; } }; } // namespace corsika -#include <corsika/detail/media/GeneratedMediaProperties.inc> +#include <corsika/media/GeneratedMediaProperties.inc> -namespace corsika::environment { +namespace corsika { constexpr MediumData const& mediumData(Medium const m) { - return detail::medium_data[static_cast<MediumIntType>(m)]; + return corsika::detail::medium_data[static_cast<MediumIntType>(m)]; } -} // namespace corsika::environment +} // namespace corsika::medium diff --git a/corsika/media/MediumPropertyModel.hpp b/corsika/media/MediumPropertyModel.hpp index bf47b9514..7833599c6 100644 --- a/corsika/media/MediumPropertyModel.hpp +++ b/corsika/media/MediumPropertyModel.hpp @@ -44,7 +44,7 @@ namespace corsika { * @param medium The medium to store. * @returns The medium type as enum environment::Medium */ - void setMedium(Medium const medium) override; + void setMedium(Medium const medium); }; // END: class MediumPropertyModel diff --git a/corsika/media/NuclearComposition.hpp b/corsika/media/NuclearComposition.hpp index feb30c3da..513ebef7e 100644 --- a/corsika/media/NuclearComposition.hpp +++ b/corsika/media/NuclearComposition.hpp @@ -1,5 +1,5 @@ /* - * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * (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 @@ -19,15 +19,21 @@ #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_; - - std::size_t hash_; + /** Describes the composition of matter + * Allowes and handles the creation of custom matter compositions + **/ + class NuclearComposition { + private: + /// TODO: replace with + /// https://www.boost.org/doc/libs/1_74_0/libs/iterator/doc/zip_iterator.html or + /// ranges zip + /** Double Iterator + * Iterator that allowes the iteration of two individual lists at the same time. The + *user needs to take care that booth lists have the same length. + * @tparam AConstIterator Iterator Type of the first list + * @tparam BConstIterator Iterator Type of the second list + **/ template <class AConstIterator, class BConstIterator> class WeightProviderIterator { AConstIterator aIter_; @@ -52,22 +58,43 @@ namespace corsika { }; public: + /** Constructor + * The constructore takes a list of elements and a list which describe the relative + * amount. Booth lists need to have the same length and the sum all of fractions + * should be 1. Otherwise an exception is thrown + * @param pComponents List of particle types + * @param pFractions List of fractions how much each particle contributes. The sum + *needs to add up to 1 + **/ NuclearComposition(std::vector<corsika::Code> pComponents, std::vector<float> pFractions); + /** Sum all all relative composition weighted by func(element) + * This function sums all relative compositions given during this classes + *construction. Each entry is weighted by the user defined function func given to this + *function. + * @tparam TFunction Type of functions for the weights. The type should be + *corsika::Code -> float + * @param func Functions for reweighting specific elements + * @retval returns the weighted sum with the type defined by the return type of func + **/ template <typename TFunction> - auto WeightedSum(TFunction func) const ; + auto weightedSum(TFunction func) const; + /** Number of elements in the composition array + * @retval returns the number of elements in the composition array + **/ auto size() const; - auto const& GetFractions() const; - auto const& GetComponents() const ; - auto const GetAverageMassNumber() const; + /// Returns a const reference to the fraction + std::vector<float> const& getFractions() const; + /// Returns a const reference to the fraction + std::vector<corsika::Code> const& getComponents() const; + auto const getAverageMassNumber() const; template <class TRNG> - Code SampleTarget( - std::vector<units::si::CrossSectionType> const& sigma, - TRNG& randomStream) const; + 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! @@ -75,6 +102,13 @@ namespace corsika { private: void updateHash(); + + std::vector<float> const numberFractions_; //!< relative fractions of number density + std::vector<corsika::Code> const components_; //!< particle codes of consitutents + + double const avgMassNumber_; + + std::size_t hash_; }; } // namespace corsika diff --git a/corsika/media/NuclearComposition_old.hpp b/corsika/media/NuclearComposition_old.hpp deleted file mode 100644 index 82577f245..000000000 --- a/corsika/media/NuclearComposition_old.hpp +++ /dev/null @@ -1,136 +0,0 @@ -/* - * (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 index 498550305..b58bc9128 100644 --- a/corsika/media/ShowerAxis.hpp +++ b/corsika/media/ShowerAxis.hpp @@ -31,8 +31,8 @@ namespace corsika { /** * \class ShowerAxis * - * The environment::ShowerAxis is created from a geometry::Point and - * a geometry::Vector and inside an Environment. It internally uses + * The environment::ShowerAxis is created from a Point and + * a 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 @@ -40,6 +40,7 @@ namespace corsika { * **/ + ///\todo documentation needs update ... class ShowerAxis { public: template <typename TEnvModel> @@ -55,7 +56,7 @@ namespace corsika { units::si::GrammageType projectedX(Point const& p) const; - units::si::GrammageType X(units::si::LengthType) const; + GrammageType getX(LengthType) const; Vector<units::si::dimensionless_d> const& getDirection() const; diff --git a/corsika/media/UniformMagneticField.hpp b/corsika/media/UniformMagneticField.hpp index cdde91eab..e85ae5d36 100644 --- a/corsika/media/UniformMagneticField.hpp +++ b/corsika/media/UniformMagneticField.hpp @@ -26,8 +26,6 @@ namespace corsika { // 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. @@ -48,7 +46,7 @@ namespace corsika { * @param point The location to evaluate the field at. * @returns The magnetic field vector. */ - MagneticFieldVector GetMagneticField( + MagneticFieldVector getMagneticField( corsika::geometry::Point const&) const final override { return B_; } @@ -59,7 +57,10 @@ namespace corsika { * @param point The location to evaluate the field at. * @returns The magnetic field vector. */ - auto SetMagneticField(MagneticFieldVector const& Bfield) -> void { B_ = Bfield; } + auto setMagneticField(MagneticFieldVector const& Bfield) -> void { B_ = Bfield; } + + private: + MagneticFieldVector B_; ///< The constant magnetic field we use. }; // END: class MagneticField diff --git a/src/media/CMakeLists.txt b/src/media/CMakeLists.txt index 1fb83c74d..94b578560 100644 --- a/src/media/CMakeLists.txt +++ b/src/media/CMakeLists.txt @@ -3,7 +3,7 @@ set (output_dir ${PROJECT_BINARY_DIR}/corsika/media) file (MAKE_DIRECTORY ${output_dir}) -add_custom_command ( +add_custom_command ( OUTPUT ${output_dir}/GeneratedMediaProperties.inc COMMAND ${input_dir}/readProperties.py ${input_dir}/properties8.dat DEPENDS ${input_dir}/readProperties.py diff --git a/src/media/readProperties.py b/src/media/readProperties.py index 552c70d70..44eb3f839 100755 --- a/src/media/readProperties.py +++ b/src/media/readProperties.py @@ -288,25 +288,25 @@ def gen_classes(media_db): public: static constexpr Medium medium() {{ return Medium::{cname}; }} - static std::string const name() {{ return data_.name(); }} - static std::string const pretty_name() {{ return data_.pretty_name(); }} - static double weight() {{ return data_.weight (); }} + static std::string const getName() {{ return data_.getName(); }} + static std::string const getPrettyName() {{ return data_.getPrettyName(); }} + static double getWeight() {{ return data_.getWeight (); }} static int weight_significant_figure() {{ return data_.weight_significant_figure (); }} static int weight_error_last_digit() {{ return data_.weight_error_last_digit(); }} static double Z_over_A() {{ return data_.Z_over_A(); }} - static double sternheimer_density() {{ return data_.sternheimer_density(); }} - static double corrected_density() {{ return data_.corrected_density(); }} - static State state() {{ return data_.state(); }} - static MediumType type() {{ return data_.type(); }} - static std::string const symbol() {{ return data_.symbol(); }} - - static double Ieff() {{ return data_.Ieff(); }} - static double Cbar() {{ return data_.Cbar(); }} - static double x0() {{ return data_.x0(); }} - static double x1() {{ return data_.x1(); }} - static double aa() {{ return data_.aa(); }} - static double sk() {{ return data_.sk(); }} - static double dlt0() {{ return data_.dlt0(); }} + static double getSternheimerDensity() {{ return data_.getSternheimerDensity(); }} + static double getCorrectedDensity() {{ return data_.getCorrectedDensity(); }} + static State getState() {{ return data_.getState(); }} + static MediumType getType() {{ return data_.getType(); }} + static std::string const getSymbol() {{ return data_.getSymbol(); }} + + static double getIeff() {{ return data_.getIeff(); }} + static double getCbar() {{ return data_.getCbar(); }} + static double getX0() {{ return data_.getX0(); }} + static double getX1() {{ return data_.getX1(); }} + static double getAA() {{ return data_.getAA(); }} + static double getSK() {{ return data_.getSK(); }} + static double getDlt0() {{ return data_.getDlt0(); }} //static constexpr Constituents constituents() {{ return {constituents}; }} //static constexpr Properties properties() {{ return {properties}; }} @@ -416,7 +416,9 @@ def inc_start(): // generated by readProperties.py // MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN. // since this is automatic code, we don't attempt to generate automatic unit testing, too: LCOV_EXCL_START -namespace corsika::environment { + +#pragma once +namespace corsika { """ return string @@ -442,7 +444,7 @@ def detail_end(): # def inc_end(): string = """ -\n} // end namespace corsika::environment +\n} // end namespace corsika // since this was automatic code, we didn't attempt to generate automatic unit testing, too: LCOV_EXCL_STOP """ return string diff --git a/tests/media/testMedium.cpp b/tests/media/testMedium.cpp index d963c5d1e..d9dcf0c86 100644 --- a/tests/media/testMedium.cpp +++ b/tests/media/testMedium.cpp @@ -16,6 +16,7 @@ #include <corsika/media/InhomogeneousMedium.hpp> #include <corsika/media/NuclearComposition.hpp> #include <corsika/media/MediumPropertyModel.hpp> +#include <corsika/media/MediumProperties.hpp> #include <catch2/catch.hpp> @@ -25,7 +26,7 @@ using namespace corsika::units::si; TEST_CASE("MediumPropertyModel w/ Homogeneous") { CoordinateSystem const& gCS = - RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); + RootCoordinateSystem::getInstance().GetRootCoordinateSystem(); Point const gOrigin(gCS, {0_m, 0_m, 0_m}); @@ -41,35 +42,35 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") { std::vector<float>{1.f}); // the refrative index that we use - const Medium type = Medium::AirDry1Atm; + const Medium type = corsika::Medium::AirDry1Atm; // create the atmospheric model AtmModel medium(type, density, protonComposition); // and require that it is constant - CHECK(type == medium.medium(Point(gCS, -10_m, 4_m, 35_km))); - CHECK(type == medium.medium(Point(gCS, +210_m, 0_m, 7_km))); - CHECK(type == medium.medium(Point(gCS, 0_m, 0_m, 0_km))); - CHECK(type == medium.medium(Point(gCS, 100_km, 400_km, 350_km))); + CHECK(type == medium.getMedium(Point(gCS, -10_m, 4_m, 35_km))); + CHECK(type == medium.getMedium(Point(gCS, +210_m, 0_m, 7_km))); + CHECK(type == medium.getMedium(Point(gCS, 0_m, 0_m, 0_km))); + CHECK(type == medium.getMedium(Point(gCS, 100_km, 400_km, 350_km))); // a new refractive index - const Medium type2 = Medium::StandardRock; + const Medium type2 = corsika::Medium::StandardRock; // update the refractive index of this atmospheric model - medium.set_medium(type2); + medium.setMedium(type2); // check that the returned refractive index is correct - CHECK(type2 == medium.medium(Point(gCS, -10_m, 4_m, 35_km))); - CHECK(type2 == medium.medium(Point(gCS, +210_m, 0_m, 7_km))); - CHECK(type2 == medium.medium(Point(gCS, 0_m, 0_m, 0_km))); - CHECK(type2 == medium.medium(Point(gCS, 100_km, 400_km, 350_km))); + CHECK(type2 == medium.getMedium(Point(gCS, -10_m, 4_m, 35_km))); + CHECK(type2 == medium.getMedium(Point(gCS, +210_m, 0_m, 7_km))); + CHECK(type2 == medium.getMedium(Point(gCS, 0_m, 0_m, 0_km))); + CHECK(type2 == medium.getMedium(Point(gCS, 100_km, 400_km, 350_km))); // define our axis vector Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1)); // check the density and nuclear composition - CHECK(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m))); - medium.GetNuclearComposition(); + CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m))); + medium.getNuclearComposition(); // create a line of length 1 m Line const line(gOrigin, Vector<SpeedType::dimension_type>( @@ -82,6 +83,6 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") { Trajectory<Line> const trajectory(line, tEnd); // and check the integrated grammage - CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1)); - CHECK((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1)); + CHECK((medium.integratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1)); + CHECK((medium.arclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1)); } diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp index 871ec4148..79215e7ac 100644 --- a/tests/media/testRefractiveIndex.cpp +++ b/tests/media/testRefractiveIndex.cpp @@ -25,7 +25,7 @@ using namespace corsika::units::si; TEST_CASE("UniformRefractiveIndex w/ Homogeneous") { CoordinateSystem const& gCS = - RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); + RootCoordinateSystem::getInstance().GetRootCoordinateSystem(); Point const gOrigin(gCS, {0_m, 0_m, 0_m}); diff --git a/tests/media/testShowerAxis.cpp b/tests/media/testShowerAxis.cpp index f3df2264a..fcf2b920d 100644 --- a/tests/media/testShowerAxis.cpp +++ b/tests/media/testShowerAxis.cpp @@ -32,7 +32,7 @@ auto setupEnvironment(Code vTargetCode) { const CoordinateSystem& cs = env->getCoordinateSystem(); auto theMedium = - Environment<IMediumModel>::CreateNode<Sphere>( + Environment<IMediumModel>::createNode<Sphere>( Point{cs, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); -- GitLab