IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 2520 deletions
Framework/Random/RNGManager.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/random/RNGManager.h>
void corsika::random::RNGManager::RegisterRandomStream(std::string const& pStreamName) {
corsika::random::RNG rng;
if (auto const& it = seeds.find(pStreamName); it != seeds.end()) {
rng.seed(it->second);
}
rngs[pStreamName] = std::move(rng);
}
corsika::random::RNG& corsika::random::RNGManager::GetRandomStream(
std::string const& pStreamName) {
return rngs.at(pStreamName);
}
std::stringstream corsika::random::RNGManager::dumpState() const {
std::stringstream buffer;
for (auto const& [streamName, rng] : rngs) {
buffer << '"' << streamName << "\" = \"" << rng << '"' << std::endl;
}
return buffer;
}
/*
void corsika::random::RNGManager::SetSeedSeq(std::string const& pStreamName,
std::seed_seq const& pSeedSeq) {
seeds[pStreamName] = pSeedSeq;
}
*/
Framework/Random/RNGManager.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_RNGManager_h_
#define _include_RNGManager_h_
#include <corsika/utl/Singleton.h>
#include <map>
#include <random>
#include <sstream>
#include <string>
/*!
* With this class modules can register streams of random numbers.
*/
namespace corsika::random {
using RNG = std::mt19937; //!< the actual RNG type that will be used
class RNGManager : public corsika::utl::Singleton<RNGManager> {
friend class corsika::utl::Singleton<RNGManager>;
std::map<std::string, RNG> rngs;
std::map<std::string, std::seed_seq> seeds;
protected:
RNGManager() {}
public:
/*!
* This function is to be called by a module requiring a random-number
* stream during its initialization.
*
* \throws sth. when stream \a pModuleName is already registered
*/
void RegisterRandomStream(std::string const& pStreamName);
/*!
* returns the pre-stored stream of given name \a pStreamName if
* available
*/
RNG& GetRandomStream(std::string const& pStreamName);
/*!
* dumps the names and states of all registered random-number streams
* into a std::stringstream.
*/
std::stringstream dumpState() const;
/**
* set seed_seq of \a pStreamName to \a pSeedSeq
*/
// void SetSeedSeq(std::string const& pStreamName, std::seed_seq& const pSeedSeq);
};
} // namespace corsika::random
#endif
Framework/Random/UniformRealDistribution.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_UniformRealDistribution_h
#define _include_UniformRealDistribution_h
#include <corsika/units/PhysicalUnits.h>
#include <random>
namespace corsika::random {
template <class TQuantity>
class UniformRealDistribution {
using RealType = typename TQuantity::value_type;
std::uniform_real_distribution<RealType> dist{RealType(0.), RealType(1.)};
TQuantity const a, b;
public:
UniformRealDistribution(TQuantity b)
: a{TQuantity(phys::units::detail::magnitude_tag, 0)}
, b(b) {}
UniformRealDistribution(TQuantity a, TQuantity b)
: a(a)
, b(b) {}
template <class Generator>
TQuantity operator()(Generator& g) {
return a + dist(g) * (b - a);
}
};
} // namespace corsika::random
#endif
Framework/Random/testRandom.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
#include <corsika/random/ExponentialDistribution.h>
#include <corsika/random/RNGManager.h>
#include <corsika/random/UniformRealDistribution.h>
#include <corsika/units/PhysicalUnits.h>
#include <iostream>
#include <limits>
#include <random>
#include <type_traits>
using namespace corsika::random;
SCENARIO("random-number streams can be registered and retrieved") {
GIVEN("a RNGManager") {
RNGManager& rngManager = RNGManager::GetInstance();
WHEN("a sequence is registered by name") {
rngManager.RegisterRandomStream("stream_A");
THEN("the sequence can be retrieved") {
REQUIRE_NOTHROW(rngManager.GetRandomStream("stream_A"));
}
THEN("an unknown sequence cannot be retrieved") {
REQUIRE_THROWS(rngManager.GetRandomStream("stream_UNKNOWN"));
}
// seeding not covered yet
}
}
}
TEST_CASE("UniformRealDistribution") {
using namespace corsika::units::si;
std::mt19937 rng;
corsika::random::UniformRealDistribution<LengthType> dist(1_m, 2_m);
SECTION("range") {
corsika::random::UniformRealDistribution<LengthType> dist(1_m, 2_m);
LengthType min =
+1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
LengthType max =
-1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
for (int i{0}; i < 1'000'000; ++i) {
LengthType x = dist(rng);
min = std::min(min, x);
max = std::max(max, x);
}
CHECK(min / 1_m == Approx(1.));
CHECK(max / 2_m == Approx(1.));
}
SECTION("range") {
corsika::random::UniformRealDistribution<LengthType> dist(18_cm);
LengthType min =
+1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
LengthType max =
-1_m * std::numeric_limits<typename LengthType::value_type>::infinity();
for (int i{0}; i < 1'000'000; ++i) {
LengthType x = dist(rng);
min = std::min(min, x);
max = std::max(max, x);
}
CHECK(min / 1_m == Approx(0.).margin(1e-3));
CHECK(max / 18_cm == Approx(1.));
}
}
TEST_CASE("ExponentialDistribution") {
using namespace corsika::units::si;
std::mt19937 rng;
auto const beta = 15_m;
corsika::random::ExponentialDistribution dist(beta);
SECTION("mean") {
std::remove_const<decltype(beta)>::type mean = beta * 0;
int constexpr N = 1'000'000;
for (int i{0}; i < N; ++i) {
decltype(beta) x = dist(rng);
mean += x / N;
}
CHECK(mean / beta == Approx(1).margin(1e-2));
}
}
Framework/StackInterface/CMakeLists.txt deleted 100644 → 0
View file @ 6f2ab230
set (
CORSIKAstackinterface_HEADERS
ParticleBase.h
StackIteratorInterface.h
Stack.h
SecondaryView.h
CombinedStack.h
)
set (
CORSIKAstackinterface_NAMESPACE
corsika/stack
)
add_library (
CORSIKAstackinterface
INTERFACE
)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAstackinterface ${CORSIKAstackinterface_NAMESPACE} ${CORSIKAstackinterface_HEADERS})
target_include_directories (
CORSIKAstackinterface
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
FILES
${CORSIKAstackinterface_HEADERS}
DESTINATION
include/${CORSIKAstackinterface_NAMESPACE}
)
#code testing
add_executable (testStackInterface testStackInterface.cc)
target_link_libraries (testStackInterface CORSIKAstackinterface CORSIKAthirdparty) # for catch2
CORSIKA_ADD_TEST(testStackInterface)
add_executable (testSecondaryView testSecondaryView.cc)
target_link_libraries (testSecondaryView CORSIKAstackinterface CORSIKAthirdparty) # for catch2
CORSIKA_ADD_TEST(testSecondaryView)
add_executable (testCombinedStack testCombinedStack.cc)
target_link_libraries (testCombinedStack CORSIKAstackinterface CORSIKAthirdparty) # for catch2
CORSIKA_ADD_TEST(testCombinedStack)
Framework/StackInterface/CombinedStack.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_stack_combinedstack_h_
#define _include_stack_combinedstack_h_
#include <corsika/particles/ParticleProperties.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::stack {
/**
* @class CombinedParticleInterface
*
* You may combine two StackData object, see class CombinedStackImpl
* below, into one Stack, using a combined StackIterator (aka
* CombinedParticleInterface) interface class.
*
* This allows to add specific information to a given Stack, could
* be special information on a subset of entries
* (e.g. NuclearStackExtension) or also (multi) thinning weights for
* all particles.
*
* Many Stacks can be combined into more complex object.
*
* The two sub-stacks must both provide their independent
* ParticleInterface classes.
*
*/
template <template <typename> typename ParticleInterfaceA,
template <typename> typename ParticleInterfaceB, typename StackIterator>
class CombinedParticleInterface
: public ParticleInterfaceB<ParticleInterfaceA<StackIterator>> {
// template<template <typename> typename _PI>
// template <typename StackDataType, template <typename> typename ParticleInterface>
// template<typename T1, template <typename> typename T2> friend class Stack<T1, T2>;
using PI_C =
CombinedParticleInterface<ParticleInterfaceA, ParticleInterfaceB, StackIterator>;
using PI_A = ParticleInterfaceA<StackIterator>;
using PI_B = ParticleInterfaceB<ParticleInterfaceA<StackIterator>>;
protected:
using PI_B::GetIndex; // choose B, A would also work
using PI_B::GetStackData; // choose B, A would also work
public:
/**
* @name wrapper for user functions
* @{
*
* In this set of functions we call the user-provide
* ParticleInterface SetParticleData(...) methods, either with
* parent particle reference, or w/o.
*
* There is one implicit assumption here: if only one data tuple
* is provided for SetParticleData, the data is passed on to
* ParticleInterfaceA and the ParticleInterfaceB is
* default-initialized. There are many occasions where this is the
* desired behaviour, e.g. for thinning etc.
*
*/
template <typename... Args1>
void SetParticleData(const std::tuple<Args1...> vA) {
PI_A::SetParticleData(vA);
PI_B::SetParticleData();
}
template <typename... Args1, typename... Args2>
void SetParticleData(const std::tuple<Args1...> vA, const std::tuple<Args2...> vB) {
PI_A::SetParticleData(vA);
PI_B::SetParticleData(vB);
}
template <typename... Args1>
void SetParticleData(PI_C& p, const std::tuple<Args1...> vA) {
// static_assert(MT<I>::has_not, "error");
PI_A::SetParticleData(static_cast<PI_A&>(p), vA); // original stack
PI_B::SetParticleData(static_cast<PI_B&>(p)); // addon stack
}
template <typename... Args1, typename... Args2>
void SetParticleData(PI_C& p, const std::tuple<Args1...> vA,
const std::tuple<Args2...> vB) {
PI_A::SetParticleData(static_cast<PI_A&>(p), vA);
PI_B::SetParticleData(static_cast<PI_B&>(p), vB);
}
///@}
};
/**
* @class CombinedStackImpl
*
* Memory implementation of a combined data stack.
*
* The two stack data user objects Stack1Impl and Stack2Impl are
* merged into one consistent Stack container object providing
* access to the combined number of data entries.
*/
template <typename Stack1Impl, typename Stack2Impl>
class CombinedStackImpl : public Stack1Impl, public Stack2Impl {
public:
void Init() {
Stack1Impl::Init();
Stack2Impl::Init();
}
void Clear() {
Stack1Impl::Clear();
Stack2Impl::Clear();
}
unsigned int GetSize() const { return Stack1Impl::GetSize(); }
unsigned int GetCapacity() const { return Stack1Impl::GetCapacity(); }
/**
* Function to copy particle at location i1 in stack to i2
*/
void Copy(const unsigned int i1, const unsigned int i2) {
if (i1 >= GetSize() || i2 >= GetSize()) {
std::ostringstream err;
err << "CombinedStack: trying to access data beyond size of stack!";
throw std::runtime_error(err.str());
}
Stack1Impl::Copy(i1, i2);
Stack2Impl::Copy(i1, i2);
}
/**
* Function to copy particle at location i2 in stack to i1
*/
void Swap(const unsigned int i1, const unsigned int i2) {
if (i1 >= GetSize() || i2 >= GetSize()) {
std::ostringstream err;
err << "CombinedStack: trying to access data beyond size of stack!";
throw std::runtime_error(err.str());
}
Stack1Impl::Swap(i1, i2);
Stack2Impl::Swap(i1, i2);
}
void IncrementSize() {
Stack1Impl::IncrementSize();
Stack2Impl::IncrementSize();
}
void DecrementSize() {
Stack1Impl::DecrementSize();
Stack2Impl::DecrementSize();
}
}; // end class CombinedStackImpl
/**
* Helper template alias `CombinedStack` to construct new combined
* stack from two stack data objects and a particle readout interface.
*
* Note that the Stack2Impl provides only /additional/ data to
* Stack1Impl. This is important (see above) since tuple data for
* initialization are forwarded to Stack1Impl (first).
*/
template <typename Stack1Impl, typename Stack2Impl, template <typename> typename _PI>
using CombinedStack = Stack<CombinedStackImpl<Stack1Impl, Stack2Impl>, _PI>;
} // namespace corsika::stack
#endif
Framework/StackInterface/ParticleBase.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_particleBase_h_
#define _include_particleBase_h_
#include <type_traits>
namespace corsika::stack {
/**
@class ParticleBase
The base class to define the readout of particle properties from a
particle stack. Every stack must implement this readout via the
ParticleBase class.
The StackIterator template argument is derived from StackIteratorInterface, which is of
type <code> template <typename StackData, template <typename> typename
ParticleInterface> class StackIteratorInterface : public
ParticleInterface<StackIteratorInterface<StackData, ParticleInterface>>
</code>
where StackData must refer to a Stack type, and
ParticleInterface<StackIteratorInterface> is the corresponding particle readout class.
Thus, StackIteratorInterface is a CRTP class, injecting the full StackIteratorInterface
machinery into the ParticleInterface (aka ParticleBase) type!
The declartion of a StackIteratorInterface type simultaneously declares the
corresponding ParticleInterface type.
Furthermore, the operator* of the StackIteratorInterface returns a
static_cast to the ParticleInterface type, allowing a direct
readout of the particle data from the iterator.
*/
template <typename StackIterator>
class ParticleBase {
public:
using StackIteratorType = StackIterator;
ParticleBase() = default;
private:
// those copy constructors and assigments should never be implemented
ParticleBase(ParticleBase&) = delete;
ParticleBase operator=(ParticleBase&) = delete;
ParticleBase(ParticleBase&&) = delete;
ParticleBase operator=(ParticleBase&&) = delete;
ParticleBase(const ParticleBase&) = delete;
ParticleBase operator=(const ParticleBase&) = delete;
public:
/**
* Delete this particle on the stack. The corresponding iterator
* will be invalidated by this operation
*/
void Delete() { GetIterator().GetStack().Delete(GetIterator()); }
/**
* Add a secondary particle based on *this on the stack @param
* args is a variadic list of input data that has to match the
* function description in the user defined ParticleInterface::AddSecondary(...)
*/
template <typename... TArgs>
StackIterator AddSecondary(const TArgs... args) {
return GetStack().AddSecondary(GetIterator(), args...);
}
// protected: // todo should [MAY]be proteced, but don't now how to 'friend Stack'
// Function to provide CRTP access to inheriting class (type)
/**
* return the corresponding StackIterator for this particle
*/
StackIterator& GetIterator() { return static_cast<StackIterator&>(*this); }
const StackIterator& GetIterator() const {
return static_cast<const StackIterator&>(*this);
}
protected:
/**
@name Access to underlying stack fData, these are service
function for user classes. User code can only rely on GetIndex
and GetStackData to retrieve data
@{
*/
auto& GetStackData() { return GetIterator().GetStackData(); }
const auto& GetStackData() const { return GetIterator().GetStackData(); }
auto& GetStack() { return GetIterator().GetStack(); }
const auto& GetStack() const { return GetIterator().GetStack(); }
/**
* return the index number of the underlying iterator object
*/
unsigned int GetIndex() const { return GetIterator().GetIndexFromIterator(); }
///@}
};
} // namespace corsika::stack
#endif
Framework/StackInterface/SecondaryView.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_stack_secondaryview_h_
#define _include_corsika_stack_secondaryview_h_
#include <corsika/stack/Stack.h>
#include <stdexcept>
#include <vector>
namespace corsika::stack {
/**
* @class SecondaryView
*
* SecondaryView can only be constructed by giving a valid
* Projectile particle, following calls to AddSecondary will
* populate the original Stack, but will be directly accessible via
* the SecondaryView, e.g.
This allows to write code like
\verbatim
auto projectileInput = mainStack.GetNextParticle();
const unsigned int nMain = mainStack.GetSize();
SecondaryView<StackData, ParticleInterface> mainStackView(projectileInput);
mainStackView.AddSecondary(...data...);
mainStackView.AddSecondary(...data...);
mainStackView.AddSecondary(...data...);
mainStackView.AddSecondary(...data...);
assert(mainStackView.GetSize() == 4);
assert(mainStack.GetSize() = nMain+4);
\endverbatim
All operations possible on a Stack object are also possible on a
SecondaryView object. This means you can add, delete, copy, swap,
iterate, etc.
*/
/* INTERNAL NOTE FOR DEVELOPERS The secondary particle indices are
stored in a std::vector fIndices, the index of the primary projectle
particle is explicitly stored in fProjectileIndex. StackIterator
indices are refering to those numbers, where
StackIterator::GetIndex()==0 refers to the fProjectileIndex and
StackIterator::GetIndex()>0 to fIndices[i+1], see function GetIndexFromIterator.
*/
template <typename StackDataType, template <typename> typename ParticleInterface>
class SecondaryView : public Stack<StackDataType&, ParticleInterface> {
using ViewType = SecondaryView<StackDataType, ParticleInterface>;
private:
/**
* Helper type for inside this class
*/
using InnerStackType = Stack<StackDataType&, ParticleInterface>;
/**
* @name We need this "special" types with non-reference StackData for
* the constructor of the SecondaryView class
* @{
*/
using InnerStackTypeValue = Stack<StackDataType, ParticleInterface>;
using StackIteratorValue =
StackIteratorInterface<typename std::remove_reference<StackDataType>::type,
ParticleInterface, InnerStackTypeValue>;
/// @}
public:
using StackIterator =
StackIteratorInterface<typename std::remove_reference<StackDataType>::type,
ParticleInterface, ViewType>;
using ConstStackIterator =
ConstStackIteratorInterface<typename std::remove_reference<StackDataType>::type,
ParticleInterface, ViewType>;
/**
* this is the full type of the declared ParticleInterface: typedef typename
*/
using ParticleType = StackIterator;
using ParticleInterfaceType = typename StackIterator::ParticleInterfaceType;
friend class StackIteratorInterface<
typename std::remove_reference<StackDataType>::type, ParticleInterface, ViewType>;
friend class ConstStackIteratorInterface<
typename std::remove_reference<StackDataType>::type, ParticleInterface, ViewType>;
private:
/**
* This is not accessible, since we don't want to allow creating a
* new stack.
*/
template <typename... Args>
SecondaryView(Args... args) = delete;
public:
/**
SecondaryView can only be constructed passing it a valid
StackIterator to another Stack object
**/
SecondaryView(StackIteratorValue& vI)
: Stack<StackDataType&, ParticleInterface>(vI.GetStackData())
, fProjectileIndex(vI.GetIndex()) {}
StackIterator GetProjectile() {
// NOTE: 0 is special marker here for PROJECTILE, see GetIndexFromIterator
return StackIterator(*this, 0);
}
template <typename... Args>
auto AddSecondary(const Args... v) {
StackIterator proj = GetProjectile();
return AddSecondary(proj, v...);
}
template <typename... Args>
auto AddSecondary(StackIterator& proj, const Args... v) {
InnerStackType::GetStackData().IncrementSize();
const unsigned int idSec = GetSize();
const unsigned int index = InnerStackType::GetStackData().GetSize() - 1;
fIndices.push_back(index);
// NOTE: "+1" is since "0" is special marker here for PROJECTILE, see
// GetIndexFromIterator
return StackIterator(*this, idSec + 1, proj, v...);
}
/**
* overwrite Stack::GetSize to return actual number of secondaries
*/
unsigned int GetSize() const { return fIndices.size(); }
/**
* @name These are functions required by std containers and std loops
* The Stack-versions must be overwritten, since here we need the correct
* SecondaryView::GetSize
* @{
*/
// NOTE: the "+1" is since "0" is special marker here for PROJECTILE, see
// GetIndexFromIterator
auto begin() { return StackIterator(*this, 0 + 1); }
auto end() { return StackIterator(*this, GetSize() + 1); }
auto last() { return StackIterator(*this, GetSize() - 1 + 1); }
auto begin() const { return ConstStackIterator(*this, 0 + 1); }
auto end() const { return ConstStackIterator(*this, GetSize() + 1); }
auto last() const { return ConstStackIterator(*this, GetSize() - 1 + 1); }
auto cbegin() const { return ConstStackIterator(*this, 0 + 1); }
auto cend() const { return ConstStackIterator(*this, GetSize() + 1); }
auto clast() const { return ConstStackIterator(*this, GetSize() - 1 + 1); }
/// @}
/**
* need overwrite Stack::Delete, since we want to call SecondaryView::DeleteLast
*/
void Delete(StackIterator p) {
if (IsEmpty()) { /* error */
throw std::runtime_error("Stack, cannot delete entry since size is zero");
}
if (p.GetIndex() < GetSize() - 1)
InnerStackType::GetStackData().Copy(GetSize() - 1, p.GetIndex());
DeleteLast();
}
/**
* need overwrite Stack::Delete, since we want to call SecondaryView::DeleteLast
*/
void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
/**
* delete last particle on stack by decrementing stack size
*/
void DeleteLast() {
fIndices.pop_back();
InnerStackType::GetStackData().DecrementSize();
}
/**
* return next particle from stack, need to overwrtie Stack::GetNextParticle to get
* right reference
*/
StackIterator GetNextParticle() { return last(); }
/**
* check if there are no further particles on stack
*/
bool IsEmpty() { return GetSize() == 0; }
protected:
/**
* We only want to 'see' secondaries indexed in fIndices. In this
* function the conversion form iterator-index to stack-index is
* performed.
*/
unsigned int GetIndexFromIterator(const unsigned int vI) const {
if (vI == 0) return fProjectileIndex;
return fIndices[vI - 1];
}
private:
unsigned int fProjectileIndex;
std::vector<unsigned int> fIndices;
};
/*
See Issue 161
unfortunately clang does not support this in the same way (yet) as
gcc, so we have to distinguish here. If clang cataches up, we
could remove the #if here and elsewhere. The gcc code is much more
generic and universal.
*/
#if not defined(__clang__) && defined(__GNUC__) || defined(__GNUG__)
template <typename S, template <typename> typename _PIType = S::template PIType>
struct MakeView {
using type = corsika::stack::SecondaryView<typename S::StackImpl, _PIType>;
};
#endif
} // namespace corsika::stack
#endif
Framework/StackInterface/Stack.dox deleted 100644 → 0
View file @ 6f2ab230
/**
@page Stack Description of particle stacks
In the CORSIKA 8 framework particle data is always stored in
particle stacks. A particle is, thus, always a reference (iterator)
to an entry on a stack, e.g.
\verbatim
ModelStack stack;
stack.begin(); // returns an iterator: StackIterator, ConstStackIterator
*stack.begin(); // return a reference to ParticleInterfaceType, which is the class provided by the user to read/write particle properties
\endverbatim
All functionality and algorithms for stack data access is located in the namespace corsika::stack
The minimal example of how to use this is shown in stack_example.cc
*/
\ No newline at end of file
Framework/StackInterface/Stack.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_Stack_h__
#define _include_Stack_h__
#include <corsika/stack/StackIteratorInterface.h>
// must be after StackIteratorInterface
#include <corsika/stack/SecondaryView.h>
#include <corsika/utl/MetaProgramming.h>
#include <stdexcept>
#include <type_traits>
/**
All classes around management of particles on a stack.
*/
namespace corsika::stack {
/**
This is just a forward declatation for the user-defined
ParticleInterface, which is one of the essential template
parameters for the Stack.
<b>Important:</b> ParticleInterface must inherit from ParticleBase !
*/
template <typename>
class ParticleInterface;
/**
The Stack class provides (and connects) the main particle data storage machinery.
The StackDataType type is the user-provided bare data storage
object. This can be of any complexity, from a simple struct
(fortran common block), to a combination of different and
distributed data sources.
The user-provided ParticleInterface template type is the base
class type of the StackIteratorInterface class (CRTP) and must
provide all functions to read single particle data from the
StackDataType, given an 'unsigned int' index.
The Stack implements the
std-type begin/end function to allow integration in normal for
loops, ranges, etc.
*/
template <typename StackDataType, template <typename> typename ParticleInterface>
class Stack {
using StackDataValueType = std::remove_reference_t<StackDataType>;
StackDataType fData; ///< this in general holds all the data and can be quite big
private:
Stack(Stack&) = delete; ///< since Stack can be very big, we don't want to copy it
Stack& operator=(Stack&) =
delete; ///< since Stack can be very big, we don't want to copy it
public:
// Stack() { Init(); }
/**
* if StackDataType is a reference member we *HAVE* to initialize
* it in the constructor, this is typically needed for SecondaryView
*/
template <typename _ = StackDataType, typename = utl::enable_if<std::is_reference<_>>>
Stack(StackDataType vD)
: fData(vD) {}
/**
* This constructor takes any argument and passes it on to the
* StackDataType user class. If the user did not provide a suited
* constructor this will fail with an error message.
*
* Furthermore, this is disabled with enable_if for SecondaryView
* stacks, where the inner data container is always a reference
* and cannot be initialized here.
*/
template <typename... Args, typename _ = StackDataType,
typename = utl::disable_if<std::is_reference<_>>>
Stack(Args... args)
: fData(args...) {}
public:
typedef StackDataType
StackImpl; ///< this is the type of the user-provided data structure
template <typename SI>
using PIType = ParticleInterface<SI>;
/**
* Via the StackIteratorInterface and ConstStackIteratorInterface
* specialization, the type of the StackIterator
* template class is declared for a particular stack data
* object. Using CRTP, this also determines the type of
* ParticleInterface template class simultaneously.
*/
using StackIterator =
StackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
using ConstStackIterator =
ConstStackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
/**
* this is the full type of the user-declared ParticleInterface
*/
using ParticleInterfaceType = typename StackIterator::ParticleInterfaceType;
/**
* In all programming context, the object to access, copy, and
* transport particle data is via the StackIterator
*/
using ParticleType = StackIterator;
// friends are needed since they need access to protected members
friend class StackIteratorInterface<StackDataValueType, ParticleInterface, Stack>;
friend class ConstStackIteratorInterface<StackDataValueType, ParticleInterface,
Stack>;
public:
/**
* @name Most generic proxy methods for StackDataType fData
* @{
*/
unsigned int GetCapacity() const { return fData.GetCapacity(); }
unsigned int GetSize() const { return fData.GetSize(); }
template <typename... Args>
auto Init(Args... args) {
return fData.Init(args...);
}
template <typename... Args>
auto Clear(Args... args) {
return fData.Clear(args...);
}
///@}
public:
/**
* @name These are functions required by std containers and std loops
* @{
*/
StackIterator begin() { return StackIterator(*this, 0); }
StackIterator end() { return StackIterator(*this, GetSize()); }
StackIterator last() { return StackIterator(*this, GetSize() - 1); }
ConstStackIterator begin() const { return ConstStackIterator(*this, 0); }
ConstStackIterator end() const { return ConstStackIterator(*this, GetSize()); }
ConstStackIterator last() const { return ConstStackIterator(*this, GetSize() - 1); }
ConstStackIterator cbegin() const { return ConstStackIterator(*this, 0); }
ConstStackIterator cend() const { return ConstStackIterator(*this, GetSize()); }
ConstStackIterator clast() const { return ConstStackIterator(*this, GetSize() - 1); }
/// @}
/**
* increase stack size, create new particle at end of stack
*/
template <typename... Args>
StackIterator AddParticle(const Args... v) {
fData.IncrementSize();
return StackIterator(*this, GetSize() - 1, v...);
}
/**
* increase stack size, create new particle at end of stack, related to parent
* particle/projectile
*/
template <typename... Args>
StackIterator AddSecondary(StackIterator& parent, const Args... v) {
fData.IncrementSize();
return StackIterator(*this, GetSize() - 1, parent, v...);
}
void Swap(StackIterator a, StackIterator b) {
fData.Swap(a.GetIndex(), b.GetIndex());
}
void Swap(ConstStackIterator a, ConstStackIterator b) {
fData.Swap(a.GetIndex(), b.GetIndex());
}
void Copy(StackIterator a, StackIterator b) {
fData.Copy(a.GetIndex(), b.GetIndex());
}
void Copy(ConstStackIterator a, StackIterator b) {
fData.Copy(a.GetIndex(), b.GetIndex());
}
/**
* delete this particle
*/
void Delete(StackIterator p) {
if (GetSize() == 0) { /*error*/
throw std::runtime_error("Stack, cannot delete entry since size is zero");
}
if (p.GetIndex() < GetSize() - 1) fData.Copy(GetSize() - 1, p.GetIndex());
DeleteLast();
}
/**
* delete this particle
*/
void Delete(ParticleInterfaceType p) { Delete(p.GetIterator()); }
/**
* delete last particle on stack by decrementing stack size
*/
void DeleteLast() { fData.DecrementSize(); }
/**
* check if there are no further particles on stack
*/
bool IsEmpty() { return GetSize() == 0; }
/**
* return next particle from stack
*/
StackIterator GetNextParticle() { return last(); }
protected:
/**
* Function to perform eventual transformation from
* StackIterator::GetIndex() to index in data stored in
* StackDataType fData. By default (and in almost all cases) this
* should just be identiy. See class SecondaryView for an alternative implementation.
*/
unsigned int GetIndexFromIterator(const unsigned int vI) const { return vI; }
/**
* @name Return reference to StackDataType object fData for data access
* @{
*/
StackDataValueType& GetStackData() { return fData; }
const StackDataValueType& GetStackData() const { return fData; }
///@}
};
} // namespace corsika::stack
#endif
Framework/StackInterface/StackIteratorInterface.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_StackIteratorinterface_h__
#define _include_StackIteratorinterface_h__
#include <corsika/stack/ParticleBase.h>
namespace corsika::stack {
template <typename StackDataType, template <typename> typename ParticleInterface>
class Stack; // forward decl
template <typename StackDataType, template <typename> typename ParticleInterface>
class SecondaryView; // forward decl
/**
@class StackIteratorInterface
The StackIteratorInterface is the main interface to iterator over
particles on a stack. At the same time StackIteratorInterface is a
Particle object by itself, thus there is no difference between
type and ref_type for convenience of the physicist.
This allows to write code like
\verbatim
for (auto& p : theStack) { p.SetEnergy(newEnergy); }
\endverbatim
The template argument Stack determines the type of Stack object
the data is stored in. A pointer to the Stack object is part of
the StackIteratorInterface. In addition to Stack the iterator only knows
the index fIndex in the Stack data.
The template argument `ParticleInterface` acts as a policy to provide
readout function of Particle data from the stack. The ParticleInterface
class must know how to retrieve information from the Stack data
for a particle entry at any index fIndex.
The ParticleInterface class must be written and provided by the
user, it contains methods like <code> auto GetData() const {
return GetStackData().GetData(GetIndex()); }</code>, where
StackIteratorInterface::GetStackData() return a reference to the
object storing the particle data of type StackDataType. And
StackIteratorInterface::GetIndex() provides the iterator index to
be readout. The StackDataType is another user-provided class to
store data and must implement functions compatible with
ParticleInterface, in this example StackDataType::GetData(const unsigned int
vIndex).
For two examples see stack_example.cc, or the
corsika::processes::sibyll::SibStack class
*/
template <typename StackDataType, template <typename> typename ParticleInterface,
typename StackType = Stack<StackDataType, ParticleInterface>>
class StackIteratorInterface
: public ParticleInterface<
StackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
public:
using ParticleInterfaceType =
ParticleInterface<corsika::stack::StackIteratorInterface<
StackDataType, ParticleInterface, StackType>>;
// friends are needed for access to protected methods
friend class Stack<StackDataType,
ParticleInterface>; // for access to GetIndex for Stack
friend class Stack<StackDataType&, ParticleInterface>; // for access to GetIndex
// SecondaryView : public Stack
friend class ParticleBase<StackIteratorInterface>; // for access to GetStackDataType
friend class SecondaryView<StackDataType,
ParticleInterface>; // access for SecondaryView
private:
unsigned int fIndex = 0;
StackType* fData = 0; // info: Particles and StackIterators become invalid when parent
// Stack is copied or deleted!
// it is not allowed to create a "dangling" stack iterator
StackIteratorInterface() = delete;
public:
StackIteratorInterface(StackIteratorInterface const& vR)
: fIndex(vR.fIndex)
, fData(vR.fData) {}
StackIteratorInterface& operator=(StackIteratorInterface const& vR) {
fIndex = vR.fIndex;
fData = vR.fData;
return *this;
}
/** iterator must always point to data, with an index:
@param data reference to the stack [rw]
@param index index on stack
*/
StackIteratorInterface(StackType& data, const unsigned int index)
: fIndex(index)
, fData(&data) {}
/** constructor that also sets new values on particle data object
@param data reference to the stack [rw]
@param index index on stack
@param args variadic list of data to initialize stack entry, this must be
consistent with the definition of the user-provided
ParticleInterfaceType::SetParticleData(...) function
*/
template <typename... Args>
StackIteratorInterface(StackType& data, const unsigned int index, const Args... args)
: fIndex(index)
, fData(&data) {
(**this).SetParticleData(args...);
}
/** constructor that also sets new values on particle data object, including reference
to parent particle
@param data reference to the stack [rw]
@param index index on stack
@param reference to parent particle [rw]. This can be used for thinning, particle
counting, history, etc.
@param args variadic list of data to initialize stack entry, this must be
consistent with the definition of the user-provided
ParticleInterfaceType::SetParticleData(...) function
*/
template <typename... Args>
StackIteratorInterface(StackType& data, const unsigned int index,
StackIteratorInterface& parent, const Args... args)
: fIndex(index)
, fData(&data) {
(**this).SetParticleData(*parent, args...);
}
public:
/** @name Iterator interface
@{
*/
StackIteratorInterface& operator++() {
++fIndex;
return *this;
}
StackIteratorInterface operator++(int) {
StackIteratorInterface tmp(*this);
++fIndex;
return tmp;
}
StackIteratorInterface operator+(int delta) {
return StackIteratorInterface(*fData, fIndex + delta);
}
bool operator==(const StackIteratorInterface& rhs) { return fIndex == rhs.fIndex; }
bool operator!=(const StackIteratorInterface& rhs) { return fIndex != rhs.fIndex; }
/**
* Convert iterator to value type, where value type is the user-provided particle
* readout class
*/
ParticleInterfaceType& operator*() {
return static_cast<ParticleInterfaceType&>(*this);
}
/**
* Convert iterator to const value type, where value type is the user-provided
* particle readout class
*/
const ParticleInterfaceType& operator*() const {
return static_cast<const ParticleInterfaceType&>(*this);
}
///@}
protected:
/**
* @name Stack data access
* @{
*/
/// Get current particle index
inline unsigned int GetIndex() const { return fIndex; }
/// Get current particle Stack object
inline StackType& GetStack() { return *fData; }
/// Get current particle const Stack object
inline const StackType& GetStack() const { return *fData; }
/// Get current user particle StackDataType object
inline StackDataType& GetStackData() { return fData->GetStackData(); }
/// Get current const user particle StackDataType object
inline const StackDataType& GetStackData() const { return fData->GetStackData(); }
/// Get data index as mapped in Stack class
inline unsigned int GetIndexFromIterator() const {
return fData->GetIndexFromIterator(fIndex);
}
///@}
}; // end class StackIterator
/**
@class ConstStackIteratorInterface
This is the iterator class for const-access to stack data
*/
template <typename StackDataType, template <typename> typename ParticleInterface,
typename StackType = Stack<StackDataType, ParticleInterface>>
class ConstStackIteratorInterface
: public ParticleInterface<
ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>> {
public:
typedef ParticleInterface<
ConstStackIteratorInterface<StackDataType, ParticleInterface, StackType>>
ParticleInterfaceType;
friend class Stack<StackDataType, ParticleInterface>; // for access to GetIndex
friend class ParticleBase<ConstStackIteratorInterface>; // for access to
// GetStackDataType
private:
unsigned int fIndex = 0;
const StackType* fData = 0; // info: Particles and StackIterators become invalid when
// parent Stack is copied or deleted!
// we don't want to allow dangling iterators to exist
ConstStackIteratorInterface() = delete;
public:
ConstStackIteratorInterface(const StackType& data, const unsigned int index)
: fIndex(index)
, fData(&data) {}
/**
@class ConstStackIteratorInterface
The const counterpart of StackIteratorInterface, which is used
for read-only iterator access on particle stack:
\verbatim
for (const auto& p : theStack) { E += p.GetEnergy(); }
\endverbatim
See documentation of StackIteratorInterface for more details.
*/
public:
/** @name Iterator interface
*/
///@{
ConstStackIteratorInterface& operator++() {
++fIndex;
return *this;
}
ConstStackIteratorInterface operator++(int) {
ConstStackIteratorInterface tmp(*this);
++fIndex;
return tmp;
}
ConstStackIteratorInterface operator+(int delta) {
return ConstStackIteratorInterface(*fData, fIndex + delta);
}
bool operator==(const ConstStackIteratorInterface& rhs) {
return fIndex == rhs.fIndex;
}
bool operator!=(const ConstStackIteratorInterface& rhs) {
return fIndex != rhs.fIndex;
}
const ParticleInterfaceType& operator*() const {
return static_cast<const ParticleInterfaceType&>(*this);
}
///@}
protected:
/** @name Stack data access
Only the const versions for read-only access
*/
///@{
inline unsigned int GetIndex() const { return fIndex; }
inline const StackType& GetStack() const { return *fData; }
inline const StackDataType& GetStackData() const { return fData->GetStackData(); }
/// Get data index as mapped in Stack class
inline unsigned int GetIndexFromIterator() const {
return fData->GetIndexFromIterator(fIndex);
}
///@}
}; // end class ConstStackIterator
} // namespace corsika::stack
#endif
Framework/StackInterface/comp deleted 100644 → 0
View file @ 6f2ab230
TestParticleInterface2<TestParticleInterface<corsika::stack::StackIteratorInterface<corsika::stack::CombinedStackImpl<corsika::stack::CombinedStackImpl<TestStackData, TestStackData2>, TestStackData3>, CombinedTestInterfaceType2, corsika::stack::Stack<corsika::stack::CombinedStackImpl<corsika::stack::CombinedStackImpl<TestStackData, TestStackData2>, TestStackData3>, CombinedTestInterfaceType2> > > >::content
Framework/StackInterface/testCombinedStack.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/stack/CombinedStack.h>
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/Stack.h>
#include <testTestStack.h> // for testing: simple stack. This is a
// test-build, and inluce file is obtained from CMAKE_CURRENT_SOURCE_DIR
#include <boost/type_index.hpp>
#include <type_traits>
using boost::typeindex::type_id_with_cvr;
#include <iomanip>
#include <iostream>
#include <vector>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::stack;
using namespace std;
////////////////////////////////////////////////////////////
// first level test: combine two stacks:
// StackTest = (TestStackData + TestStackData2)
// definition of stack-data object
class TestStackData2 {
public:
// these functions are needed for the Stack interface
void Init() {}
void Clear() { fData2.clear(); }
unsigned int GetSize() const { return fData2.size(); }
unsigned int GetCapacity() const { return fData2.size(); }
void Copy(const int i1, const int i2) { fData2[i2] = fData2[i1]; }
void Swap(const int i1, const int i2) {
double tmp0 = fData2[i1];
fData2[i1] = fData2[i2];
fData2[i2] = tmp0;
}
// custom data access function
void SetData2(const int i, const double v) { fData2[i] = v; }
double GetData2(const int i) const { return fData2[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fData2.push_back(0.); }
void DecrementSize() {
if (fData2.size() > 0) { fData2.pop_back(); }
}
// custom private data section
private:
std::vector<double> fData2;
};
// defintion of a stack-readout object, the iteractor dereference
// operator will deliver access to these function
template <typename T>
class TestParticleInterface2 : public T {
public:
using T::GetIndex;
using T::GetStackData;
using T::SetParticleData;
// default version for particle-creation from input data
void SetParticleData(const std::tuple<double> v = {0.}) { SetData2(std::get<0>(v)); }
void SetParticleData(TestParticleInterface2<T>& parent,
const std::tuple<double> v = {0.}) {
SetData2(parent.GetData2() + std::get<0>(v));
}
void SetData2(const double v) { GetStackData().SetData2(GetIndex(), v); }
double GetData2() const { return GetStackData().GetData2(GetIndex()); }
};
// combined stack: StackTest = (TestStackData + TestStackData2)
template <typename StackIter>
using CombinedTestInterfaceType =
corsika::stack::CombinedParticleInterface<TestParticleInterface,
TestParticleInterface2, StackIter>;
using StackTest = CombinedStack<TestStackData, TestStackData2, CombinedTestInterfaceType>;
TEST_CASE("Combined Stack", "[stack]") {
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
double v = 0;
for (const auto& p : stack) v += p.GetData();
return v;
};
auto sum2 = [](const StackTest& stack) {
double v = 0;
for (const auto& p : stack) v += p.GetData2();
return v;
};
SECTION("StackInterface") {
// construct a valid Stack object
StackTest s;
s.Init();
s.Clear();
s.AddParticle(std::tuple{0.});
s.Copy(s.cbegin(), s.begin());
s.Swap(s.begin(), s.begin());
REQUIRE(s.GetSize() == 1);
}
SECTION("construct") {
// construct a valid, empty Stack object
StackTest s;
}
SECTION("write and read") {
StackTest s;
s.AddParticle(std::tuple{9.9});
REQUIRE(sum2(s) == 0.);
REQUIRE(sum(s) == 9.9);
}
SECTION("delete from stack") {
StackTest s;
REQUIRE(s.GetSize() == 0);
StackTest::StackIterator p =
s.AddParticle(std::tuple{0.}); // valid way to access particle data
p.SetData(8.9);
p.SetData2(3.);
REQUIRE(sum2(s) == 3.);
REQUIRE(sum(s) == 8.9);
REQUIRE(s.GetSize() == 1);
s.Delete(p);
REQUIRE(s.GetSize() == 0);
}
SECTION("delete particle") {
StackTest s;
REQUIRE(s.GetSize() == 0);
auto p = s.AddParticle(
std::tuple{9.9}); // also valid way to access particle data, identical to above
REQUIRE(s.GetSize() == 1);
p.Delete();
REQUIRE(s.GetSize() == 0);
}
SECTION("create secondaries") {
StackTest s;
REQUIRE(s.GetSize() == 0);
auto iter = s.AddParticle(std::tuple{9.9});
iter.SetData2(2);
REQUIRE(s.GetSize() == 1);
iter.AddSecondary(std::tuple{4.4});
REQUIRE(s.GetSize() == 2);
// p.AddSecondary(3.3, 2.2, 1.);
// REQUIRE(s.GetSize() == 3);
double v = 0;
for (const auto& i : s) {
v += i.GetData();
REQUIRE(i.GetData2() == 2);
}
REQUIRE(v == 9.9 + 4.4);
}
SECTION("get next particle") {
StackTest s;
REQUIRE(s.GetSize() == 0);
auto p1 = s.AddParticle(std::tuple{9.9});
auto p2 = s.AddParticle(std::tuple{8.8});
p1.SetData2(20.2);
p2.SetData2(20.3);
auto particle = s.GetNextParticle(); // first particle
REQUIRE(particle.GetData() == 8.8);
REQUIRE(particle.GetData2() == 20.3);
particle.Delete();
auto particle2 = s.GetNextParticle(); // first particle
REQUIRE(particle2.GetData() == 9.9);
REQUIRE(particle2.GetData2() == 20.2);
particle2.Delete();
REQUIRE(s.GetSize() == 0);
}
}
////////////////////////////////////////////////////////////
// next level: combine three stacks:
// combined stack: StackTest2 = ((TestStackData + TestStackData2) + TestStackData3)
// definition of stack-data object
class TestStackData3 {
public:
// these functions are needed for the Stack interface
void Init() {}
void Clear() { fData3.clear(); }
unsigned int GetSize() const { return fData3.size(); }
unsigned int GetCapacity() const { return fData3.size(); }
void Copy(const int i1, const int i2) { fData3[i2] = fData3[i1]; }
void Swap(const int i1, const int i2) {
double tmp0 = fData3[i1];
fData3[i1] = fData3[i2];
fData3[i2] = tmp0;
}
// custom data access function
void SetData3(const int i, const double v) { fData3[i] = v; }
double GetData3(const int i) const { return fData3[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fData3.push_back(0.); }
void DecrementSize() {
if (fData3.size() > 0) { fData3.pop_back(); }
}
// custom private data section
private:
std::vector<double> fData3;
};
// ---------------------------------------
// defintion of a stack-readout object, the iteractor dereference
// operator will deliver access to these function
template <typename T>
class TestParticleInterface3 : public T {
public:
using T::GetIndex;
using T::GetStackData;
using T::SetParticleData;
// default version for particle-creation from input data
void SetParticleData(const std::tuple<double> v = {0.}) { SetData3(std::get<0>(v)); }
void SetParticleData(TestParticleInterface3<T>& parent,
const std::tuple<double> v = {0.}) {
SetData3(parent.GetData3() + std::get<0>(v));
}
void SetData3(const double v) { GetStackData().SetData3(GetIndex(), v); }
double GetData3() const { return GetStackData().GetData3(GetIndex()); }
};
// double combined stack:
// combined stack
template <typename StackIter>
using CombinedTestInterfaceType2 =
corsika::stack::CombinedParticleInterface<StackTest::PIType, TestParticleInterface3,
StackIter>;
using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
CombinedTestInterfaceType2>;
TEST_CASE("Combined Stack - multi", "[stack]") {
SECTION("create secondaries") {
StackTest2 s;
REQUIRE(s.GetSize() == 0);
// add new particle, only provide tuple data for StackTest
auto p1 = s.AddParticle(std::tuple{9.9});
// add new particle, provide tuple data for both StackTest and TestStackData3
auto p2 = s.AddParticle(std::tuple{8.8}, std::tuple{0.1});
// examples to explicitly change data on stack
p2.SetData2(0.1); // not clear why this is needed, need to check
// SetParticleData workflow for more complicated
// settings
p1.SetData3(20.2);
p2.SetData3(10.3);
REQUIRE(p1.GetData() == 9.9);
REQUIRE(p1.GetData2() == 0.);
p1.SetData2(10.2);
REQUIRE(p1.GetData2() == 10.2);
REQUIRE(p1.GetData3() == 20.2);
REQUIRE(p2.GetData() == 8.8);
REQUIRE(p2.GetData2() == 0.1);
REQUIRE(p2.GetData3() == 10.3);
auto particle = s.GetNextParticle(); // first particle
REQUIRE(particle.GetData() == 8.8);
REQUIRE(particle.GetData2() == 0.1);
REQUIRE(particle.GetData3() == 10.3);
REQUIRE(s.GetSize() == 2);
auto sec = particle.AddSecondary(std::tuple{4.4});
REQUIRE(s.GetSize() == 3);
REQUIRE(sec.GetData() == 4.4);
REQUIRE(sec.GetData2() == 0.1);
REQUIRE(sec.GetData3() == 10.3);
sec.Delete();
s.DeleteLast();
s.GetNextParticle().Delete();
REQUIRE(s.GetSize() == 0);
}
}
////////////////////////////////////////////////////////////
// final level test, create SecondaryView on StackTest2
/*
See Issue 161
unfortunately clang does not support this in the same way (yet) as
gcc, so we have to distinguish here. If clang cataches up, we could
remove the clang branch here and also in corsika::Cascade. The gcc
code is much more generic and universal.
*/
template <typename StackIter>
using CombinedTestInterfaceType2 =
corsika::stack::CombinedParticleInterface<StackTest::PIType, TestParticleInterface3,
StackIter>;
using StackTest2 = CombinedStack<typename StackTest::StackImpl, TestStackData3,
CombinedTestInterfaceType2>;
#if defined(__clang__)
using StackTestView = SecondaryView<typename StackTest2::StackImpl, CombinedTestInterfaceType2>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackTestView = corsika::stack::MakeView<StackTest2>::type;
#endif
using Particle2 = typename StackTest2::ParticleType;
TEST_CASE("Combined Stack - secondary view") {
SECTION("create secondaries via secondaryview") {
StackTest2 stack;
auto particle = stack.AddParticle(std::tuple{9.9});
// cout << boost::typeindex::type_id_runtime(particle).pretty_name() << endl;
StackTestView view(particle);
auto projectile = view.GetProjectile();
projectile.AddSecondary(std::tuple{8.8});
REQUIRE(stack.GetSize() == 2);
}
}
Framework/StackInterface/testSecondaryView.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/Stack.h>
#include <testTestStack.h> // for testing: simple stack. This is a
// test-build, and inluce file is obtained from CMAKE_CURRENT_SOURCE_DIR
#include <boost/type_index.hpp>
#include <type_traits>
using boost::typeindex::type_id_with_cvr;
#include <iomanip>
#include <iostream>
#include <vector>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::stack;
using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest;
/*
See Issue 161
unfortunately clang does not support this in the same way (yet) as
gcc, so we have to distinguish here. If clang cataches up, we could
remove the clang branch here and also in corsika::Cascade. The gcc
code is much more generic and universal.
*/
#if defined(__clang__)
using StackTestView = SecondaryView<TestStackData, TestParticleInterface>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackTestView = MakeView<StackTest>::type;
#endif
using Particle = typename StackTest::ParticleType;
TEST_CASE("SecondaryStack", "[stack]") {
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
double v = 0;
for (const auto& p : stack) v += p.GetData();
return v;
};
SECTION("secondary view") {
StackTest s;
REQUIRE(s.GetSize() == 0);
s.AddParticle(std::tuple{9.9});
s.AddParticle(std::tuple{8.8});
const double sumS = 9.9 + 8.8;
auto particle = s.GetNextParticle();
StackTestView view(particle);
REQUIRE(view.GetSize() == 0);
{
auto proj = view.GetProjectile();
REQUIRE(proj.GetData() == particle.GetData());
proj.AddSecondary(std::tuple{4.4});
}
view.AddSecondary(std::tuple{4.5});
view.AddSecondary(std::tuple{4.6});
REQUIRE(view.GetSize() == 3);
REQUIRE(s.GetSize() == 5);
REQUIRE(!view.IsEmpty());
auto sumView = [](const StackTestView& stack) {
double value = 0;
for (const auto& p : stack) { value += p.GetData(); }
return value;
};
REQUIRE(sum(s) == sumS + 4.4 + 4.5 + 4.6);
REQUIRE(sumView(view) == 4.4 + 4.5 + 4.6);
view.DeleteLast();
REQUIRE(view.GetSize() == 2);
REQUIRE(s.GetSize() == 4);
REQUIRE(sum(s) == sumS + 4.4 + 4.5);
REQUIRE(sumView(view) == 4.4 + 4.5);
auto pDel = view.GetNextParticle();
view.Delete(pDel);
REQUIRE(view.GetSize() == 1);
REQUIRE(s.GetSize() == 3);
REQUIRE(sum(s) == sumS + 4.4 + 4.5 - pDel.GetData());
REQUIRE(sumView(view) == 4.4 + 4.5 - pDel.GetData());
view.Delete(view.GetNextParticle());
REQUIRE(sum(s) == sumS);
REQUIRE(sumView(view) == 0);
REQUIRE(view.IsEmpty());
{
auto proj = view.GetProjectile();
REQUIRE(proj.GetData() == particle.GetData());
}
}
SECTION("secondary view, construct from ParticleType") {
StackTest s;
REQUIRE(s.GetSize() == 0);
s.AddParticle(std::tuple{9.9});
s.AddParticle(std::tuple{8.8});
auto iterator = s.GetNextParticle();
typename StackTest::ParticleType& particle = iterator; // as in corsika::Cascade
StackTestView view(particle);
REQUIRE(view.GetSize() == 0);
view.AddSecondary(std::tuple{4.4});
REQUIRE(view.GetSize() == 1);
}
}
Framework/StackInterface/testStackInterface.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/stack/Stack.h>
#include <testTestStack.h> // simple test-stack for testing. This is
// for testing only: include from
// CMAKE_CURRENT_SOURCE_DIR
#include <boost/type_index.hpp>
#include <type_traits>
using boost::typeindex::type_id_with_cvr;
#include <iomanip>
#include <iostream>
#include <tuple>
#include <vector>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::stack;
using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest;
TEST_CASE("Stack", "[Stack]") {
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
double v = 0;
for (const auto& p : stack) v += p.GetData();
return v;
};
SECTION("StackInterface") {
// construct a valid Stack object
StackTest s;
s.Clear();
s.AddParticle(std::tuple{0.});
s.Copy(s.cbegin(), s.begin());
s.Swap(s.begin(), s.begin());
REQUIRE(s.GetSize() == 1);
}
SECTION("construct") {
// construct a valid, empty Stack object
StackTest s;
}
SECTION("write and read") {
StackTest s;
s.AddParticle(std::tuple{9.9});
const double v = sum(s);
REQUIRE(v == 9.9);
}
SECTION("delete from stack") {
StackTest s;
REQUIRE(s.GetSize() == 0);
StackTest::StackIterator p =
s.AddParticle(std::tuple{0.}); // valid way to access particle data
p.SetData(9.9);
REQUIRE(s.GetSize() == 1);
s.Delete(p);
REQUIRE(s.GetSize() == 0);
}
SECTION("delete particle") {
StackTest s;
REQUIRE(s.GetSize() == 0);
auto p = s.AddParticle(
std::tuple{9.9}); // also valid way to access particle data, identical to above
REQUIRE(s.GetSize() == 1);
p.Delete();
REQUIRE(s.GetSize() == 0);
}
SECTION("create secondaries") {
StackTest s;
REQUIRE(s.GetSize() == 0);
auto iter = s.AddParticle(std::tuple{9.9});
StackTest::ParticleInterfaceType& p =
*iter; // also this is valid to access particle data
REQUIRE(s.GetSize() == 1);
p.AddSecondary(std::tuple{4.4});
REQUIRE(s.GetSize() == 2);
/*p.AddSecondary(3.3, 2.2);
REQUIRE(s.GetSize() == 3);
double v = 0;
for (auto& p : s) { v += p.GetData(); }
REQUIRE(v == 9.9 + 4.4 + 3.3 + 2.2);*/
}
SECTION("get next particle") {
StackTest s;
REQUIRE(s.GetSize() == 0);
s.AddParticle(std::tuple{9.9});
s.AddParticle(std::tuple{8.8});
auto particle = s.GetNextParticle(); // first particle
REQUIRE(particle.GetData() == 8.8);
particle.Delete();
auto particle2 = s.GetNextParticle(); // first particle
REQUIRE(particle2.GetData() == 9.9);
particle2.Delete();
REQUIRE(s.GetSize() == 0);
}
}
Framework/StackInterface/testTestStack.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_stack_testTestStack_h_
#define _include_corsika_stack_testTestStack_h_
#include <corsika/stack/Stack.h>
#include <tuple>
#include <vector>
/**
* definition of stack-data object for unit tests
*
* TestStackData contain only a single variable "Data" stored in fData
* with Get/SetData functions.
*/
class TestStackData {
public:
// these functions are needed for the Stack interface
void Init() {}
void Clear() { fData.clear(); }
unsigned int GetSize() const { return fData.size(); }
unsigned int GetCapacity() const { return fData.size(); }
void Copy(const unsigned int i1, const unsigned int i2) { fData[i2] = fData[i1]; }
void Swap(const unsigned int i1, const unsigned int i2) {
double tmp0 = fData[i1];
fData[i1] = fData[i2];
fData[i2] = tmp0;
}
// custom data access function
void SetData(const unsigned int i, const double v) { fData[i] = v; }
double GetData(const unsigned int i) const { return fData[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fData.push_back(0.); }
void DecrementSize() {
if (fData.size() > 0) { fData.pop_back(); }
}
// custom private data section
private:
std::vector<double> fData;
};
/**
* From static_cast of a StackIterator over entries in the
* TestStackData class you get and object of type
* TestParticleInterface defined here
*
* It provides Get/Set methods to read and write data to the "Data"
* storage of TestStackData obtained via
* "StackIteratorInterface::GetStackData()", given the index of the
* iterator "StackIteratorInterface::GetIndex()"
*
*/
template <typename StackIteratorInterface>
class TestParticleInterface
: public corsika::stack::ParticleBase<StackIteratorInterface> {
public:
using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
/*
The SetParticleData methods are called for creating new entries
on the stack. You can specifiy various parametric versions to
perform this task:
*/
// default version for particle-creation from input data
void SetParticleData(const std::tuple<double> v) { SetData(std::get<0>(v)); }
void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
std::tuple<double> v) {
SetData(std::get<0>(v));
}
/// alternative set-particle data for non-standard construction from different inputs
/*
void SetParticleData(const double v, const double p) { SetData(v + p); }
void SetParticleData(TestParticleInterface<StackIteratorInterface>&,
const double v, const double p) {
SetData(v + p);
}*/
// here are the fundamental methods for access to TestStackData data
void SetData(const double v) { GetStackData().SetData(GetIndex(), v); }
double GetData() const { return GetStackData().GetData(GetIndex()); }
};
#endif
Framework/Units/CMakeLists.txt deleted 100644 → 0
View file @ 6f2ab230
add_library (CORSIKAunits INTERFACE)
set (CORSIKAunits_NAMESPACE corsika/units)
set (
CORSIKAunits_HEADERS
PhysicalUnits.h
PhysicalConstants.h
)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAunits ${CORSIKAunits_NAMESPACE} ${CORSIKAunits_HEADERS})
target_include_directories (CORSIKAunits
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/ThirdParty>
$<INSTALL_INTERFACE:include>
)
install (FILES ${CORSIKAunits_HEADERS} DESTINATION include/${CORSIKAunits_NAMESPACE})
# code testing
add_executable (testUnits testUnits.cc)
target_link_libraries (testUnits CORSIKAunits CORSIKAthirdparty) # for catch2
CORSIKA_ADD_TEST(testUnits)
Framework/Units/PhysicalConstants.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef INCLUDE_PHYSICAL_CONSTANTS_H
#define INCLUDE_PHYSICAL_CONSTANTS_H
#include <phys/units/quantity.hpp>
namespace corsika::units::constants {
using namespace phys::units;
// acceleration of free-fall, standard
constexpr phys::units::quantity<phys::units::acceleration_d> g_sub_n{
phys::units::Rep(9.80665L) * phys::units::meter /
phys::units::square(phys::units::second)};
// Avogadro constant
constexpr quantity<dimensions<0, 0, 0, 0, 0, -1>> N_sub_A{Rep(6.02214199e+23L) / mole};
// elementary charge
constexpr quantity<electric_charge_d> e{Rep(1.6021766208e-19L) * coulomb};
// electronvolt
// constexpr quantity<hepenergy_d> eV{e / coulomb * joule};
// Planck constant
constexpr quantity<dimensions<2, 1, -1>> h{Rep(6.62606876e-34L) * joule * second};
constexpr quantity<dimensions<2, 1, -1>> hBar{h / (2 * M_PI)};
// speed of light in a vacuum
constexpr quantity<speed_d> c{Rep(299792458L) * meter / second};
constexpr auto cSquared = c * c;
// hbar * c
constexpr quantity<dimensions<1, 0, 0, 0, 0, 0, 0, 1>> hBarC{
Rep(1.973'269'78e-7L) * electronvolt * meter}; // from RPP 2018
auto constexpr invGeVsq = 1e-18 / (electronvolt * electronvolt);
// unified atomic mass unit
constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
auto constexpr nucleonMass = 0.5 * (0.93827 + 0.93957) * 1e9 * electronvolt;
// etc.
} // namespace corsika::units::constants
#endif // PHYS_UNITS_PHYSICAL_CONSTANTS_HPP_INCLUDED
Framework/Units/PhysicalUnits.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_PhysicalUnits_h_
#define _include_PhysicalUnits_h_
#include <corsika/units/PhysicalConstants.h>
#include <phys/units/io.hpp>
#include <phys/units/quantity.hpp>
/*
It is essentially a bug of the phys/units package to define the
operator<< not in the same namespace as the types it is working
on. This breaks ADL (argument-dependent lookup). Here we "fix" this:
*/
namespace phys::units {
using phys::units::io::operator<<;
} // namespace phys::units
/**
* @file PhysicalUnits
*
* Add new units and types we need. Units are compile-time. Literals are
* used for optimal coding style.
*
*/
namespace corsika::units::si {
using namespace phys::units;
using namespace phys::units::literals;
using namespace phys::units::io;
using phys::units::io::operator<<;
/// defining momentum you suckers
/// dimensions, i.e. composition in base SI dimensions
using hepmomentum_d = phys::units::hepenergy_d;
using hepmass_d = phys::units::hepenergy_d;
/// defining cross section as area
using sigma_d = phys::units::area_d;
// constexpr quantity<area_d> barn{Rep(1e-28L) * square(meter)};
/// add the unit-types
using LengthType = phys::units::quantity<phys::units::length_d, double>;
using TimeType = phys::units::quantity<phys::units::time_interval_d, double>;
using SpeedType = phys::units::quantity<phys::units::speed_d, double>;
using FrequencyType = phys::units::quantity<phys::units::frequency_d, double>;
using ElectricChargeType =
phys::units::quantity<phys::units::electric_charge_d, double>;
using HEPEnergyType = phys::units::quantity<phys::units::hepenergy_d, double>;
using MassType = phys::units::quantity<phys::units::mass_d, double>;
using HEPMassType = phys::units::quantity<hepmass_d, double>;
using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
using GrammageType = phys::units::quantity<phys::units::dimensions<-2, 1, 0>, double>;
using HEPMomentumType = phys::units::quantity<hepmomentum_d, double>;
using CrossSectionType = phys::units::quantity<area_d, double>;
using InverseLengthType =
phys::units::quantity<phys::units::dimensions<-1, 0, 0>, double>;
using InverseTimeType =
phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
using InverseGrammageType =
phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
namespace detail {
template <int N, typename T>
auto constexpr static_pow([[maybe_unused]] T x) {
if constexpr (N == 0) {
return 1;
} else if constexpr (N > 0) {
return x * static_pow<N - 1, T>(x);
} else {
return 1 / static_pow<-N, T>(x);
}
}
} // namespace detail
template <typename DimFrom, typename DimTo>
auto constexpr ConversionFactorHEPToSI() {
static_assert(DimFrom::dim1 == 0 && DimFrom::dim2 == 0 && DimFrom::dim3 == 0 &&
DimFrom::dim4 == 0 && DimFrom::dim5 == 0 && DimFrom::dim6 == 0 &&
DimFrom::dim7 == 0,
"must be a pure HEP type");
static_assert(
DimTo::dim4 == 0 && DimTo::dim5 == 0 && DimTo::dim6 == 0 && DimTo::dim7 == 0,
"conversion possible only into L, M, T dimensions");
int constexpr e = DimFrom::dim8; // HEP dim.
int constexpr l = DimTo::dim1; // SI length dim.
int constexpr m = DimTo::dim2; // SI mass dim.
int constexpr t = DimTo::dim3; // SI time dim.
int constexpr p = m;
int constexpr q = -m - t;
static_assert(q == l + e - 2 * m, "HEP/SI dimension mismatch!");
using namespace detail;
return static_pow<-e>(corsika::units::constants::hBarC) *
static_pow<p>(corsika::units::constants::hBar) *
static_pow<q>(corsika::units::constants::c);
}
template <typename DimFrom>
auto constexpr ConversionFactorSIToHEP() {
static_assert(DimFrom::dim4 == 0 && DimFrom::dim5 == 0 && DimFrom::dim6 == 0 &&
DimFrom::dim7 == 0 && DimFrom::dim8 == 0,
"must be pure L, M, T type");
int constexpr l = DimFrom::dim1; // SI length dim.
int constexpr m = DimFrom::dim2; // SI mass dim.
int constexpr t = DimFrom::dim3; // SI time dim.
int constexpr p = -m;
int constexpr q = m + t;
int constexpr e = m - t - l;
using namespace detail;
return static_pow<e>(corsika::units::constants::hBarC) *
static_pow<p>(corsika::units::constants::hBar) *
static_pow<q>(corsika::units::constants::c);
}
template <typename DimTo, typename DimFrom>
auto constexpr ConvertHEPToSI(quantity<DimFrom> q) {
return ConversionFactorHEPToSI<DimFrom, DimTo>() * q;
}
template <typename DimFrom>
auto constexpr ConvertSIToHEP(quantity<DimFrom> q) {
return ConversionFactorSIToHEP<DimFrom>() * q;
}
} // end namespace corsika::units::si
/**
* @file PhysicalUnits
*
*/
namespace phys {
namespace units {
namespace literals {
/**
* Define new _XeV literals, alowing 10_GeV in the code.
* Define new _barn literal
*/
QUANTITY_DEFINE_SCALING_LITERALS(eV, hepenergy_d, 1)
QUANTITY_DEFINE_SCALING_LITERALS(b, corsika::units::si::sigma_d,
magnitude(corsika::units::constants::barn))
} // namespace literals
} // namespace units
} // namespace phys
#endif
Framework/Units/testUnits.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
#include <corsika/units/PhysicalUnits.h>
#include <array>
#include <sstream>
using namespace corsika;
using namespace corsika::units::si;
TEST_CASE("PhysicalUnits", "[Units]") {
SECTION("Consistency") {
REQUIRE(1_m / 1_m == Approx(1));
// REQUIRE_FALSE( 1_m/1_s == 1 ); // static assert
}
SECTION("Constructors") {
[[maybe_unused]] auto E1 = 10_GeV;
REQUIRE(E1 == 10_GeV);
LengthType l1 = 10_nm;
[[maybe_unused]] auto l2 = l1;
REQUIRE(l2 == l1);
LengthType arr0[5];
arr0[0] = 5_m;
[[maybe_unused]] LengthType arr1[2] = {{1_mm}, {2_cm}};
[[maybe_unused]] std::array<HEPEnergyType, 4> arr2; // empty array
[[maybe_unused]] std::array<HEPEnergyType, 4> arr3 = {1_GeV, 1_eV, 5_MeV};
auto p1 = 10_s * newton;
REQUIRE(p1 == 10_s * newton);
}
SECTION("Powers in literal units") {
REQUIRE(1_s / 1_ds == Approx(1e1));
REQUIRE(1_m / 1_cm == Approx(1e2));
REQUIRE(1_m / 1_mm == Approx(1e3));
REQUIRE(1_V / 1_uV == Approx(1e6));
REQUIRE(1_s / 1_ns == Approx(1e9));
REQUIRE(1_eV / 1_peV == Approx(1e12));
REQUIRE(1_A / 1_fA == Approx(1e15));
REQUIRE(1_mol / 1_amol == Approx(1e18));
REQUIRE(1_K / 1_zK == Approx(1e21));
REQUIRE(1_K / 1_yK == Approx(1e24));
REQUIRE(1_b / 1_mb == Approx(1e3));
REQUIRE(1_A / 1_hA == Approx(1e-2));
REQUIRE(1_m / 1_km == Approx(1e-3));
REQUIRE(1_m / 1_Mm == Approx(1e-6));
REQUIRE(1_V / 1_GV == Approx(1e-9));
REQUIRE(1_s / 1_Ts == Approx(1e-12));
REQUIRE(1_eV / 1_PeV == Approx(1e-15));
REQUIRE(1_A / 1_EA == Approx(1e-18));
REQUIRE(1_K / 1_ZK == Approx(1e-21));
REQUIRE(1_mol / 1_Ymol == Approx(1e-24));
REQUIRE(std::min(1_A, 2_A) == 1_A);
}
SECTION("Powers and units") {
REQUIRE(1 * ampere / 1_A == Approx(1e0));
REQUIRE(mega * bar / bar == Approx(1e6));
}
SECTION("Formulas") {
const HEPEnergyType E2 = 20_GeV * 2;
REQUIRE(E2 == 40_GeV);
REQUIRE(E2 / 1_GeV == Approx(40));
const MassType m = 1_kg;
const SpeedType v = 1_m / 1_s;
REQUIRE(m * v == 1_s * newton);
const double lgE = log10(E2 / 1_GeV);
REQUIRE(lgE == Approx(log10(40.)));
const auto E3 = E2 + 100_GeV + pow(10, lgE) * 1_GeV;
REQUIRE(E3 == 180_GeV);
}
SECTION("Output") {
{
const HEPEnergyType E = 5_eV;
std::stringstream stream;
stream << E;
REQUIRE(stream.str() == std::string("5 eV"));
}
{
const HEPEnergyType E = 5_EeV;
std::stringstream stream;
stream << E;
REQUIRE(stream.str() == std::string("5e+18 eV"));
}
}
SECTION("Special") {
const LengthType farAway = std::numeric_limits<double>::infinity() * meter;
REQUIRE(farAway > 100000_m);
REQUIRE_FALSE(farAway < 1e19 * meter);
}
SECTION("static_pow") {
using namespace corsika::units::si::detail;
double x = 235.7913;
REQUIRE(1 == static_pow<0, double>(x));
REQUIRE(x == static_pow<1, double>(x));
REQUIRE(x * x == static_pow<2, double>(x));
REQUIRE(1 / x == static_pow<-1, double>(x));
REQUIRE(1 / x / x == static_pow<-2, double>(x));
}
SECTION("HEP/SI conversion") {
auto const invEnergy = 1 / 197.326978_MeV; // should be convertible to length or time
LengthType const length = ConvertHEPToSI<LengthType::dimension_type>(invEnergy);
REQUIRE((length / 1_fm) == Approx(1));
TimeType const time = ConvertHEPToSI<TimeType::dimension_type>(invEnergy);
REQUIRE((time / (1_fm / corsika::units::constants::c)) == Approx(1));
auto const protonMass = 938.272'081'3_MeV; // convertible to mass or SI energy
MassType protonMassSI = ConvertHEPToSI<MassType::dimension_type>(protonMass);
REQUIRE((protonMassSI / 1.672'621'898e-27_kg) == Approx(1));
REQUIRE((protonMassSI / (1.007'276 * corsika::units::constants::u)) == Approx(1));
}
SECTION("SI/HEP conversion") {
REQUIRE(ConvertSIToHEP(units::constants::c) == Approx(1));
REQUIRE(ConvertSIToHEP(units::constants::hBar) == Approx(1));
{
auto const invLength = 1 / 197.326978_fm; // should be convertible to HEPEnergy
HEPEnergyType const energy = ConvertSIToHEP(invLength);
REQUIRE(energy / 1_MeV == Approx(1));
}
REQUIRE(ConvertSIToHEP(6.5823e-25_s) * 1_GeV == Approx(1).epsilon(1e-4));
REQUIRE(ConvertSIToHEP(3.8938e-32 * meter * meter) * 1_GeV * 1_GeV ==
Approx(1).epsilon(1e-4));
}
}

data privacy - impressum