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 2319 deletions
Framework/ProcessSequence/DecayProcess.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/process/BaseProcess.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class DecayProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type DecayProcess<T>
*/
template <typename TDerived>
struct DecayProcess : BaseProcess<TDerived> {
public:
using BaseProcess<TDerived>::GetRef;
/// here starts the interface-definition part
// -> enforce TDerived to implement DoDecay...
template <typename TParticle>
EProcessReturn DoDecay(TParticle&);
template <typename TParticle>
corsika::units::si::TimeType GetLifetime(const TParticle&);
template <typename TParticle>
corsika::units::si::InverseTimeType GetInverseLifetime(const TParticle& particle) {
return 1. / GetRef().GetLifetime(particle);
}
/* template <typename TParticle>
corsika::units::si::InverseTimeType GetInverseInteractionLength(TParticle&& particle)
{ auto p = std::move(particle); return 1. / GetRef().GetLifetime(p);
}*/
};
} // namespace corsika::process
Framework/ProcessSequence/InteractionProcess.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/process/BaseProcess.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class InteractionProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type InteractionProcess<T>
*/
template <typename TDerived>
class InteractionProcess : public BaseProcess<TDerived> {
public:
using BaseProcess<TDerived>::GetRef;
/// here starts the interface-definition part
// -> enforce TDerived to implement DoInteraction...
template <typename TParticle>
EProcessReturn DoInteraction(TParticle&);
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(const TParticle&);
template <typename TParticle>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(
const TParticle& particle) {
return 1. / GetRef().GetInteractionLength(particle);
}
/*
template <typename TParticle>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle&&
particle) { auto p = std::move(particle); return 1. /
GetRef().GetInteractionLength(p);
}*/
};
} // namespace corsika::process
Framework/ProcessSequence/ProcessReturn.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
namespace corsika::process {
/**
since in a process sequence many status updates can accumulate
for a single particle, this enum should define only bit-flags
that can be accumulated easily with "|="
*/
enum class EProcessReturn : int {
eOk = (1 << 0),
eParticleAbsorbed = (1 << 2),
eInteracted = (1 << 3),
eDecayed = (1 << 4),
};
inline EProcessReturn operator|(EProcessReturn a, EProcessReturn b) {
return static_cast<EProcessReturn>(static_cast<int>(a) | static_cast<int>(b));
}
inline EProcessReturn& operator|=(EProcessReturn& a, const EProcessReturn b) {
return a = a | b;
}
inline EProcessReturn operator&(const EProcessReturn a, const EProcessReturn b) {
return static_cast<EProcessReturn>(static_cast<int>(a) & static_cast<int>(b));
}
inline bool operator==(const EProcessReturn a, const EProcessReturn b) {
return (static_cast<int>(a) & static_cast<int>(b)) != 0;
}
inline bool isOk(const EProcessReturn a) {
return static_cast<int>(a & EProcessReturn::eOk);
}
inline bool isAbsorbed(const EProcessReturn a) {
return static_cast<int>(a & EProcessReturn::eParticleAbsorbed);
}
inline bool isDecayed(const EProcessReturn a) {
return static_cast<int>(a & EProcessReturn::eDecayed);
}
inline bool isInteracted(const EProcessReturn a) {
return static_cast<int>(a & EProcessReturn::eInteracted);
}
} // namespace corsika::process
Framework/ProcessSequence/ProcessSequence.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/process/BaseProcess.h>
#include <corsika/process/ProcessTraits.h>
#include <corsika/process/BoundaryCrossingProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/process/StackProcess.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/logging/Logging.h>
#include <cmath>
#include <limits>
namespace corsika::process {
/**
\class ProcessSequence
A compile time static list of processes. The compiler will
generate a new type based on template logic containing all the
elements.
TProcess1 and TProcess2 must both be derived from BaseProcess,
and are both references if possible (lvalue), otherwise (rvalue)
they are just classes. This allows us to handle both, rvalue as
well as lvalue Processes in the ProcessSequence.
\comment Using CRTP pattern,
https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
*/
template <typename TProcess1, typename TProcess2>
class ProcessSequence : public BaseProcess<ProcessSequence<TProcess1, TProcess2>> {
using TProcess1type = typename std::decay<TProcess1>::type;
using TProcess2type = typename std::decay<TProcess2>::type;
static bool constexpr t1ProcSeq = is_process_sequence_v<TProcess1type>;
static bool constexpr t2ProcSeq = is_process_sequence_v<TProcess2type>;
static bool constexpr t1SwitchProcSeq = is_switch_process_sequence_v<TProcess1type>;
static bool constexpr t2SwitchProcSeq = is_switch_process_sequence_v<TProcess2type>;
protected:
TProcess1 A; // this is a reference, if possible
TProcess2 B; // this is a reference, if possible
public:
ProcessSequence(TProcess1 in_A, TProcess2 in_B)
: A(in_A)
, B(in_B) {}
template <typename Particle, typename VTNType>
EProcessReturn DoBoundaryCrossing(Particle& particle, VTNType const& from,
VTNType const& to) {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
ret |= A.DoBoundaryCrossing(particle, from, to);
}
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
ret |= B.DoBoundaryCrossing(particle, from, to);
}
return ret;
}
template <typename TParticle, typename TTrack>
inline EProcessReturn DoContinuous(TParticle& particle, TTrack& vT) {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>, TProcess1type> ||
t1ProcSeq) {
ret |= A.DoContinuous(particle, vT);
}
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>, TProcess2type> ||
t2ProcSeq) {
if (!isAbsorbed(ret)) { ret |= B.DoContinuous(particle, vT); }
}
return ret;
}
template <typename TSecondaries>
inline void DoSecondaries(TSecondaries& vS) {
if constexpr (std::is_base_of_v<SecondariesProcess<TProcess1type>, TProcess1type> ||
t1ProcSeq) {
A.DoSecondaries(vS);
}
if constexpr (std::is_base_of_v<SecondariesProcess<TProcess2type>, TProcess2type> ||
t2ProcSeq) {
B.DoSecondaries(vS);
}
}
/**
The processes of type StackProcess do have an internal counter,
so they can be exectuted only each N steps. Often these are
"maintenacne processes" that do not need to run after each
single step of the simulations. In the CheckStep function it is
tested if either A or B are StackProcess and if they are due
for execution.
*/
inline bool CheckStep() {
bool ret = false;
if constexpr (std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type> ||
(t1ProcSeq && !t1SwitchProcSeq)) {
ret |= A.CheckStep();
}
if constexpr (std::is_base_of_v<StackProcess<TProcess2type>, TProcess2type> ||
(t2ProcSeq && !t2SwitchProcSeq)) {
ret |= B.CheckStep();
}
return ret;
}
/**
Execute the StackProcess-es in the ProcessSequence
*/
template <typename TStack>
inline void DoStack(TStack& stack) {
if constexpr (std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type> ||
(t1ProcSeq && !t1SwitchProcSeq)) {
if (A.CheckStep()) { A.DoStack(stack); }
}
if constexpr (std::is_base_of_v<StackProcess<TProcess2type>, TProcess2type> ||
(t2ProcSeq && !t2SwitchProcSeq)) {
if (B.CheckStep()) { B.DoStack(stack); }
}
}
template <typename TParticle, typename TTrack>
inline corsika::units::si::LengthType MaxStepLength(TParticle& particle,
TTrack& vTrack) {
corsika::units::si::LengthType
max_length = // if no other process in the sequence implements it
std::numeric_limits<double>::infinity() * corsika::units::si::meter;
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>, TProcess1type> ||
t1ProcSeq) {
corsika::units::si::LengthType const len = A.MaxStepLength(particle, vTrack);
max_length = std::min(max_length, len);
}
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>, TProcess2type> ||
t2ProcSeq) {
corsika::units::si::LengthType const len = B.MaxStepLength(particle, vTrack);
max_length = std::min(max_length, len);
}
return max_length;
}
template <typename TParticle>
inline corsika::units::si::GrammageType GetInteractionLength(TParticle&& particle) {
return 1. / GetInverseInteractionLength(particle);
}
template <typename TParticle>
inline corsika::units::si::InverseGrammageType GetInverseInteractionLength(
TParticle&& particle) {
using namespace corsika::units::si;
InverseGrammageType tot = 0 * meter * meter / gram; // default value
if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>, TProcess1type> ||
t1ProcSeq) {
tot += A.GetInverseInteractionLength(particle);
}
if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>, TProcess2type> ||
t2ProcSeq) {
tot += B.GetInverseInteractionLength(particle);
}
return tot;
}
template <typename TSecondaryView>
inline EProcessReturn SelectInteraction(
TSecondaryView& view,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_select,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_sum =
corsika::units::si::InverseGrammageType::zero()) {
// TODO: add check for lambda_inv_select>lambda_inv_tot
if constexpr (t1ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
A.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
// if A did succeed, stop routine. Not checking other static branch B.
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>,
TProcess1type>) {
// if this is not a ContinuousProcess --> evaluate probability
const auto particle = view.parent();
lambda_inv_sum += A.GetInverseInteractionLength(particle);
// check if we should execute THIS process and then EXIT
if (lambda_inv_select < lambda_inv_sum) {
A.DoInteraction(view);
return EProcessReturn::eInteracted;
}
} // end branch A
if constexpr (t2ProcSeq) {
// if B is a process sequence --> check inside
return B.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
} else if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>,
TProcess2type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B.GetInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
if (lambda_inv_select < lambda_inv_sum) {
B.DoInteraction(view);
return EProcessReturn::eInteracted;
}
} // end branch B
return EProcessReturn::eOk;
}
template <typename TParticle>
inline corsika::units::si::TimeType GetLifetime(TParticle& particle) {
return 1. / GetInverseLifetime(particle);
}
template <typename TParticle>
inline corsika::units::si::InverseTimeType GetInverseLifetime(TParticle&& particle) {
using namespace corsika::units::si;
corsika::units::si::InverseTimeType tot = 0 / second; // default value
if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, TProcess1type> ||
t1ProcSeq) {
tot += A.GetInverseLifetime(particle);
}
if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, TProcess2type> ||
t2ProcSeq) {
tot += B.GetInverseLifetime(particle);
}
return tot;
}
// select decay process
template <typename TSecondaryView>
inline EProcessReturn SelectDecay(
TSecondaryView& view,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_select,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_sum =
corsika::units::si::InverseTimeType::zero()) {
// TODO: add check for decay_inv_select>decay_inv_tot
if constexpr (t1ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret = A.SelectDecay(view, decay_inv_select, decay_inv_sum);
// if A did succeed, stop routine here (not checking other static branch B)
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>,
TProcess1type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select <
// decay_inv_sum / decay_inv_tot
A.DoDecay(view);
return EProcessReturn::eDecayed;
}
} // end branch A
if constexpr (t2ProcSeq) {
// if B is a process sequence --> check inside
return B.SelectDecay(view, decay_inv_select, decay_inv_sum);
} else if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>,
TProcess2type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
if (decay_inv_select < decay_inv_sum) {
B.DoDecay(view);
return EProcessReturn::eDecayed;
}
} // end branch B
return EProcessReturn::eOk;
}
};
// the % operator assembles many BaseProcess, ContinuousProcess, and
// Interaction/DecayProcess objects into a ProcessSequence, all combinatorics
// must be allowed, this is why we define a macro to define all
// combinations here:
// enable the % operator to construct ProcessSequence from two
// Processes, only if both, Processes1 and Processes2, derive from
// BaseProcesses
template <typename TProcess1, typename TProcess2>
inline typename std::enable_if<
std::is_base_of<BaseProcess<typename std::decay<TProcess1>::type>,
typename std::decay<TProcess1>::type>::value &&
std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>,
typename std::decay<TProcess2>::type>::value,
ProcessSequence<TProcess1, TProcess2>>::type
operator%(TProcess1&& vA, TProcess2&& vB) {
return ProcessSequence<TProcess1, TProcess2>(vA, vB);
}
/// traits marker to identify objectas ProcessSequence
template <typename TProcess1, typename TProcess2>
struct is_process_sequence<corsika::process::ProcessSequence<TProcess1, TProcess2>>
: std::true_type {};
} // namespace corsika::process
Framework/ProcessSequence/ProcessSignature.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#define FORCE_SIGNATURE(nameTrait, nameMethod, signatureMethod) \
template <typename U> \
class nameTrait { \
private: \
template <typename T, T> \
struct helper; \
template <typename T> \
static std::uint8_t check(helper<signatureMethod, &nameMethod>*); \
template <typename T> \
static std::uint16_t check(...); \
\
public: \
static constexpr bool value = sizeof(check<U>(0)) == sizeof(std::uint8_t); \
}
// FORCE_SIGNATURE(thisMustBeDefined, T::thisMustBeDefined, int(*)(void));
Framework/ProcessSequence/ProcessTraits.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
namespace corsika::process {
/**
* A traits marker to track which BaseProcess is also a ProcessSequence
**/
template <typename TClass>
struct is_process_sequence : std::false_type {};
template <typename TClass>
bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value;
/**
* A traits marker to identiy a BaseProcess that is also SwitchProcessesSequence
**/
template <typename TClass>
struct is_switch_process_sequence : std::false_type {};
template <typename TClass>
bool constexpr is_switch_process_sequence_v = is_switch_process_sequence<TClass>::value;
} // namespace corsika::process
Framework/ProcessSequence/SecondariesProcess.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/process/BaseProcess.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class SecondariesProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type SecondariesProcess<T>
*/
template <typename TDerived>
class SecondariesProcess : public BaseProcess<TDerived> {
public:
/// here starts the interface-definition part
// -> enforce TDerived to implement DoSecondaries...
template <typename TSecondaries>
inline void DoSecondaries(TSecondaries&);
};
} // namespace corsika::process
Framework/ProcessSequence/StackProcess.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/process/BaseProcess.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
/**
\class StackProcess
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type StackProcess<T>
*/
template <typename TDerived>
class StackProcess : public BaseProcess<TDerived> {
private:
protected:
public:
StackProcess() = delete;
StackProcess(const unsigned int nStep)
: fNStep(nStep) {}
/// here starts the interface-definition part
// -> enforce TDerived to implement DoStack...
template <typename TStack>
inline void DoStack(TStack&);
int GetStep() const { return fIStep; }
bool CheckStep() { return !((++fIStep) % fNStep); }
private:
/**
@name The number of "steps" during the cascade processing after
which this StackProcess is going to be executed. The logic is
"fIStep modulo fNStep"
@{
*/
unsigned int fNStep = 0;
unsigned long int fIStep = 0;
//! @}
};
} // namespace corsika::process
Framework/ProcessSequence/SwitchProcessSequence.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/process/BaseProcess.h>
#include <corsika/process/ProcessTraits.h>
#include <corsika/process/BoundaryCrossingProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/process/StackProcess.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
#include <limits>
namespace corsika::process {
// enum for the process switch selection: identify if First or
// Second process branch should be used.
enum class SwitchResult { First, Second };
/**
\class SwitchProcessSequence
A compile time static list of processes that uses an internal
TSelect class to switch between different versions of processes
(or process sequence).
TProcess1 and TProcess2 must be derived from BaseProcess and are
both references if possible (lvalue), otherwise (rvalue) they are
just classes. This allows us to handle both, rvalue as well as
lvalue Processes in the SwitchProcessSequence.
TSelect has to implement a `select(const Particle&)` and has to
return either SwitchResult::First or SwitchResult::Second. Note:
TSelect may absolutely also use random numbers to sample between
its results. This can be used to achieve arbitrarily smooth
transition or mixtures of processes.
Warning: do not put StackProcess into a SwitchProcessSequence
since this makes no sense. The StackProcess acts on an entire
particle stack and not on indiviidual particles.
\comment See also class ProcessSequence
**/
template <typename TProcess1, typename TProcess2, typename TSelect>
class SwitchProcessSequence
: public BaseProcess<SwitchProcessSequence<TProcess1, TProcess2, TSelect>> {
using TProcess1type = typename std::decay<TProcess1>::type;
using TProcess2type = typename std::decay<TProcess2>::type;
static bool constexpr t1ProcSeq = is_process_sequence_v<TProcess1type>;
static bool constexpr t2ProcSeq = is_process_sequence_v<TProcess2type>;
TSelect select_; // this is a reference, if possible
public:
TProcess1 A; // this is a reference, if possible
TProcess2 B; // this is a reference, if possible
SwitchProcessSequence(TProcess1 in_A, TProcess2 in_B, TSelect sel)
: select_(sel)
, A(in_A)
, B(in_B) {}
template <typename Particle, typename VTNType>
EProcessReturn DoBoundaryCrossing(Particle& particle, VTNType const& from,
VTNType const& to) {
switch (select_.select(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
return A.DoBoundaryCrossing(particle, from, to);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
return B.DoBoundaryCrossing(particle, from, to);
}
break;
}
}
return EProcessReturn::eOk;
}
template <typename TParticle, typename TTrack>
inline EProcessReturn DoContinuous(TParticle& particle, TTrack& vT) {
switch (select_.select(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
return A.DoContinuous(particle, vT);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
return B.DoContinuous(particle, vT);
}
break;
}
}
return EProcessReturn::eOk;
}
template <typename TSecondaries>
inline void DoSecondaries(TSecondaries& vS) {
const auto& particle = vS.parent();
switch (select_.select(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<SecondariesProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
A.DoSecondaries(vS);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<SecondariesProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
B.DoSecondaries(vS);
}
break;
}
}
}
template <typename TParticle, typename TTrack>
inline corsika::units::si::LengthType MaxStepLength(TParticle& particle,
TTrack& vTrack) {
switch (select_.select(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
return A.MaxStepLength(particle, vTrack);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
return B.MaxStepLength(particle, vTrack);
}
break;
}
}
// if no other process in the sequence implements it
return std::numeric_limits<double>::infinity() * corsika::units::si::meter;
}
template <typename TParticle>
inline corsika::units::si::GrammageType GetInteractionLength(TParticle&& particle) {
return 1. / GetInverseInteractionLength(particle);
}
template <typename TParticle>
inline corsika::units::si::InverseGrammageType GetInverseInteractionLength(
TParticle&& particle) {
using namespace corsika::units::si;
switch (select_.select(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
return A.GetInverseInteractionLength(particle);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
return B.GetInverseInteractionLength(particle);
}
break;
}
}
return 0 * meter * meter / gram; // default value
}
template <typename TSecondaryView>
inline EProcessReturn SelectInteraction(
TSecondaryView& view,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_select,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_sum =
corsika::units::si::InverseGrammageType::zero()) {
switch (select_.select(view.parent())) {
case SwitchResult::First: {
if constexpr (t1ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
A.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
// if A did succeed, stop routine. Not checking other static branch B.
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>,
TProcess1type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A.GetInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
if (lambda_inv_select < lambda_inv_sum) {
A.DoInteraction(view);
return EProcessReturn::eInteracted;
}
} // end branch A
break;
}
case SwitchResult::Second: {
if constexpr (t2ProcSeq) {
// if B is a process sequence --> check inside
return B.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
} else if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>,
TProcess2type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B.GetInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
if (lambda_inv_select < lambda_inv_sum) {
B.DoInteraction(view);
return EProcessReturn::eInteracted;
}
} // end branch B
break;
}
}
return EProcessReturn::eOk;
}
template <typename TParticle>
inline corsika::units::si::TimeType GetLifetime(TParticle&& particle) {
return 1. / GetInverseLifetime(particle);
}
template <typename TParticle>
inline corsika::units::si::InverseTimeType GetInverseLifetime(TParticle&& particle) {
using namespace corsika::units::si;
switch (select_.select(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, TProcess1type> ||
t1ProcSeq) {
return A.GetInverseLifetime(particle);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, TProcess2type> ||
t2ProcSeq) {
return B.GetInverseLifetime(particle);
}
break;
}
}
return 0 / second; // default value
}
// select decay process
template <typename TSecondaryView>
inline EProcessReturn SelectDecay(
TSecondaryView& view,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_select,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_sum =
corsika::units::si::InverseTimeType::zero()) {
switch (select_.select(view.parent())) {
case SwitchResult::First: {
if constexpr (t1ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
A.SelectDecay(view, decay_inv_select, decay_inv_sum);
// if A did succeed, stop routine here (not checking other static branch B)
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>,
TProcess1type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
if (decay_inv_select < decay_inv_sum) {
// more pedagogical: rndm_select < decay_inv_sum / decay_inv_tot
A.DoDecay(view);
return EProcessReturn::eDecayed;
}
} // end branch A
break;
}
case SwitchResult::Second: {
if constexpr (t2ProcSeq) {
// if B is a process sequence --> check inside
return B.SelectDecay(view, decay_inv_select, decay_inv_sum);
} else if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>,
TProcess2type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
if (decay_inv_select < decay_inv_sum) {
B.DoDecay(view);
return EProcessReturn::eDecayed;
}
} // end branch B
break;
}
}
return EProcessReturn::eOk;
}
};
// the method `select(proc1,proc1,selector)` assembles many
// BaseProcesses, and ProcessSequences into a SwitchProcessSequence,
// all combinatorics must be allowed, this is why we define a macro
// to define all combinations here:
// Both, Processes1 and Processes2, must derive from BaseProcesses
template <typename TProcess1, typename TProcess2, typename TSelect>
inline typename std::enable_if<
std::is_base_of<BaseProcess<typename std::decay<TProcess1>::type>,
typename std::decay<TProcess1>::type>::value &&
std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>,
typename std::decay<TProcess2>::type>::value,
SwitchProcessSequence<TProcess1, TProcess2, TSelect>>::type
select(TProcess1&& vA, TProcess2&& vB, TSelect selector) {
return SwitchProcessSequence<TProcess1, TProcess2, TSelect>(vA, vB, selector);
}
/// traits marker to identify objectas ProcessSequence
template <typename TProcess1, typename TProcess2, typename TSelect>
struct is_process_sequence<
corsika::process::SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
: std::true_type {};
/// traits marker to identify objectas SwitchProcessSequence
template <typename TProcess1, typename TProcess2, typename TSelect>
struct is_switch_process_sequence<
corsika::process::SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
: std::true_type {};
} // namespace corsika::process
Framework/ProcessSequence/testProcessSequence.cc deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <catch2/catch.hpp>
#include <array>
#include <iomanip>
#include <iostream>
#include <typeinfo>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/SwitchProcessSequence.h>
#include <corsika/units/PhysicalUnits.h>
using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::process;
using namespace std;
static const int nData = 10;
int globalCount = 0; // simple counter
int checkDecay = 0; // use this as a bit field
int checkInteract = 0; // use this as a bit field
int checkSec = 0; // use this as a bit field
int checkCont = 0; // use this as a bit field
class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
int fV = 0;
public:
ContinuousProcess1(const int v)
: fV(v) {
cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
globalCount++;
}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess1::DoContinuous" << endl;
checkCont |= 1;
for (int i = 0; i < nData; ++i) d.data_[i] += 0.933;
return EProcessReturn::eOk;
}
};
class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> {
int fV = 0;
public:
ContinuousProcess2(const int v)
: fV(v) {
cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
globalCount++;
}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess2::DoContinuous" << endl;
checkCont |= 2;
for (int i = 0; i < nData; ++i) d.data_[i] += 0.933;
return EProcessReturn::eOk;
}
};
class ContinuousProcess3 : public ContinuousProcess<ContinuousProcess3> {
int fV = 0;
public:
ContinuousProcess3(const int v)
: fV(v) {
cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
globalCount++;
}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess3::DoContinuous" << endl;
checkCont |= 4;
for (int i = 0; i < nData; ++i) d.data_[i] += 0.933;
return EProcessReturn::eOk;
}
};
class Process1 : public InteractionProcess<Process1> {
public:
Process1(const int v)
: fV(v) {
cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
globalCount++;
}
template <typename TView>
inline EProcessReturn DoInteraction(TView& v) const {
checkInteract |= 1;
for (int i = 0; i < nData; ++i) v.parent().data_[i] += 1 + i;
return EProcessReturn::eOk;
}
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle&) const {
return 10_g / square(1_cm);
}
private:
int fV;
};
class Process2 : public InteractionProcess<Process2> {
int fV = 0;
public:
Process2(const int v)
: fV(v) {
cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
globalCount++;
}
template <typename TView>
inline EProcessReturn DoInteraction(TView&) const {
checkInteract |= 2;
cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle>
GrammageType GetInteractionLength(Particle&) const {
cout << "Process2::GetInteractionLength" << endl;
return 20_g / (1_cm * 1_cm);
}
};
class Process3 : public InteractionProcess<Process3> {
int fV = 0;
public:
Process3(const int v)
: fV(v) {
cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
globalCount++;
}
template <typename TView>
inline EProcessReturn DoInteraction(TView&) const {
checkInteract |= 4;
cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle>
GrammageType GetInteractionLength(Particle&) const {
cout << "Process3::GetInteractionLength" << endl;
return 30_g / (1_cm * 1_cm);
}
};
class Process4 : public BaseProcess<Process4> {
int fV = 0;
public:
Process4(const int v)
: fV(v) {
cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
globalCount++;
}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
std::cout << "Base::DoContinuous" << std::endl;
checkCont |= 8;
for (int i = 0; i < nData; ++i) { d.data_[i] /= 1.2; }
return EProcessReturn::eOk;
}
template <typename TView>
EProcessReturn DoInteraction(TView&) const {
checkInteract |= 8;
return EProcessReturn::eOk;
}
};
class Decay1 : public DecayProcess<Decay1> {
public:
Decay1(const int) {
cout << "Decay1()" << endl;
globalCount++;
}
template <typename Particle>
TimeType GetLifetime(Particle&) const {
return 1_s;
}
template <typename TView>
EProcessReturn DoDecay(TView&) const {
checkDecay |= 1;
return EProcessReturn::eOk;
}
};
class Decay2 : public DecayProcess<Decay2> {
public:
Decay2(const int) {
cout << "Decay2()" << endl;
globalCount++;
}
template <typename Particle>
TimeType GetLifetime(Particle&) const {
return 2_s;
}
template <typename TView>
EProcessReturn DoDecay(TView&) const {
checkDecay |= 2;
return EProcessReturn::eOk;
}
};
class Stack1 : public StackProcess<Stack1> {
int fCount = 0;
public:
Stack1(const int n)
: StackProcess(n) {}
template <typename TStack>
EProcessReturn DoStack(TStack&) {
fCount++;
return EProcessReturn::eOk;
}
int GetCount() const { return fCount; }
};
struct DummyStack {};
struct DummyData {
double data_[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
struct DummyTrajectory {};
struct DummyView {
DummyView(DummyData& p)
: p_(p) {}
DummyData& p_;
DummyData& parent() { return p_; }
};
TEST_CASE("Process Sequence", "[Process Sequence]") {
SECTION("Check construction") {
globalCount = 0;
Process1 m1(0);
CHECK(globalCount == 1);
Process2 m2(1);
CHECK(globalCount == 2);
Process3 m3(2);
CHECK(globalCount == 3);
Process4 m4(3);
CHECK(globalCount == 4);
auto sequence = m1 % m2 % m3 % m4;
CHECK(is_process_sequence_v<decltype(sequence)> == true);
CHECK(is_process_sequence_v<decltype(m2)> == false);
}
SECTION("interaction length") {
globalCount = 0;
ContinuousProcess1 cp1(0);
Process2 m2(1);
Process3 m3(2);
DummyData particle;
auto sequence2 = cp1 % m2 % m3;
GrammageType const tot = sequence2.GetInteractionLength(particle);
InverseGrammageType const tot_inv = sequence2.GetInverseInteractionLength(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
globalCount = 0;
}
SECTION("lifetime") {
globalCount = 0;
ContinuousProcess1 cp1(0);
Process2 m2(1);
Process3 m3(2);
Decay1 d3(3);
DummyData particle;
auto sequence2 = cp1 % m2 % m3 % d3;
TimeType const tot = sequence2.GetLifetime(particle);
InverseTimeType const tot_inv = sequence2.GetInverseLifetime(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
globalCount = 0;
}
SECTION("sectionTwo") {
globalCount = 0;
ContinuousProcess1 cp1(0);
ContinuousProcess2 cp2(1);
Process2 m2(2);
Process3 m3(3);
auto sequence2 = cp1 % m2 % m3 % cp2;
DummyData particle;
DummyTrajectory track;
cout << "-->init sequence2" << endl;
globalCount = 0;
cout << "-->docont" << endl;
sequence2.DoContinuous(particle, track);
cout << "-->dodisc" << endl;
cout << "-->done" << endl;
const int nLoop = 5;
cout << "Running loop with n=" << nLoop << endl;
for (int i = 0; i < nLoop; ++i) { sequence2.DoContinuous(particle, track); }
for (int i = 0; i < nData; i++) {
cout << "data_[" << i << "]=" << particle.data_[i] << endl;
}
cout << "done" << endl;
}
SECTION("StackProcess") {
globalCount = 0;
Stack1 s1(1);
Stack1 s2(2);
auto sequence = s1 % s2;
DummyStack stack;
const int nLoop = 20;
for (int i = 0; i < nLoop; ++i) { sequence.DoStack(stack); }
CHECK(s1.GetCount() == 20);
CHECK(s2.GetCount() == 10);
}
}
TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
SECTION("Check construction") {
struct TestSelect {
corsika::process::SwitchResult select(const DummyData& p) const {
std::cout << "TestSelect data=" << p.data_[0] << std::endl;
if (p.data_[0] > 0) return corsika::process::SwitchResult::First;
return corsika::process::SwitchResult::Second;
}
};
TestSelect select;
auto sequence1 = Process1(0) % ContinuousProcess2(0) % Decay1(0);
auto sequence2 = ContinuousProcess3(0) % Process2(0) % Decay2(0);
auto sequence = ContinuousProcess1(0) % Process3(0) %
SwitchProcessSequence(sequence1, sequence2, select);
auto sequence_alt =
(ContinuousProcess1(0) % Process3(0)) %
process::select(Process1(0) % ContinuousProcess2(0) % Decay1(0),
ContinuousProcess3(0) % Process2(0) % Decay2(0), select);
// check that same process sequence can be build in different ways
CHECK(typeid(sequence) == typeid(sequence_alt));
CHECK(is_process_sequence_v<decltype(sequence)> == true);
CHECK(is_process_sequence_v<decltype(
SwitchProcessSequence(sequence1, sequence2, select))> == true);
DummyData particle;
DummyTrajectory track;
DummyView view(particle);
checkDecay = 0;
checkInteract = 0;
checkSec = 0;
checkCont = 0;
particle.data_[0] = 100; // data positive
sequence.DoContinuous(particle, track);
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
CHECK(checkCont == 0b011);
CHECK(checkSec == 0);
checkDecay = 0;
checkInteract = 0;
checkSec = 0;
checkCont = 0;
particle.data_[0] = -100; // data negative
sequence_alt.DoContinuous(particle, track);
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
CHECK(checkCont == 0b101);
CHECK(checkSec == 0);
// 1/(30g/cm2) is Process3
corsika::units::si::InverseGrammageType lambda_select = .9 / 30. * square(1_cm) / 1_g;
corsika::units::si::InverseTimeType time_select = 0.1 / second;
checkDecay = 0;
checkInteract = 0;
checkSec = 0;
checkCont = 0;
particle.data_[0] = 100; // data positive
sequence.SelectInteraction(view, lambda_select);
sequence.SelectDecay(view, time_select);
CHECK(checkInteract == 0b100); // this is Process3
CHECK(checkDecay == 0b001); // this is Decay1
CHECK(checkCont == 0);
CHECK(checkSec == 0);
lambda_select = 1.01 / 30. * square(1_cm) / 1_g;
checkInteract = 0;
sequence.SelectInteraction(view, lambda_select);
CHECK(checkInteract == 0b001); // this is Process1
checkDecay = 0;
checkInteract = 0;
checkSec = 0;
checkCont = 0;
particle.data_[0] = -100; // data negative
sequence.SelectInteraction(view, lambda_select);
sequence.SelectDecay(view, time_select);
CHECK(checkInteract == 0b010); // this is Process2
CHECK(checkDecay == 0b010); // this is Decay2
CHECK(checkCont == 0);
CHECK(checkSec == 0);
checkDecay = 0;
checkInteract = 0;
checkSec = 0;
checkCont = 0;
particle.data_[0] = -100; // data negative
sequence.DoSecondaries(view);
Stack1 stack(0);
sequence.DoStack(stack);
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
CHECK(checkCont == 0);
CHECK(checkSec == 0);
}
}
Framework/ProcessSequence/testSwitchProcessSequence.cc deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
#include <catch2/catch.hpp>
#include <algorithm>
#include <random>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units::si;
class TestStackData {
public:
// these functions are needed for the Stack interface
void Clear() { fData.clear(); }
unsigned int GetSize() const { return fData.size(); }
unsigned int GetCapacity() const { return fData.size(); }
void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
// custom data access function
void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fData.resize(fData.size() + 1); }
void DecrementSize() {
if (fData.size() > 0) { fData.pop_back(); }
}
// custom private data section
private:
std::vector<HEPEnergyType> 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<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
std::tuple<HEPEnergyType> v) {
SetEnergy(std::get<0>(v));
}
// here are the fundamental methods for access to TestStackData data
void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
};
using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
// see issue 161
#if defined(__clang__)
using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
#endif
auto constexpr kgMSq = 1_kg / (1_m * 1_m);
template <int N>
struct DummyProcess : InteractionProcess<DummyProcess<N>> {
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
return N * kgMSq;
}
template <typename TSecondaries>
corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
// to figure out which process was selected in the end, we produce N
// secondaries for DummyProcess<N>
for (int i = 0; i < N; ++i) {
vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N});
}
return EProcessReturn::eOk;
}
};
using DummyLowEnergyProcess = DummyProcess<1>;
using DummyHighEnergyProcess = DummyProcess<2>;
using DummyAdditionalProcess = DummyProcess<3>;
TEST_CASE("SwitchProcess from InteractionProcess") {
DummyLowEnergyProcess low;
DummyHighEnergyProcess high;
DummyAdditionalProcess proc;
switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
auto seq = switchProcess << proc;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
// low energy process returns 1 kg/m²
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
REQUIRE(seq.GetInteractionLength(p) / kgMSq == Approx(3. / 4));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
auto p = stack.GetNextParticle();
// high energy process returns 2 kg/m²
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
REQUIRE(seq.GetInteractionLength(p) / kgMSq == Approx(6. / 5));
}
// high energy process creates 2 secondaries
SECTION("SelectInteraction") {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
InverseGrammageType invLambda = 0 / kgMSq;
switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
REQUIRE(view.getSize() == 2);
}
}
}
TEST_CASE("SwitchProcess from ProcessSequence") {
DummyProcess<1> innerA;
DummyProcess<2> innerB;
DummyProcess<3> outer;
DummyProcess<4> additional;
auto seq = innerA << innerB;
switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
auto completeSeq = switchProcess << additional;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
REQUIRE(completeSeq.GetInteractionLength(p) / kgMSq == Approx(4. / 7));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 4 / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.getSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
REQUIRE(mean == Approx(12. / 7.).margin(0.01));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
REQUIRE(completeSeq.GetInteractionLength(p) / kgMSq == Approx(12. / 7.));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.getSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
REQUIRE(mean == Approx(24. / 7.).margin(0.01));
}
}
}
Framework/ProcessSequence/testSwitchProcessSequence.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
#include <catch2/catch.hpp>
#include <algorithm>
#include <random>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units::si;
class TestStackData {
public:
// these functions are needed for the Stack interface
void Clear() { fData.clear(); }
unsigned int GetSize() const { return fData.size(); }
unsigned int GetCapacity() const { return fData.size(); }
void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
// custom data access function
void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fData.resize(fData.size() + 1); }
void DecrementSize() {
if (fData.size() > 0) { fData.pop_back(); }
}
// custom private data section
private:
std::vector<HEPEnergyType> 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<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
std::tuple<HEPEnergyType> v) {
SetEnergy(std::get<0>(v));
}
// here are the fundamental methods for access to TestStackData data
void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
};
using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
// see issue 161
#if defined(__clang__)
using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
#endif
auto constexpr kgMSq = 1_kg / (1_m * 1_m);
template <int N>
struct DummyProcess : InteractionProcess<DummyProcess<N>> {
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
return N * kgMSq;
}
template <typename TSecondaries>
corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
// to figure out which process was selected in the end, we produce N
// secondaries for DummyProcess<N>
for (int i = 0; i < N; ++i) {
vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N});
}
return EProcessReturn::eOk;
}
};
using DummyLowEnergyProcess = DummyProcess<1>;
using DummyHighEnergyProcess = DummyProcess<2>;
using DummyAdditionalProcess = DummyProcess<3>;
TEST_CASE("SwitchProcess from InteractionProcess") {
DummyLowEnergyProcess low;
DummyHighEnergyProcess high;
DummyAdditionalProcess proc;
switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
auto seq = switchProcess << proc;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
// low energy process returns 1 kg/m²
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(3. / 4));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
auto p = stack.GetNextParticle();
// high energy process returns 2 kg/m²
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(6. / 5));
}
// high energy process creates 2 secondaries
SECTION("SelectInteraction") {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
InverseGrammageType invLambda = 0 / kgMSq;
switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
REQUIRE(view.getSize() == 2);
}
}
}
TEST_CASE("SwitchProcess from ProcessSequence") {
DummyProcess<1> innerA;
DummyProcess<2> innerB;
DummyProcess<3> outer;
DummyProcess<4> additional;
auto seq = innerA << innerB;
switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
auto completeSeq = switchProcess << additional;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(4. / 7));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 4 / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.getSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
REQUIRE(mean == Approx(12. / 7.).margin(0.01));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(12. / 7.));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.getSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
REQUIRE(mean == Approx(24. / 7.).margin(0.01));
}
}
}
Framework/Random/CMakeLists.txt deleted 100644 → 0
View file @ 17d0edcc
set (
CORSIKArandom_SOURCES
RNGManager.cc
)
set (
CORSIKArandom_HEADERS
RNGManager.h
UniformRealDistribution.h
ExponentialDistribution.h
)
set (
CORSIKArandom_NAMESPACE
corsika/random
)
add_library (CORSIKArandom STATIC ${CORSIKArandom_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKArandom ${CORSIKArandom_NAMESPACE} ${CORSIKArandom_HEADERS})
target_link_libraries (
CORSIKArandom
CORSIKAutilities
CORSIKAlogging
)
target_include_directories (
CORSIKArandom
PUBLIC
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/>
)
# target dependencies on other libraries (also the header onlys)
# none
install (
TARGETS CORSIKArandom
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/${CORSIKArandom_NAMESPACE}
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testRandom)
target_link_libraries (
testRandom
CORSIKArandom
CORSIKAtesting
CORSIKAunits
)
Framework/Random/ExponentialDistribution.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/units/PhysicalUnits.h>
#include <random>
namespace corsika::random {
template <class TQuantity>
class ExponentialDistribution {
using RealType = typename TQuantity::value_type;
std::exponential_distribution<RealType> dist{1.};
TQuantity const fBeta;
public:
ExponentialDistribution(TQuantity beta)
: fBeta(beta) {}
template <class Generator>
TQuantity operator()(Generator& g) {
return fBeta * dist(g);
}
};
} // namespace corsika::random
Framework/Random/RNGManager.cc deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/random/RNGManager.h>
#include <corsika/logging/Logging.h>
#include <sstream>
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) {
if (IsRegistered(pStreamName)) {
return rngs.at(pStreamName);
} else { // this stream name is not in the map
throw std::runtime_error("'" + pStreamName + "' is not a registered stream.");
}
}
bool corsika::random::RNGManager::IsRegistered(std::string const& pStreamName) const {
return rngs.count(pStreamName) > 0;
}
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::SeedAll(uint64_t vSeed) {
for (auto& entry : rngs) {
auto seed = vSeed++;
C8LOG_TRACE("Random seed stream {} seed {}", entry.first, seed);
entry.second.seed(seed);
}
}
void corsika::random::RNGManager::SeedAll() {
std::random_device rd;
std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()};
for (auto& entry : rngs) {
std::vector<std::uint32_t> seeds(1);
sseq.generate(seeds.begin(), seeds.end());
std::uint32_t seed = seeds[0];
C8LOG_TRACE("Random seed stream {} seed {}", entry.first, seed);
entry.second.seed(seed);
}
}
Framework/Random/RNGManager.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/utl/Singleton.h>
#include <map>
#include <random>
#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
/*!
* Manage random number generators.
*/
class RNGManager final : 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);
/*!
* Check whether a stream has been registered.
*/
bool IsRegistered(std::string const& pStreamName) const;
/*!
* dumps the names and states of all registered random-number streams
* into a std::stringstream.
*/
std::stringstream dumpState() const;
/**
* Set explicit seeds for all currently registered streams. The actual seed values
* are incremented from \a vSeed.
*/
void SeedAll(uint64_t vSeed);
void SeedAll(); //!< seed all currently registered streams with "real" randomness
};
} // namespace corsika::random
Framework/Random/UniformRealDistribution.h deleted 100644 → 0
View file @ 17d0edcc
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <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
Framework/StackInterface/CMakeLists.txt deleted 100644 → 0
View file @ 17d0edcc
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_link_libraries (
CORSIKAstackinterface
INTERFACE
CORSIKAlogging
CORSIKAsetup
)
target_include_directories (
CORSIKAstackinterface
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
FILES ${CORSIKAstackinterface_HEADERS}
DESTINATION include/${CORSIKAstackinterface_NAMESPACE}
)
#code testing
CORSIKA_ADD_TEST(testStackInterface)
target_link_libraries (testStackInterface CORSIKAstackinterface CORSIKAtesting)
CORSIKA_ADD_TEST(testSecondaryView)
target_link_libraries (testSecondaryView CORSIKAstackinterface CORSIKAtesting)
CORSIKA_ADD_TEST(testCombinedStack)
target_link_libraries (testCombinedStack CORSIKAstackinterface CORSIKAtesting)
Framework/StackInterface/Stack.dox deleted 100644 → 0
View file @ 17d0edcc
/**
@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/Testing/CMakeLists.txt deleted 100644 → 0
View file @ 17d0edcc
set (
TESTING_SOURCES
TestMain.cc
)
set (
TESTING_HEADERS
)
set (
TESTING_NAMESPACE
corsika/testing
)
add_library (CORSIKAtesting STATIC ${TESTING_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAtesting ${TESTING_NAMESPACE} ${TESTING_HEADERS})
set_target_properties (
CORSIKAtesting
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
PUBLIC_HEADER "${TESTING_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAtesting
CORSIKAthirdparty
)
# target_include_directories (
# CORSIKAtesting
# SYSTEM
# PUBLIC ${CATCH2_INCLUDE_DIR}
# INTERFACE ${CATCH2_INCLUDE_DIR}
# )
target_include_directories (
CORSIKAtesting
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS CORSIKAtesting
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/${TESTING_NAMESPACE}
)

data privacy - impressum