diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 51d7099500a4606ad818e91c6437c6872a6b9e8c..cad21e7420e61a9fba3225d62d93e8b1a122b857 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -71,8 +71,9 @@ namespace corsika::cascade {
fProcessSequence.GetTotalInverseInteractionLength(particle, step);
// sample random exponential step length in grammage
- std::exponential_distribution expDist((1_m * 1_m / 1_g) / total_inv_lambda);
- GrammageType const next_interact = (1_g / (1_m * 1_m)) * expDist(fRNG);
+ auto constexpr grammageConversion = 1_g / (1_m * 1_m);
+ std::exponential_distribution expDist(1 / (grammageConversion * total_inv_lambda));
+ GrammageType const next_interact = grammageConversion * expDist(fRNG);
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index 86779739c3e7f2ce480b95501fb29907a47659ce..19effade79572ef87d0091f9e86f416dbfdac30e 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/Framework/ProcessSequence/ContinuousProcess.h
@@ -13,6 +13,7 @@
#define _include_corsika_continuousprocess_h_
#include <corsika/process/ProcessReturn.h> // for convenience
+#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -32,8 +33,12 @@ namespace corsika::process {
// here starts the interface part
// -> enforce derived to implement DoContinuous...
- template <typename P, typename T, typename S>
- inline EProcessReturn DoContinuous(P&, T&, S&) const;
+ template <typename Particle, typename Track, typename Stack>
+ EProcessReturn DoContinuous(Particle&, Track&, Stack&) const;
+
+ // -> enforce derived to implement MaxStepLength...
+ template <typename Particle, typename Track>
+ corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const;
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h
index 076a32410495f8e7524b265284f7aa1d4f8417d5..0175cee60f8a4b8f8e0e4fd033cf07d64c239209 100644
--- a/Framework/ProcessSequence/DecayProcess.h
+++ b/Framework/ProcessSequence/DecayProcess.h
@@ -29,14 +29,13 @@ namespace corsika::process {
template <typename derived>
struct DecayProcess {
-
derived& GetRef() { return static_cast<derived&>(*this); }
const derived& GetRef() const { return static_cast<const derived&>(*this); }
/// here starts the interface-definition part
// -> enforce derived to implement DoDecay...
template <typename Particle, typename Stack>
- inline EProcessReturn DoDecay(Particle&, Stack&) const;
+ EProcessReturn DoDecay(Particle&, Stack&) const;
template <typename Particle>
corsika::units::si::TimeType GetLifetime(Particle& p) const;
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index bc4bf679e4c92bf5d5dc60188ecb3404eef05a35..b2c4fed30cbc13e69ad71633784e332d0a2af93e 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -71,8 +71,9 @@ namespace corsika::process {
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle& p, Track& track) const {
- corsika::units::si::LengthType max_length =
+ 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<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
corsika::units::si::LengthType const len = A.MaxStepLength(p, track);
diff --git a/Processes/NullModel/NullModel.cc b/Processes/NullModel/NullModel.cc
index 313d38034d3c90925eb461aa59b1399af8308e4d..4130bf0c4d5519408b272322dd45d60ee3b27b0e 100644
--- a/Processes/NullModel/NullModel.cc
+++ b/Processes/NullModel/NullModel.cc
@@ -1,4 +1,3 @@
-
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -11,37 +10,4 @@
#include <corsika/process/null_model/NullModel.h>
-#include <corsika/logging/Logger.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <iostream>
-#include <limits>
-using namespace std;
-
-using namespace corsika;
-using namespace corsika::units::si;
-using namespace corsika::process::null_model;
-
-template <typename Stack>
-NullModel<Stack>::NullModel() {}
-
-template <typename Stack>
-NullModel<Stack>::~NullModel() {}
-
-template <typename Stack>
-process::EProcessReturn NullModel<Stack>::DoContinuous(Particle&, setup::Trajectory&,
- Stack&) const {
- return EProcessReturn::eOk;
-}
-
-template <typename Stack>
-double NullModel<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const {
- return std::numeric_limits<double>::infinity();
-}
-
-template <typename Stack>
-void NullModel<Stack>::Init() {}
-
-#include <corsika/setup/SetupStack.h>
-
-template class process::null_model::NullModel<setup::Stack>;
+void corsika::process::null_model::NullModel::Init() {}
diff --git a/Processes/NullModel/NullModel.h b/Processes/NullModel/NullModel.h
index 367c907f38df5ab6b04bbc23ccaefc24a54df4ad..e1fc5b3ff8eaeeffdc5b6703f9dcc9f51516d87b 100644
--- a/Processes/NullModel/NullModel.h
+++ b/Processes/NullModel/NullModel.h
@@ -1,4 +1,3 @@
-
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -15,26 +14,29 @@
#include <corsika/process/ContinuousProcess.h>
#include <corsika/setup/SetupTrajectory.h>
-namespace corsika::process {
-
- namespace null_model {
+namespace corsika::process::null_model {
- template <typename Stack>
- class NullModel {
+ class NullModel : public corsika::process::ContinuousProcess<NullModel> {
+ corsika::units::si::LengthType fMaxStepLength{
+ corsika::units::si::meter * std::numeric_limits<double>::infinity()};
- typedef typename Stack::ParticleType Particle;
+ public:
+ NullModel(corsika::units::si::LengthType maxStepLength)
+ : fMaxStepLength(maxStepLength) {}
- public:
- NullModel();
- ~NullModel();
+ void Init();
- void Init();
- EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s) const;
- double MaxStepLength(Particle&, corsika::setup::Trajectory&) const;
- };
+ template <typename Particle, typename Track, typename Stack>
+ process::EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
- } // namespace null_model
+ template <typename Particle, typename Track>
+ corsika::units::si::LengthType MaxStepLength(Particle&, Track&) const {
+ return fMaxStepLength;
+ }
+ };
-} // namespace corsika::process
+} // namespace corsika::process::null_model
#endif
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
index ffcba5e0351acefa519f6a28e003be7999741f68..56d635d1fafbf1806ed0a940985040934aa6a31d 100644
--- a/Processes/NullModel/testNullModel.cc
+++ b/Processes/NullModel/testNullModel.cc
@@ -43,11 +43,13 @@ TEST_CASE("NullModel", "[processes]") {
SECTION("interface") {
- NullModel<setup::Stack> model;
+ NullModel model(10_m);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoContinuous(particle, track, stack);
- [[maybe_unused]] const double length = model.MaxStepLength(particle, track);
+ LengthType const length = model.MaxStepLength(particle, track);
+
+ CHECK((length / 10_m) == Approx(1));
}
}
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 24fd807d2723343be6a827f0783fa8fd5fdd0383..a8c586061bb7e119588bab7ba8fb3b1f67ef7bec 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -86,9 +86,8 @@ namespace corsika::process {
// i.e. corsika::particles::ListOfParticles()
std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
// make ALL particles unstable, then set EM stable
- for (auto& p : corsika2sibyll) {
+ for (int sibCode : corsika2sibyll) {
// std::cout << (int)p << std::endl;
- const int sibCode = (int)p;
// skip unknown and antiparticles
if (sibCode < 1) continue;
// std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
@@ -128,19 +127,14 @@ namespace corsika::process {
const corsika::units::si::TimeType t0 =
corsika::particles::GetLifetime(p.GetPID());
+ auto const lifetime = gamma * t0;
+
cout << "Decay: code: " << p.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl;
- cout << "Decay: MinStep: energy: " << E / 1_GeV << endl;
+ cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl;
- // cout << "Decay: MinStep: density: " << density << endl;
- // return as column density
- // const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
- // cout << "Decay: MinStep: x0: " << x0 << endl;
- corsika::units::si::TimeType const lifetime = gamma * t0;
cout << "Decay: MinStep: tau: " << lifetime << endl;
- // int a = 1;
- // const double x = -x0 * log(s_rndm_(a));
- // cout << "Decay: next decay: " << x << endl;
+
return lifetime;
}
diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h
index cd96f83d5af25e4be078d693458ad3239d741aaf..604988a1da2a1d1b72d6bd1b4beb819924a7afc3 100644
--- a/Processes/Sibyll/ParticleConversion.h
+++ b/Processes/Sibyll/ParticleConversion.h
@@ -34,9 +34,8 @@ namespace corsika::process::sibyll {
}
SibyllCode constexpr ConvertToSibyll(corsika::particles::Code pCode) {
- //~ assert(handledBySibyll(pCode));
return static_cast<SibyllCode>(
- corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
+ corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
}
corsika::particles::Code constexpr ConvertFromSibyll(SibyllCode pCode) {
@@ -44,8 +43,7 @@ namespace corsika::process::sibyll {
}
int constexpr ConvertToSibyllRaw(corsika::particles::Code pCode) {
- return (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>(
- corsika::process::sibyll::ConvertToSibyll(pCode));
+ return static_cast<int>(ConvertToSibyll(pCode));
}
int constexpr GetSibyllXSCode(corsika::particles::Code pCode) {
diff --git a/Processes/Sibyll/sibyll2.3c.cc b/Processes/Sibyll/sibyll2.3c.cc
index 77ad3052724c066fe474c1d41d44cb013f0786a4..fefea6d0b1acbac71e804cc83234bfe5bb335248 100644
--- a/Processes/Sibyll/sibyll2.3c.cc
+++ b/Processes/Sibyll/sibyll2.3c.cc
@@ -1,10 +1,12 @@
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/random/RNGManager.h>
+#include <random>
double s_rndm_(int&) {
- static corsika::random::RNG& rmng =
+ static corsika::random::RNG& rng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
- ;
- return rmng() / (double)rmng.max();
+
+ std::uniform_real_distribution<double> dist;
+ return dist(rng);
}
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 3e9bfe831d1dfbe3d3d16fb1f67347dcaf0e51fa..f07bac485f36e6b1144430d249178759dda695a4 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -45,15 +45,19 @@ namespace corsika::stack {
void SetPID(const corsika::particles::Code id) {
GetStackData().SetPID(GetIndex(), id);
}
+
void SetEnergy(const corsika::units::hep::EnergyType& e) {
GetStackData().SetEnergy(GetIndex(), e);
}
+
void SetMomentum(const MomentumVector& v) {
GetStackData().SetMomentum(GetIndex(), v);
}
+
void SetPosition(const corsika::geometry::Point& v) {
GetStackData().SetPosition(GetIndex(), v);
}
+
void SetTime(const corsika::units::si::TimeType& v) {
GetStackData().SetTime(GetIndex(), v);
}
@@ -61,15 +65,19 @@ namespace corsika::stack {
corsika::particles::Code GetPID() const {
return GetStackData().GetPID(GetIndex());
}
+
corsika::units::hep::EnergyType GetEnergy() const {
return GetStackData().GetEnergy(GetIndex());
}
+
MomentumVector GetMomentum() const {
return GetStackData().GetMomentum(GetIndex());
}
+
corsika::geometry::Point GetPosition() const {
return GetStackData().GetPosition(GetIndex());
}
+
corsika::units::si::TimeType GetTime() const {
return GetStackData().GetTime(GetIndex());
}
@@ -99,22 +107,30 @@ namespace corsika::stack {
}
int GetSize() const { return fDataPID.size(); }
+
int GetCapacity() const { return fDataPID.size(); }
void SetPID(const int i, const corsika::particles::Code id) { fDataPID[i] = id; }
+
void SetEnergy(const int i, const corsika::units::hep::EnergyType e) {
fDataE[i] = e;
}
void SetMomentum(const int i, const MomentumVector& v) { fMomentum[i] = v; }
+
void SetPosition(const int i, const corsika::geometry::Point& v) {
fPosition[i] = v;
}
+
void SetTime(const int i, const corsika::units::si::TimeType& v) { fTime[i] = v; }
corsika::particles::Code GetPID(const int i) const { return fDataPID[i]; }
+
corsika::units::hep::EnergyType GetEnergy(const int i) const { return fDataE[i]; }
+
MomentumVector GetMomentum(const int i) const { return fMomentum[i]; }
+
corsika::geometry::Point GetPosition(const int i) const { return fPosition[i]; }
+
corsika::units::si::TimeType GetTime(const int i) const { return fTime[i]; }
/**
@@ -156,6 +172,7 @@ namespace corsika::stack {
0 * corsika::units::si::meter}));
fTime.push_back(0 * corsika::units::si::second);
}
+
void DecrementSize() {
if (fDataE.size() > 0) {
fDataPID.pop_back();