IAP GITLAB

Commit 91194780 authored by Dominik Baack's avatar Dominik Baack Committed by Ralf Ulrich

Several modifications to make the code more standard conform and document it

parent d19881b5
......@@ -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();
......
......@@ -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\"");
......
......@@ -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_);
......
......@@ -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_;
}
......
......@@ -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
......@@ -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>
......
......@@ -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(
......
......@@ -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_;
}
......
/*
* (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;
}
......
......@@ -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_);
......
......@@ -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_;
......
......@@ -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
......@@ -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
......@@ -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
......
/*
* (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
......
/*
* (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
......@@ -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;