diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index b7b884cca012130b4590c26c3c83877b598092ab..232e7d002945f41a613266c0738015a2dd34717e 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -143,18 +143,18 @@ public:
cout << "removing em. particle..." << endl;
fEmEnergy += energy;
fEmCount += 1;
- p.Delete();
+ // p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += energy;
fInvCount += 1;
- p.Delete();
+ // p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += energy;
- p.Delete();
+ // p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
}
return ret;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 6693f764c5a7ddacce02ca7a20f720bcf1a0ddaf..225b58d64db655bb7467fd3698b2c122d3a8cc32 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -119,7 +119,10 @@ namespace corsika::cascade {
fProcessSequence.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
- // fStack.Delete(particle); // TODO: check if this is really needed
+ std::cout << "Cascade: delete absorbed particle " << particle.GetPID() << " "
+ << particle.GetEnergy() / 1_GeV << "GeV" << std::endl;
+ particle.Delete();
+ return;
} else {
std::cout << "sth. happening before geometric limit ?"
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 80ebf4b21b1999754a2ccf2d244e4dfd7e5e0b4f..081325f96038dc8c733f4599603e0ee97a22d04f 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -80,7 +80,7 @@ public:
}
template <typename Particle, typename T, typename Stack>
- void DoContinuous(Particle& p, T&, Stack& s) const {
+ EProcessReturn DoContinuous(Particle& p, T&, Stack& s) const {
EnergyType E = p.GetEnergy();
if (E < 85_MeV) {
p.Delete();
@@ -95,6 +95,7 @@ public:
pnew.SetPosition(p.GetPosition());
pnew.SetMomentum(p.GetMomentum());
}
+ return EProcessReturn::eOk;
}
void Init() { fCount = 0; }
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index 86779739c3e7f2ce480b95501fb29907a47659ce..dd9572ac94cbb8b24352df1c77a22227679ea1c2 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/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h
index ae4869ea84ba228e1b59b2fe89778f2e932efc25..215994b42fb74daaba3b382513d94e50a45d430b 100644
--- a/Framework/ProcessSequence/ProcessReturn.h
+++ b/Framework/ProcessSequence/ProcessReturn.h
@@ -20,13 +20,25 @@ namespace corsika::process {
that can be accumulated easily with "|="
*/
- enum class EProcessReturn {
- eOk = 1,
- eParticleAbsorbed = 2,
- eInteracted = 3,
- eDecayed = 4,
+ 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, EProcessReturn b) {
+ return a = a | b;
+ }
+
+ inline bool operator==(EProcessReturn a, EProcessReturn b) {
+ return (static_cast<int>(a) & static_cast<int>(b)) != 0;
+ }
+
} // namespace corsika::process
#endif
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index a1b92b0e3f9c27dd896012a16b56d0e063c3954f..a75ad6adc3d9aeae6159028f0d905683ff8d88d0 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -60,19 +60,21 @@ namespace corsika::process {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of<ContinuousProcess<T1>, T1>::value ||
is_process_sequence<T1>::value) {
- A.DoContinuous(p, t, s);
+ ret |= A.DoContinuous(p, t, s);
}
if constexpr (std::is_base_of<ContinuousProcess<T2>, T2>::value ||
is_process_sequence<T2>::value) {
- B.DoContinuous(p, t, s);
+ ret |= B.DoContinuous(p, t, s);
}
return ret;
}
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 4d3eca70f0df8eadde09469ad70ea326cdf7ece7..93516c141ad13888e47b57d72d247c897c2ebdbf 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -89,9 +89,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(
@@ -131,13 +130,14 @@ namespace corsika::process {
const corsika::units::si::TimeType t0 =
corsika::particles::GetLifetime(p.GetPID());
- cout << "Decay: GetLifetime: \n"
- << " code: " << p.GetPID() << endl;
- cout << " t0: " << t0 << endl;
- cout << " energy: " << E / 1_GeV << endl;
- cout << " gamma: " << gamma << endl;
- corsika::units::si::TimeType const lifetime = gamma * t0;
- cout << " -> tau: " << lifetime << endl;
+ auto const lifetime = gamma * t0;
+
+ cout << "Decay: code: " << p.GetPID() << endl;
+ cout << "Decay: MinStep: t0: " << t0 << endl;
+ cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
+ cout << "Decay: MinStep: gamma: " << gamma << endl;
+ cout << "Decay: MinStep: tau: " << lifetime << endl;
+
return lifetime;
}
diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h
index cd96f83d5af25e4be078d693458ad3239d741aaf..9bd417da2398d6eaec29b24d85242c1ad0d867ab 100644
--- a/Processes/Sibyll/ParticleConversion.h
+++ b/Processes/Sibyll/ParticleConversion.h
@@ -34,7 +34,6 @@ 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)]);
}
@@ -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..de48630e207e9a1831bd4ec9f34c4603333f8cfe 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..481f8d48dc1abab991086e63d1f3b0962313844f 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();