From ea7b3832ceda7c2f0217ed077b862e03ab56d004 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 5 May 2021 14:26:19 +0200
Subject: [PATCH] kinetic energy
---
corsika/detail/framework/stack/Stack.inl | 15 +--
.../modules/energy_loss/BetheBlochPDG.inl | 16 +--
.../detail/modules/proposal/Interaction.inl | 2 +-
corsika/detail/modules/pythia8/Decay.inl | 8 +-
.../detail/modules/pythia8/Interaction.inl | 3 +-
.../detail/modules/qgsjetII/Interaction.inl | 13 +-
corsika/detail/modules/sibyll/Decay.inl | 3 +-
corsika/detail/modules/sibyll/Interaction.inl | 5 +-
.../modules/sibyll/NuclearInteraction.inl | 23 ++--
corsika/detail/modules/urqmd/UrQMD.inl | 4 +-
.../detail/stack/NuclearStackExtension.inl | 119 +++++++++++-------
corsika/detail/stack/VectorStack.inl | 61 +++++----
corsika/framework/stack/Stack.hpp | 2 +-
.../stack/StackIteratorInterface.hpp | 7 +-
corsika/stack/NuclearStackExtension.hpp | 92 ++++++++++----
corsika/stack/VectorStack.hpp | 103 +++++++++------
examples/boundary_example.cpp | 2 +-
examples/cascade_example.cpp | 3 +-
examples/cascade_proton_example.cpp | 2 +-
examples/em_shower.cpp | 2 +-
examples/hybrid_MC.cpp | 4 +-
examples/stack_example.cpp | 4 +-
examples/stopping_power.cpp | 2 +-
examples/vertical_EAS.cpp | 15 +--
tests/common/SetupTestStack.hpp | 10 +-
tests/framework/testCascade.cpp | 10 +-
tests/modules/testParticleCut.cpp | 81 ++++++------
tests/modules/testPythia8.cpp | 11 +-
tests/modules/testStackInspector.cpp | 4 +-
tests/stack/testHistoryView.cpp | 36 +++---
tests/stack/testNuclearStackExtension.cpp | 84 ++++++-------
tests/stack/testVectorStack.cpp | 24 ++--
32 files changed, 421 insertions(+), 349 deletions(-)
diff --git a/corsika/detail/framework/stack/Stack.inl b/corsika/detail/framework/stack/Stack.inl
index 4569a23a8..2929f93fe 100644
--- a/corsika/detail/framework/stack/Stack.inl
+++ b/corsika/detail/framework/stack/Stack.inl
@@ -141,7 +141,6 @@ namespace corsika {
template <typename... TArgs>
typename Stack<StackData, MParticleInterface>::stack_iterator_type inline Stack<
StackData, MParticleInterface>::addParticle(const TArgs... v) {
- CORSIKA_LOG_TRACE("Stack::AddParticle");
data_.incrementSize();
deleted_.push_back(false);
return stack_iterator_type(*this, getSize() - 1, v...);
@@ -150,21 +149,18 @@ namespace corsika {
template <typename StackData, template <typename> typename MParticleInterface>
inline void Stack<StackData, MParticleInterface>::swap(stack_iterator_type a,
stack_iterator_type b) {
- CORSIKA_LOG_TRACE("Stack::Swap");
swap(a.getIndex(), b.getIndex());
}
template <typename StackData, template <typename> typename MParticleInterface>
inline void Stack<StackData, MParticleInterface>::copy(stack_iterator_type a,
stack_iterator_type b) {
- CORSIKA_LOG_TRACE("Stack::Copy");
copy(a.getIndex(), b.getIndex());
}
template <typename StackData, template <typename> typename MParticleInterface>
inline void Stack<StackData, MParticleInterface>::copy(const_stack_iterator_type a,
stack_iterator_type b) {
- CORSIKA_LOG_TRACE("Stack::Copy");
data_.copy(a.getIndex(), b.getIndex());
if (deleted_[b.getIndex()] && !deleted_[a.getIndex()]) nDeleted_--;
if (!deleted_[b.getIndex()] && deleted_[a.getIndex()]) nDeleted_++;
@@ -173,11 +169,10 @@ namespace corsika {
template <typename StackData, template <typename> typename MParticleInterface>
inline void Stack<StackData, MParticleInterface>::erase(stack_iterator_type p) {
- CORSIKA_LOG_TRACE("Stack::Delete");
- if (this->isEmpty()) { /*error*/
+ if (this->isEmpty()) {
throw std::runtime_error("Stack, cannot delete entry since size is zero");
}
- if (deleted_[p.getIndex()]) { /*error*/
+ if (deleted_[p.getIndex()]) {
throw std::runtime_error("Stack, cannot delete entry since already deleted");
}
this->erase(p.getIndex());
@@ -231,7 +226,7 @@ namespace corsika {
if (!deleted_.back())
return false; // the last particle is not marked for deletion. Do nothing.
- CORSIKA_LOG_TRACE("Stack::purgeLastIfDeleted: yes");
+ CORSIKA_LOG_TRACE("stack: purgeLastIfDeleted: yes");
data_.decrementSize();
nDeleted_--;
deleted_.pop_back();
@@ -290,7 +285,6 @@ namespace corsika {
inline typename Stack<StackData, MParticleInterface>::stack_iterator_type
Stack<StackData, MParticleInterface>::addSecondary(stack_iterator_type& parent,
const TArgs... v) {
- CORSIKA_LOG_TRACE("Stack::AddSecondary");
data_.incrementSize();
deleted_.push_back(false);
return stack_iterator_type(*this, getSize() - 1, parent, v...);
@@ -299,7 +293,6 @@ namespace corsika {
template <typename StackData, template <typename> typename MParticleInterface>
inline void Stack<StackData, MParticleInterface>::swap(unsigned int const a,
unsigned int const b) {
- CORSIKA_LOG_TRACE("Stack::Swap(unsigned int)");
data_.swap(a, b);
std::swap(deleted_[a], deleted_[b]);
}
@@ -307,7 +300,6 @@ namespace corsika {
template <typename StackData, template <typename> typename MParticleInterface>
inline void Stack<StackData, MParticleInterface>::copy(unsigned int const a,
unsigned int const b) {
- CORSIKA_LOG_TRACE("Stack::Copy");
data_.copy(a, b);
if (deleted_[b] && !deleted_[a]) nDeleted_--;
if (!deleted_[b] && deleted_[a]) nDeleted_++;
@@ -346,7 +338,6 @@ namespace corsika {
template <typename StackData, template <typename> typename MParticleInterface>
inline unsigned int Stack<StackData, MParticleInterface>::getIndexFromIterator(
const unsigned int vI) const {
- // this is too much: CORSIKA_LOG_TRACE("Stack::getIndexFromIterator({})={}", vI, vI);
return vI;
}
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 693814c4d..31c34ce7b 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -179,17 +179,17 @@ namespace corsika {
auto constexpr dX = 1_g / square(1_cm);
auto const dEdX = -getTotalEnergyLoss(vParticle, dX) / dX;
- auto const energy = vParticle.getEnergy();
+ auto const energy = vParticle.getKineticEnergy();
auto const energy_lim =
std::max(energy * 0.9, // either 10% relative loss max., or
- get_kinetic_energy_threshold(vParticle.getPID()) +
- get_mass(vParticle.getPID()) // energy thresholds globally defined
- // for individual particles
- * 0.99999 // need to go slightly below global e-cut to assure removal in
- // ParticleCut. The 1% does not matter since at cut-time
- // the entire energy is removed.
+ get_kinetic_energy_threshold(
+ vParticle.getPID()) // energy thresholds globally defined
+ // for individual particles
+ * 0.99999 // need to go slightly below global e-cut to assure
+ // removal in ParticleCut. The 1% does not matter since
+ // at cut-time the entire energy is removed.
);
- auto const maxGrammage = (vParticle.getEnergy() - energy_lim) / dEdX;
+ auto const maxGrammage = (energy - energy_lim) / dEdX;
return vParticle.getNode()->getModelProperties().getArclengthFromGrammage(
vTrack, maxGrammage);
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index 868b00810..26aa7904c 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -97,7 +97,7 @@ namespace corsika::proposal {
auto sec_code =
convert_from_PDG(static_cast<PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
view.addSecondary(
- make_tuple(sec_code, E, p, projectile.getPosition(), projectile.getTime()));
+ make_tuple(sec_code, p, projectile.getPosition(), projectile.getTime()));
}
}
return ProcessReturn::Ok;
diff --git a/corsika/detail/modules/pythia8/Decay.inl b/corsika/detail/modules/pythia8/Decay.inl
index eb8cbd2e5..27fd2a31f 100644
--- a/corsika/detail/modules/pythia8/Decay.inl
+++ b/corsika/detail/modules/pythia8/Decay.inl
@@ -204,7 +204,8 @@ namespace corsika::pythia8 {
double const m = en;
// add particle to pythia stack
- event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
+ event.append(pdgCode, 1, 0, 0, // PID, status, col, acol
+ px, py, pz, en, m);
if (!Pythia8::Pythia::next())
throw std::runtime_error("Pythia::Decay: decay failed!");
@@ -232,9 +233,8 @@ namespace corsika::pythia8 {
fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
fourMomLab.getTimeLikeComponent());
- view.addSecondary(std::make_tuple(pyId, fourMomLab.getTimeLikeComponent(),
- fourMomLab.getSpaceLikeComponents(), decayPoint,
- t0));
+ view.addSecondary(
+ std::make_tuple(pyId, fourMomLab.getSpaceLikeComponents(), decayPoint, t0));
}
// set particle stable
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 4b74ed67f..4f288918b 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -364,11 +364,10 @@ namespace corsika::pythia8 {
MomentumVector const pyPlab(
labCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
- HEPEnergyType const pyEn = p8p.e() * 1_GeV;
// add to corsika stack
auto pnew =
- projectile.addSecondary(std::make_tuple(pyId, pyEn, pyPlab, pOrig, tOrig));
+ projectile.addSecondary(std::make_tuple(pyId, pyPlab, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl
index 9a86aa880..f442b239c 100644
--- a/corsika/detail/modules/qgsjetII/Interaction.inl
+++ b/corsika/detail/modules/qgsjetII/Interaction.inl
@@ -329,14 +329,13 @@ namespace corsika::qgsjetII {
sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
(projectileEnergyLabPerNucleon - nucleonMass))});
- auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleonMass));
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(
"secondary fragment> id= {}"
" p={}",
idFragm, momentum.getComponents());
- auto pnew = view.addSecondary(
- std::make_tuple(idFragm, energy, momentum, pOrig, tOrig));
+ auto pnew =
+ view.addSecondary(std::make_tuple(idFragm, momentum, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
} break;
@@ -363,7 +362,6 @@ namespace corsika::qgsjetII {
sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
(projectileEnergyLabPerNucleon * A - nucleusMass))});
- auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleusMass));
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(
"secondary fragment> id={}"
@@ -372,8 +370,8 @@ namespace corsika::qgsjetII {
" Z= {}",
idFragm, momentum.getComponents(), A, Z);
- auto pnew = view.addSecondary(
- std::make_tuple(idFragm, energy, momentum, pOrig, tOrig, A, Z));
+ auto pnew =
+ view.addSecondary(std::make_tuple(idFragm, momentum, pOrig, tOrig, A, Z));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
@@ -384,7 +382,6 @@ namespace corsika::qgsjetII {
for (auto& psec : qs) {
auto momentum = psec.getMomentum(zAxisFrame);
- auto const energy = psec.getEnergy();
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(
@@ -393,7 +390,7 @@ namespace corsika::qgsjetII {
corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
momentum.getComponents());
auto pnew = view.addSecondary(
- std::make_tuple(corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), energy,
+ std::make_tuple(corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
momentum, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
diff --git a/corsika/detail/modules/sibyll/Decay.inl b/corsika/detail/modules/sibyll/Decay.inl
index 5ca2edb85..1474df5c3 100644
--- a/corsika/detail/modules/sibyll/Decay.inl
+++ b/corsika/detail/modules/sibyll/Decay.inl
@@ -215,8 +215,7 @@ namespace corsika::sibyll {
if (psib.hasDecayed()) continue;
// add to corsika stack
projectile.addSecondary(std::make_tuple(sibyll::convertFromSibyll(psib.getPID()),
- psib.getEnergy(), psib.getMomentum(),
- decayPoint, t0));
+ psib.getMomentum(), decayPoint, t0));
}
// empty sibyll stack
ss.clear();
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index 69c4d4339..78b237e94 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -354,9 +354,8 @@ namespace corsika::sibyll {
assert(p3lab.getCoordinateSystem() == originalCS); // just to be sure!
// add to corsika stack
- auto pnew = view.addSecondary(
- std::make_tuple(corsika::sibyll::convertFromSibyll(psib.getPID()),
- Plab.getTimeLikeComponent(), p3lab, pOrig, tOrig));
+ auto pnew = view.addSecondary(std::make_tuple(
+ corsika::sibyll::convertFromSibyll(psib.getPID()), p3lab, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
diff --git a/corsika/detail/modules/sibyll/NuclearInteraction.inl b/corsika/detail/modules/sibyll/NuclearInteraction.inl
index 91b85c681..6e212db02 100644
--- a/corsika/detail/modules/sibyll/NuclearInteraction.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteraction.inl
@@ -546,14 +546,12 @@ namespace corsika::sibyll {
if (nuclA == 1)
// add nucleon
- projectile.addSecondary(std::make_tuple(specCode, Plab.getTimeLikeComponent(),
- Plab.getSpaceLikeComponents(), pOrig,
- tOrig));
+ projectile.addSecondary(
+ std::make_tuple(specCode, Plab.getSpaceLikeComponents(), pOrig, tOrig));
else
// add nucleus
- projectile.addSecondary(std::make_tuple(specCode, Plab.getTimeLikeComponent(),
- Plab.getSpaceLikeComponents(), pOrig,
- tOrig, nuclA, nuclZ));
+ projectile.addSecondary(std::make_tuple(specCode, Plab.getSpaceLikeComponents(),
+ pOrig, tOrig, nuclA, nuclZ));
}
// add elastic nucleons to corsika stack
@@ -570,9 +568,8 @@ namespace corsika::sibyll {
const double mass_ratio = get_mass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
- projectile.addSecondary(std::make_tuple(elaNucCode, Plab.getTimeLikeComponent(),
- Plab.getSpaceLikeComponents(), pOrig,
- tOrig));
+ projectile.addSecondary(
+ std::make_tuple(elaNucCode, Plab.getSpaceLikeComponents(), pOrig, tOrig));
}
// add inelastic interactions
@@ -584,8 +581,7 @@ namespace corsika::sibyll {
CORSIKA_LOG_DEBUG("inelastic interaction no. {}", j);
setup::Stack nucleonStack;
auto inelasticNucleon = nucleonStack.addParticle(
- std::make_tuple(pCode, PprojNucLab.getTimeLikeComponent(),
- PprojNucLab.getSpaceLikeComponents(), pOrig, tOrig));
+ std::make_tuple(pCode, PprojNucLab.getSpaceLikeComponents(), pOrig, tOrig));
inelasticNucleon.setNode(projectile.getNode());
// create inelastic interaction for each nucleon
CORSIKA_LOG_TRACE("calling HadronicInteraction...");
@@ -595,9 +591,8 @@ namespace corsika::sibyll {
hadronicInteraction_.doInteraction(nucleon_secondaries);
// inelasticNucleon.Delete(); // this is just a temporary object
for (const auto& pSec : nucleon_secondaries) {
- projectile.addSecondary(std::make_tuple(pSec.getPID(), pSec.getEnergy(),
- pSec.getMomentum(), pSec.getPosition(),
- pSec.getTime()));
+ projectile.addSecondary(std::make_tuple(pSec.getPID(), pSec.getMomentum(),
+ pSec.getPosition(), pSec.getTime()));
}
}
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index e6572b37b..e80662686 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -313,13 +313,11 @@ namespace corsika::urqmd {
::urqmd::coor_.px[i], ::urqmd::coor_.py[i], ::urqmd::coor_.pz[i]} *
1_GeV);
- auto const energy = sqrt(momentum.getSquaredNorm() + square(get_mass(code)));
-
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(" {} {} {} ", i, code, momentum.getComponents());
projectile.addSecondary(
- std::make_tuple(code, energy, momentum, projectilePosition, projectileTime));
+ std::make_tuple(code, momentum, projectilePosition, projectileTime));
}
CORSIKA_LOG_DEBUG("UrQMD generated {} secondaries!", ::urqmd::sys_.npart);
}
diff --git a/corsika/detail/stack/NuclearStackExtension.inl b/corsika/detail/stack/NuclearStackExtension.inl
index ca39c60eb..9bb5843d1 100644
--- a/corsika/detail/stack/NuclearStackExtension.inl
+++ b/corsika/detail/stack/NuclearStackExtension.inl
@@ -24,7 +24,7 @@ namespace corsika::nuclear_stack {
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(particle_data_type const& v) {
+ setParticleData(typename super_type::particle_data_type const& v) {
if (std::get<0>(v) == Code::Nucleus) {
std::ostringstream err;
@@ -39,10 +39,11 @@ namespace corsika::nuclear_stack {
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(altenative_particle_data_type const& v) {
- const unsigned short A = std::get<5>(v);
- const unsigned short Z = std::get<6>(v);
- if (std::get<0>(v) != Code::Nucleus || A == 0 || Z == 0) {
+ setParticleData(nuclear_particle_data_type const& v) {
+ unsigned short const A = std::get<5>(v);
+ unsigned short const Z = std::get<6>(v);
+ Code const PID = std::get<0>(v);
+ if (PID != Code::Nucleus || A == 0 || Z == 0) {
std::ostringstream err;
err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
throw std::runtime_error(err.str());
@@ -51,34 +52,31 @@ namespace corsika::nuclear_stack {
super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
setNuclearA(A);
setNuclearZ(Z);
- super_type::setParticleData(particle_data_type{
- std::get<0>(v), std::get<1>(v), std::get<2>(v), std::get<3>(v), std::get<4>(v)});
+ super_type::setParticleData(typename super_type::particle_data_type{
+ PID, std::get<1>(v), std::get<2>(v), std::get<3>(v), std::get<4>(v)});
}
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(super_type& p, particle_data_type const& v) {
+ setParticleData(super_type& p, typename super_type::particle_data_type const& v) {
if (std::get<0>(v) == Code::Nucleus) {
std::ostringstream err;
err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
throw std::runtime_error(err.str());
}
- super_type::setParticleData(
- p, particle_data_type{std::get<0>(v), std::get<1>(v), std::get<2>(v),
- std::get<3>(v), std::get<4>(v)});
-
+ super_type::setParticleData(p, v);
setNucleusRef(-1); // this is not a nucleus
}
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(super_type& p, altenative_particle_data_type const& v) {
+ setParticleData(super_type& p, nuclear_particle_data_type const& v) {
- const unsigned short A = std::get<5>(v);
- const unsigned short Z = std::get<6>(v);
+ unsigned short const A = std::get<5>(v);
+ unsigned short const Z = std::get<6>(v);
if (std::get<0>(v) != Code::Nucleus || A == 0 || Z == 0) {
std::ostringstream err;
@@ -90,15 +88,15 @@ namespace corsika::nuclear_stack {
super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
setNuclearA(A);
setNuclearZ(Z);
- super_type::setParticleData(
- p, particle_data_type{std::get<0>(v), std::get<1>(v), std::get<2>(v),
- std::get<3>(v), std::get<4>(v)});
+ super_type::setParticleData(p, typename super_type::particle_data_type{
+ std::get<0>(v), std::get<1>(v), std::get<2>(v),
+ std::get<3>(v), std::get<4>(v)});
}
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(particle_data_momentum_type const& v) {
+ setParticleData(typename super_type::particle_data_momentum_type const& v) {
if (std::get<0>(v) == Code::Nucleus) {
std::ostringstream err;
@@ -113,10 +111,11 @@ namespace corsika::nuclear_stack {
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(altenative_particle_data_momentum_type const& v) {
- const unsigned short A = std::get<5>(v);
- const unsigned short Z = std::get<6>(v);
- if (std::get<0>(v) != Code::Nucleus || A == 0 || Z == 0) {
+ setParticleData(nuclear_particle_data_momentum_type const& v) {
+ unsigned short const A = std::get<4>(v);
+ unsigned short const Z = std::get<5>(v);
+ Code const PID = std::get<0>(v);
+ if (PID != Code::Nucleus || A == 0 || Z == 0) {
std::ostringstream err;
err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
throw std::runtime_error(err.str());
@@ -125,36 +124,43 @@ namespace corsika::nuclear_stack {
super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
setNuclearA(A);
setNuclearZ(Z);
- super_type::setParticleData(particle_data_momentum_type{
- std::get<0>(v), std::get<1>(v), std::get<2>(v), std::get<3>(v), std::get<4>(v)});
+ HEPMassType m = 0_GeV;
+ if (PID == Code::Nucleus) {
+ m = get_nucleus_mass(A, Z);
+ } else {
+ m = get_mass(PID);
+ }
+ MomentumVector const& p = std::get<1>(v);
+ auto const P2 = p.getSquaredNorm();
+ HEPEnergyType Ekin = sqrt(P2 + square(m)) - m;
+ super_type::setParticleData(typename super_type::particle_data_type{
+ PID, Ekin, p / sqrt(P2), std::get<2>(v), std::get<3>(v)});
}
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(super_type& p, particle_data_momentum_type const& v) {
+ setParticleData(super_type& p,
+ typename super_type::particle_data_momentum_type const& v) {
if (std::get<0>(v) == Code::Nucleus) {
std::ostringstream err;
err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
throw std::runtime_error(err.str());
}
- super_type::setParticleData(
- p, particle_data_momentum_type{std::get<0>(v), std::get<1>(v), std::get<2>(v),
- std::get<3>(v), std::get<4>(v)});
-
+ super_type::setParticleData(p, v);
setNucleusRef(-1); // this is not a nucleus
}
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(super_type& p, altenative_particle_data_momentum_type const& v) {
-
- const unsigned short A = std::get<5>(v);
- const unsigned short Z = std::get<6>(v);
+ setParticleData(super_type& parent, nuclear_particle_data_momentum_type const& v) {
- if (std::get<0>(v) != Code::Nucleus || A == 0 || Z == 0) {
+ unsigned short const A = std::get<4>(v);
+ unsigned short const Z = std::get<5>(v);
+ Code const PID = std::get<0>(v);
+ if (PID != Code::Nucleus || A == 0 || Z == 0) {
std::ostringstream err;
err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
throw std::runtime_error(err.str());
@@ -164,9 +170,18 @@ namespace corsika::nuclear_stack {
super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
setNuclearA(A);
setNuclearZ(Z);
+ HEPMassType m = 0_GeV;
+ if (PID == Code::Nucleus) {
+ m = get_nucleus_mass(A, Z);
+ } else {
+ m = get_mass(PID);
+ }
+ MomentumVector const& p = std::get<1>(v);
+ auto const P2 = p.getSquaredNorm();
+ HEPEnergyType Ekin = sqrt(P2 + square(m)) - m;
super_type::setParticleData(
- p, particle_data_momentum_type{std::get<0>(v), std::get<1>(v), std::get<2>(v),
- std::get<3>(v), std::get<4>(v)});
+ parent, typename super_type::particle_data_type{PID, Ekin, p / sqrt(P2),
+ std::get<2>(v), std::get<3>(v)});
}
template <template <typename> class InnerParticleInterface,
@@ -185,12 +200,13 @@ namespace corsika::nuclear_stack {
MomentumVector const& v) {
HEPMomentumType const P = v.getNorm();
if (P == 0_eV) {
- super_type::getStackData().setEnergy(super_type::getIndex(), this->getMass());
+ super_type::getStackData().setKineticEnergy(super_type::getIndex(), 0_eV);
super_type::getStackData().setDirection(
super_type::getIndex(), DirectionVector(v.getCoordinateSystem(), {0, 0, 0}));
} else {
- super_type::getStackData().setEnergy(super_type::getIndex(),
- sqrt(square(this->getMass()) + square(P)));
+ super_type::getStackData().setKineticEnergy(
+ super_type::getIndex(),
+ sqrt(square(this->getMass()) + square(P)) - this->getMass());
super_type::getStackData().setDirection(super_type::getIndex(), v / P);
}
}
@@ -199,10 +215,18 @@ namespace corsika::nuclear_stack {
typename StackIteratorInterface>
inline MomentumVector NuclearParticleInterface<
InnerParticleInterface, StackIteratorInterface>::getMomentum() const {
- auto const P = sqrt(square(this->getEnergy()) - square(this->getMass()));
+ auto const P = sqrt(square(this->getKineticEnergy() + this->getMass()) -
+ square(this->getMass()));
return super_type::getStackData().getDirection(super_type::getIndex()) * P;
}
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline VelocityVector NuclearParticleInterface<
+ InnerParticleInterface, StackIteratorInterface>::getVelocity() const {
+ return this->getMomentum() / this->getEnergy() * constants::c;
+ }
+
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline HEPMassType NuclearParticleInterface<InnerParticleInterface,
@@ -223,8 +247,17 @@ namespace corsika::nuclear_stack {
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline HEPEnergyType NuclearParticleInterface<
- InnerParticleInterface, StackIteratorInterface>::getKineticEnergy() const {
- return this->getEnergy() - this->getMass();
+ InnerParticleInterface, StackIteratorInterface>::getEnergy() const {
+ return this->getKineticEnergy() + this->getMass();
+ }
+
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline void
+ NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::setEnergy(
+ HEPEnergyType const& e) {
+ super_type::getStackData().setKineticEnergy(super_type::getIndex(),
+ e - this->getMass());
}
template <template <typename> class InnerParticleInterface,
diff --git a/corsika/detail/stack/VectorStack.inl b/corsika/detail/stack/VectorStack.inl
index b950182bb..40fb4e61c 100644
--- a/corsika/detail/stack/VectorStack.inl
+++ b/corsika/detail/stack/VectorStack.inl
@@ -24,42 +24,44 @@ namespace corsika {
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v) {
+ particle_data_momentum_type const& v) {
this->setPID(std::get<0>(v));
- this->setEnergy(std::get<1>(v));
- MomentumVector const p = std::get<2>(v);
- HEPMomentumType const P = p.getNorm();
- if (P == 0_eV) {
+ MomentumVector const p = std::get<1>(v);
+ auto const P2 = p.getSquaredNorm();
+ HEPMassType const M = getMass();
+ this->setKineticEnergy(sqrt(P2 + square(M)) - M);
+ if (P2 == static_pow<2>(0_eV)) {
this->setDirection(DirectionVector(p.getCoordinateSystem(), {0, 0, 0}));
} else {
- this->setDirection(p / P);
+ this->setDirection(p / sqrt(P2));
}
- this->setPosition(std::get<3>(v));
- this->setTime(std::get<4>(v));
+ this->setPosition(std::get<2>(v));
+ this->setTime(std::get<3>(v));
}
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
ParticleInterface<StackIteratorInterface> const&,
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v) {
+ particle_data_momentum_type const& v) {
this->setPID(std::get<0>(v));
- this->setEnergy(std::get<1>(v));
- MomentumVector const p = std::get<2>(v);
- HEPMomentumType const P = p.getNorm();
- if (P == 0_eV) {
+ MomentumVector const p = std::get<1>(v);
+ auto const P2 = p.getSquaredNorm();
+ HEPMassType const M = getMass();
+ this->setKineticEnergy(sqrt(P2 + square(M)) - M);
+ if (P2 == static_pow<2>(0_eV)) {
this->setDirection(DirectionVector(p.getCoordinateSystem(), {0, 0, 0}));
} else {
- this->setDirection(p / P);
+ this->setDirection(p / sqrt(P2));
}
- this->setPosition(std::get<3>(v));
- this->setTime(std::get<4>(v));
+ this->setPosition(std::get<2>(v));
+ this->setTime(std::get<3>(v));
}
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
- std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType> const& v) {
+ particle_data_type const& v) {
this->setPID(std::get<0>(v));
- this->setEnergy(std::get<1>(v));
+ this->setKineticEnergy(std::get<1>(v));
this->setDirection(std::get<2>(v));
this->setPosition(std::get<3>(v));
this->setTime(std::get<4>(v));
@@ -67,18 +69,23 @@ namespace corsika {
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
- ParticleInterface<StackIteratorInterface> const&,
- std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType> const& v) {
+ ParticleInterface<StackIteratorInterface> const&, particle_data_type const& v) {
this->setPID(std::get<0>(v));
- this->setEnergy(std::get<1>(v));
+ this->setKineticEnergy(std::get<1>(v));
this->setDirection(std::get<2>(v));
this->setPosition(std::get<3>(v));
this->setTime(std::get<4>(v));
}
+ template <typename StackIteratorInterface>
+ inline std::string ParticleInterface<StackIteratorInterface>::asString() const {
+ return fmt::format("particle: i={}, PID={}, Ekin={}GeV", super_type::getIndex(),
+ get_name(this->getPID()), this->getKineticEnergy() / 1_GeV);
+ }
+
inline void VectorStackImpl::clear() {
dataPID_.clear();
- dataE_.clear();
+ dataEkin_.clear();
direction_.clear();
position_.clear();
time_.clear();
@@ -86,7 +93,7 @@ namespace corsika {
inline void VectorStackImpl::copy(size_t i1, size_t i2) {
dataPID_[i2] = dataPID_[i1];
- dataE_[i2] = dataE_[i1];
+ dataEkin_[i2] = dataEkin_[i1];
direction_[i2] = direction_[i1];
position_[i2] = position_[i1];
time_[i2] = time_[i1];
@@ -94,7 +101,7 @@ namespace corsika {
inline void VectorStackImpl::swap(size_t i1, size_t i2) {
std::swap(dataPID_[i2], dataPID_[i1]);
- std::swap(dataE_[i2], dataE_[i1]);
+ std::swap(dataEkin_[i2], dataEkin_[i1]);
std::swap(direction_[i2], direction_[i1]);
std::swap(position_[i2], position_[i1]);
std::swap(time_[i2], time_[i1]);
@@ -102,7 +109,7 @@ namespace corsika {
inline void VectorStackImpl::incrementSize() {
dataPID_.push_back(Code::Unknown);
- dataE_.push_back(0 * electronvolt);
+ dataEkin_.push_back(0 * electronvolt);
CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
@@ -113,9 +120,9 @@ namespace corsika {
}
inline void VectorStackImpl::decrementSize() {
- if (dataE_.size() > 0) {
+ if (dataEkin_.size() > 0) {
dataPID_.pop_back();
- dataE_.pop_back();
+ dataEkin_.pop_back();
direction_.pop_back();
position_.pop_back();
time_.pop_back();
diff --git a/corsika/framework/stack/Stack.hpp b/corsika/framework/stack/Stack.hpp
index 2ff173fb3..caaf9d065 100644
--- a/corsika/framework/stack/Stack.hpp
+++ b/corsika/framework/stack/Stack.hpp
@@ -255,7 +255,7 @@ namespace corsika {
* particle/projectile
*
* This should only get internally called from a
- * StackIterator::AddSecondary via ParticleBase
+ * StackIterator::addSecondary via ParticleBase
*/
template <typename... TArgs>
stack_iterator_type addSecondary(stack_iterator_type& parent, const TArgs... v);
diff --git a/corsika/framework/stack/StackIteratorInterface.hpp b/corsika/framework/stack/StackIteratorInterface.hpp
index 76815746a..02ef7e542 100644
--- a/corsika/framework/stack/StackIteratorInterface.hpp
+++ b/corsika/framework/stack/StackIteratorInterface.hpp
@@ -77,9 +77,6 @@ namespace corsika {
corsika::StackIteratorInterface<TStackData, TParticleInterface, TStackType>>
particle_interface_type;
- // using ParticleInterfaceType = TParticleInterface<
- // corsika::StackIteratorInterface<TStackData, TParticleInterface, TStackType>>;
-
// it is not allowed to create a "dangling" stack iterator
StackIteratorInterface() = delete; //! \todo check rule of five
@@ -112,7 +109,7 @@ namespace corsika {
@param index index on stack
@param args variadic list of data to initialize stack entry, this must be
consistent with the definition of the user-provided
- particle_interface_type::SetParticleData(...) function
+ particle_interface_type::setParticleData(...) function
*/
template <typename... TArgs>
StackIteratorInterface(TStackType& data, unsigned int const index,
@@ -130,7 +127,7 @@ namespace corsika {
counting, history, etc.
@param args variadic list of data to initialize stack entry, this must be
consistent with the definition of the user-provided
- particle_interface_type::SetParticleData(...) function
+ particle_interface_type::setParticleData(...) function
*/
template <typename... Args>
StackIteratorInterface(TStackType& data, unsigned int const index,
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index 350a343ca..1139d440b 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -47,35 +47,74 @@ namespace corsika::nuclear_stack {
typedef InnerParticleInterface<StackIteratorInterface> super_type;
public:
- typedef std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType>
- particle_data_type;
-
typedef std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType,
unsigned short, unsigned short>
- altenative_particle_data_type;
-
- void setParticleData(particle_data_type const& v);
+ nuclear_particle_data_type;
- void setParticleData(altenative_particle_data_type const& v);
+ typedef std::tuple<Code, MomentumVector, Point, TimeType, unsigned short,
+ unsigned short>
+ nuclear_particle_data_momentum_type;
- void setParticleData(super_type& p, particle_data_type const& v);
+ /**
+ *
+ * @param v which is a tuple containing: PID, kinetic Energy, DirectionVector,
+ * Position, Time
+ */
+ void setParticleData(typename super_type::particle_data_type const& v);
- void setParticleData(super_type& p, altenative_particle_data_type const& v);
+ /**
+ *
+ * @param v which is a tuple containing: PID, kinetic Energy, DirectionVector,
+ * Position, Time, A, Z
+ */
+ void setParticleData(nuclear_particle_data_type const& v);
- typedef std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType>
- particle_data_momentum_type;
+ /**
+ *
+ * @param p the parent particle
+ * @param v which is a tuple containing: PID, Momentum Vector, Position,
+ * Time
+ */
+ void setParticleData(super_type& p, typename super_type::particle_data_type const& v);
- typedef std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType,
- unsigned short, unsigned short>
- altenative_particle_data_momentum_type;
+ /**
+ *
+ * @param p the parent particle
+ * @param v which is a tuple containing: PID, Momentum Vector, Position,
+ * Time, A, Z
+ */
+ void setParticleData(super_type& p, nuclear_particle_data_type const& v);
- void setParticleData(particle_data_momentum_type const& v);
+ /**
+ *
+ * @param v which is a tuple containing: PID, Total Energy, MomentumVector, Position,
+ * Time
+ */
+ void setParticleData(typename super_type::particle_data_momentum_type const& v);
- void setParticleData(altenative_particle_data_momentum_type const& v);
+ /**
+ *
+ * @param v which is a tuple containing: PID, Total Energy, MomentumVector, Position,
+ * Time, A, Z
+ */
+ void setParticleData(nuclear_particle_data_momentum_type const& v);
- void setParticleData(super_type& p, particle_data_momentum_type const& v);
+ /**
+ *
+ * @param p parent particle
+ * @param v which is a tuple containing: PID, Total Energy, MomentumVector, Position,
+ * Time
+ */
+ void setParticleData(super_type& p,
+ typename super_type::particle_data_momentum_type const& v);
- void setParticleData(super_type& p, altenative_particle_data_momentum_type const& v);
+ /**
+ *
+ * @param p parent particle
+ * @param v which is a tuple containing: PID, Total Energy, MomentumVector, Position,
+ * Time, A, Z
+ */
+ void setParticleData(super_type& p, nuclear_particle_data_momentum_type const& v);
std::string asString() const;
@@ -107,11 +146,21 @@ namespace corsika::nuclear_stack {
* Overwrite normal setMomentum function with nuclear version
*/
void setMomentum(MomentumVector const& v);
-
+
/**
* Overwrite normal getMomentum function with nuclear version
*/
MomentumVector getMomentum() const;
+
+ /**
+ * Overwrite normal getEnergy function with nuclear version
+ */
+ void setEnergy(HEPEnergyType const& e);
+
+ /**
+ * Overwrite normal getVelocity function with nuclear version
+ */
+ VelocityVector getVelocity() const;
/**
* Overwrite normal getMass function with nuclear version
@@ -124,9 +173,9 @@ namespace corsika::nuclear_stack {
ElectricChargeType getCharge() const;
/**
- * Overwrite normal getKineticEnergy function with nuclear version
+ * Overwrite normal getEnergy function with nuclear version
*/
- HEPEnergyType getKineticEnergy() const;
+ HEPEnergyType getEnergy() const;
/**
* Overwirte normal getChargeNumber function with nuclear version
@@ -204,6 +253,7 @@ namespace corsika::nuclear_stack {
int getNuclearA(const unsigned int i) const { return nuclearA_[getNucleusRef(i)]; }
int getNuclearZ(const unsigned int i) const { return nuclearZ_[getNucleusRef(i)]; }
+
// this function will create new storage for Nuclear Properties, and return the
// reference to it
int getNucleusNextRef();
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index e6bc6a3a4..0c20a2361 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -33,41 +33,68 @@ namespace corsika {
typedef ParticleBase<StackIteratorInterface> super_type;
public:
- std::string asString() const {
- return fmt::format("particle: i={}, PID={}, E={}GeV", super_type::getIndex(),
- get_name(this->getPID()), this->getEnergy() / 1_GeV);
- }
+ typedef std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType>
+ particle_data_type;
+
+ typedef std::tuple<Code, MomentumVector, Point, TimeType> particle_data_momentum_type;
+
+ std::string asString() const;
/**
- MomentumVector is only used to determine the DirectionVector, the normalization is
- lost.
+ * Set data of new particle.
+ *
+ * @param v tuple containing: PID, Momentum Vector, Position, Time
+ *
+ * MomentumVector is only used to determine the DirectionVector, the normalization
+ * is lost.
*/
- void setParticleData(
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v);
+ void setParticleData(particle_data_type const& v);
/**
- MomentumVector is only used to determine the DirectionVector, the normalization is
- lost.
+ * Set data of new particle.
+ *
+ * @param p parent particle
+ * @param v tuple containing: PID, Momentum Vector, Position, Time
+ *
+ * MomentumVector is only used to determine the DirectionVector, the normalization
+ * is lost.
*/
- void setParticleData(
- ParticleInterface<StackIteratorInterface> const&,
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v);
+ void setParticleData(ParticleInterface<StackIteratorInterface> const& p,
+ particle_data_type const& v);
- void setParticleData(
- std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType> const& v);
+ /**
+ * Set data of new particle.
+ *
+ * @param v tuple containing: PID, kinetic Energy, Direction Vector, Position, Time
+ *
+ */
+ void setParticleData(particle_data_momentum_type const& v);
- void setParticleData(
- ParticleInterface<StackIteratorInterface> const&,
- std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType> const& v);
+ /**
+ * Set data of new particle.
+ *
+ * @param p parent particle
+ * @param v tuple containing: PID, kinetic Energy, Direction Vector, Position, Time
+ *
+ */
+ void setParticleData(ParticleInterface<StackIteratorInterface> const& p,
+ particle_data_momentum_type const& v);
///! Set particle corsika::Code
void setPID(Code const id) {
super_type::getStackData().setPID(super_type::getIndex(), id);
}
+
///! Set energy
void setEnergy(HEPEnergyType const& e) {
- super_type::getStackData().setEnergy(super_type::getIndex(), e);
+ super_type::getStackData().setKineticEnergy(super_type::getIndex(), e - this->getMass());
+ }
+
+ ///! Set kinetic energy
+ void setKineticEnergy(HEPEnergyType const& ekin) {
+ super_type::getStackData().setKineticEnergy(super_type::getIndex(), ekin);
}
+
/**
The MomentumVector v is used to determine the DirectionVector, and to update the
particle energy.
@@ -75,12 +102,12 @@ namespace corsika {
void setMomentum(MomentumVector const& v) {
HEPMomentumType const P = v.getNorm();
if (P == 0_eV) {
- super_type::getStackData().setEnergy(super_type::getIndex(), getMass());
+ super_type::getStackData().setKineticEnergy(super_type::getIndex(), 0_eV);
super_type::getStackData().setDirection(
super_type::getIndex(), DirectionVector(v.getCoordinateSystem(), {0, 0, 0}));
} else {
- super_type::getStackData().setEnergy(super_type::getIndex(),
- sqrt(square(getMass()) + square(P)));
+ super_type::getStackData().setKineticEnergy(
+ super_type::getIndex(), sqrt(square(getMass()) + square(P)) - this->getMass());
super_type::getStackData().setDirection(super_type::getIndex(), v / P);
}
}
@@ -101,14 +128,9 @@ namespace corsika {
Code getPID() const {
return super_type::getStackData().getPID(super_type::getIndex());
}
- //! Get energy
- HEPEnergyType getEnergy() const {
- return super_type::getStackData().getEnergy(super_type::getIndex());
- }
- //! Get momentum
- MomentumVector getMomentum() const {
- auto const P = sqrt(square(getEnergy()) - square(getMass()));
- return super_type::getStackData().getDirection(super_type::getIndex()) * P;
+ //! Get kinetic energy
+ HEPEnergyType getKineticEnergy() const {
+ return super_type::getStackData().getKineticEnergy(super_type::getIndex());
}
//! Get direction
DirectionVector getDirection() const {
@@ -127,15 +149,24 @@ namespace corsika {
*
* @{
*/
+ //! Get velocity
VelocityVector getVelocity() const {
return this->getMomentum() / this->getEnergy() * constants::c;
}
-
+ //! Get momentum
+ MomentumVector getMomentum() const {
+ auto const P = sqrt(square(getEnergy()) - square(this->getMass()));
+ return super_type::getStackData().getDirection(super_type::getIndex()) * P;
+ }
+ //! Get mass of particle
HEPMassType getMass() const { return get_mass(this->getPID()); }
+ //! Get electric charge
ElectricChargeType getCharge() const { return get_charge(this->getPID()); }
- HEPEnergyType getKineticEnergy() const { return this->getEnergy() - this->getMass(); }
+ //! Get kinetic energy
+ HEPEnergyType getEnergy() const { return this->getKineticEnergy() + this->getMass(); }
+
//! Get charge number
int16_t getChargeNumber() const { return get_charge_number(this->getPID()); }
///@}
@@ -152,7 +183,7 @@ namespace corsika {
public:
typedef std::vector<Code> code_vector_type;
- typedef std::vector<HEPEnergyType> energy_vector_type;
+ typedef std::vector<HEPEnergyType> kinetic_energy_vector_type;
typedef std::vector<Point> point_vector_type;
typedef std::vector<TimeType> time_vector_type;
typedef std::vector<DirectionVector> direction_vector_type;
@@ -175,13 +206,13 @@ namespace corsika {
unsigned int getCapacity() const { return dataPID_.size(); }
void setPID(size_t i, Code const id) { dataPID_[i] = id; }
- void setEnergy(size_t i, HEPEnergyType const& e) { dataE_[i] = e; }
+ void setKineticEnergy(size_t i, HEPEnergyType const& e) { dataEkin_[i] = e; }
void setDirection(size_t i, DirectionVector const& v) { direction_[i] = v; }
void setPosition(size_t i, Point const& v) { position_[i] = v; }
void setTime(size_t i, TimeType const& v) { time_[i] = v; }
Code getPID(size_t i) const { return dataPID_[i]; }
- HEPEnergyType getEnergy(size_t i) const { return dataE_[i]; }
+ HEPEnergyType getKineticEnergy(size_t i) const { return dataEkin_[i]; }
DirectionVector getDirection(size_t i) const { return direction_[i]; }
Point getPosition(size_t i) const { return position_[i]; }
TimeType getTime(size_t i) const { return time_[i]; }
@@ -202,7 +233,7 @@ namespace corsika {
private:
/// the actual memory to store particle data
code_vector_type dataPID_;
- energy_vector_type dataE_;
+ kinetic_energy_vector_type dataEkin_;
direction_vector_type direction_;
point_vector_type position_;
time_vector_type time_;
diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp
index ebb202597..3ce8df868 100644
--- a/examples/boundary_example.cpp
+++ b/examples/boundary_example.cpp
@@ -167,7 +167,7 @@ int main() {
beamCode, theta, phi, plab.getComponents() / 1_GeV);
// shoot particles from inside target out
Point pos(rootCS, 0_m, 0_m, 0_m);
- stack.addParticle(std::make_tuple(beamCode, E0, plab, pos, 0_ns));
+ stack.addParticle(std::make_tuple(beamCode, plab, pos, 0_ns));
}
// define air shower object, run simulation
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index f7758db0a..2852428d2 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -126,8 +126,7 @@ int main() {
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
- stack.addParticle(
- std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, nuclA, nuclZ));
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, nuclA, nuclZ));
}
// setup processes, decays and interactions
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index a02c914e4..71be6f4f9 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -113,7 +113,7 @@ int main() {
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
- stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
}
// setup processes, decays and interactions
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 606d72b41..12400723b 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -133,7 +133,7 @@ int main(int argc, char** argv) {
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
- stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
<< std::endl;
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index f1602f5ce..25603ad8a 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -163,10 +163,10 @@ int main(int argc, char** argv) {
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A != 1) {
- stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
} else {
- stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
}
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
diff --git a/examples/stack_example.cpp b/examples/stack_example.cpp
index 30187c520..41beda4f9 100644
--- a/examples/stack_example.cpp
+++ b/examples/stack_example.cpp
@@ -22,8 +22,8 @@ using namespace std;
void fill(VectorStack& s) {
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
for (int i = 0; i < 11; ++i) {
- s.addParticle(std::make_tuple(Code::Electron, 1.5_GeV * i,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
+ s.addParticle(std::make_tuple(Code::Electron, 1.5_GeV * i - Electron::mass,
+ DirectionVector(rootCS, {1, 0, 0}),
Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
}
}
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index d419efbb8..4c9e43efc 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -78,7 +78,7 @@ int main() {
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
- stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
auto const p = stack.getNextParticle();
HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm));
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index c652a2db8..d6f5e98f5 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -99,12 +99,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
// 3.total energy in GeV, 4.number of showers,
// 5.seed (0 by default to generate random values for all)
-int Point::count_construct_ = 0;
-int Point::count_copy_construct_ = 0;
-int Point::count_move_ = 0;
-int Point::count_assign_ = 0;
-int Point::count_destruct_ = 0;
-
int main(int argc, char** argv) {
logging::set_level(logging::level::info);
@@ -326,20 +320,20 @@ int main(int argc, char** argv) {
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A > 1) {
- stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
} else {
if (A == 1) {
if (Z == 1) {
- stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
} else if (Z == 0) {
- stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
} else {
std::cerr << "illegal parameters" << std::endl;
return EXIT_FAILURE;
}
} else {
- stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
}
}
@@ -436,5 +430,4 @@ int main(int argc, char** argv) {
output.endOfLibrary();
}
- Point::trace();
}
diff --git a/tests/common/SetupTestStack.hpp b/tests/common/SetupTestStack.hpp
index 3d5322d4a..d41f1aa27 100644
--- a/tests/common/SetupTestStack.hpp
+++ b/tests/common/SetupTestStack.hpp
@@ -45,18 +45,14 @@ namespace corsika::setup::testing {
MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
if (vProjectileType == Code::Nucleus) {
- auto constexpr mN = constants::nucleonMass;
- HEPEnergyType const E0 = sqrt(static_pow<2>(mN * vA) + pLab.getSquaredNorm());
- auto particle = stack->addParticle(
- std::make_tuple(Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ));
+ auto particle =
+ stack->addParticle(std::make_tuple(Code::Nucleus, pLab, origin, 0_ns, vA, vZ));
particle.setNode(vNodePtr);
return std::make_tuple(std::move(stack),
std::make_unique<setup::StackView>(particle));
} else { // not a nucleus
- HEPEnergyType const E0 =
- sqrt(static_pow<2>(get_mass(vProjectileType)) + pLab.getSquaredNorm());
auto particle =
- stack->addParticle(std::make_tuple(vProjectileType, E0, pLab, origin, 0_ns));
+ stack->addParticle(std::make_tuple(vProjectileType, pLab, origin, 0_ns));
particle.setNode(vNodePtr);
return std::make_tuple(std::move(stack),
std::make_unique<setup::StackView>(particle));
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index 44d277bae..1d5cb6be2 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -95,12 +95,12 @@ public:
template <typename TView>
void doInteraction(TView& view) {
- calls_++;
+ ++calls_;
auto vP = view.getProjectile();
- const HEPEnergyType E = vP.getEnergy();
- vP.addSecondary(std::make_tuple(vP.getPID(), E / 2, vP.getMomentum(),
+ const HEPEnergyType Ekin = vP.getKineticEnergy();
+ vP.addSecondary(std::make_tuple(vP.getPID(), Ekin / 2, vP.getMomentum().normalized(),
vP.getPosition(), vP.getTime()));
- vP.addSecondary(std::make_tuple(vP.getPID(), E / 2, vP.getMomentum(),
+ vP.addSecondary(std::make_tuple(vP.getPID(), Ekin / 2, vP.getMomentum().normalized(),
vP.getPosition(), vP.getTime()));
}
@@ -159,7 +159,7 @@ TEST_CASE("Cascade", "[Cascade]") {
TestCascadeStack stack;
stack.clear();
stack.addParticle(std::make_tuple(
- Code::Electron, E0,
+ Code::Electron,
MomentumVector(rootCS, {0_GeV, 0_GeV,
-sqrt(E0 * E0 - static_pow<2>(get_mass(Code::Electron)))}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index 475764c96..94ec3f27f 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -55,9 +55,9 @@ TEST_CASE("ParticleCut", "processes") {
CHECK(cut.getHadronKineticECut() == 20_GeV);
// add primary particle to stack
- auto particle = stack.addParticle(std::make_tuple(
- Code::Proton, Eabove, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+ auto particle = stack.addParticle(
+ std::make_tuple(Code::Proton, Eabove, DirectionVector(rootCS, {1, 0, 0}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
// view on secondary particles
setup::StackView view(particle);
// ref. to primary particle through the secondary view.
@@ -67,7 +67,7 @@ TEST_CASE("ParticleCut", "processes") {
// only cut is by species
for (auto proType : particleList)
projectile.addSecondary(std::make_tuple(
- proType, Eabove, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ proType, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
CHECK(view.getEntries() == 11);
CHECK(stack.getEntries() == 12);
@@ -83,9 +83,8 @@ TEST_CASE("ParticleCut", "processes") {
ParticleCut cut(20_GeV, true, false);
// add primary particle to stack
- auto particle = stack.addParticle(
- std::make_tuple(Code::Proton, Eabove,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ auto particle = stack.addParticle(std::make_tuple(
+ Code::Proton, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
// view on secondary particles
setup::StackView view(particle);
// ref. to primary particle through the secondary view.
@@ -95,22 +94,22 @@ TEST_CASE("ParticleCut", "processes") {
// only cut is by species
for (auto proType : particleList) {
projectile.addSecondary(std::make_tuple(
- proType, Eabove, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ proType, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
}
cut.doSecondaries(view);
CHECK(view.getEntries() == 10);
CHECK(cut.getNumberEmParticles() == 1);
- CHECK(cut.getEmEnergy() / 1_GeV == 1000);
+ CHECK(cut.getEmEnergy() / 1_GeV == Approx(1000. + 0.000511));
}
SECTION("cut low energy") {
ParticleCut cut(20_GeV, true, true);
// add primary particle to stack
- auto particle = stack.addParticle(
- std::make_tuple(Code::Proton, Eabove,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ auto particle = stack.addParticle(std::make_tuple(Code::Proton, Eabove - Proton::mass,
+ DirectionVector(rootCS, {1, 0, 0}),
+ point0, 0_ns));
// view on secondary particles
setup::StackView view(particle);
// ref. to primary particle through the secondary view.
@@ -120,15 +119,15 @@ TEST_CASE("ParticleCut", "processes") {
// only cut is by species
for (auto proType : particleList)
projectile.addSecondary(std::make_tuple(
- proType, Ebelow, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ proType, Ebelow, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
unsigned short A = 18;
unsigned short Z = 8;
projectile.addSecondary(std::make_tuple(Code::Nucleus, Eabove * A,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- point0, 0_ns, A, Z));
+ DirectionVector(rootCS, {1, 0, 0}), point0,
+ 0_ns, A, Z));
projectile.addSecondary(std::make_tuple(Code::Nucleus, Ebelow * A,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- point0, 0_ns, A, Z));
+ DirectionVector(rootCS, {1, 0, 0}), point0,
+ 0_ns, A, Z));
cut.doSecondaries(view);
@@ -140,32 +139,32 @@ TEST_CASE("ParticleCut", "processes") {
ParticleCut cut(5_MeV, 5_MeV, 5_GeV, 5_GeV, true);
// add primary particle to stack
- auto particle = stack.addParticle(
- std::make_tuple(Code::Proton, Eabove,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ auto particle = stack.addParticle(std::make_tuple(Code::Proton, Eabove - Proton::mass,
+ DirectionVector(rootCS, {1, 0, 0}),
+ point0, 0_ns));
// view on secondary particles
setup::StackView view(particle);
// ref. to primary particle through the secondary view.
// only this way the secondary view is populated
auto projectile = view.getProjectile();
// add secondaries
- projectile.addSecondary(std::make_tuple(Code::Photon, 3_MeV,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- point0, 0_ns));
- projectile.addSecondary(std::make_tuple(Code::Electron, 3_MeV,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- point0, 0_ns));
- projectile.addSecondary(std::make_tuple(Code::PiPlus, 4_GeV,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- point0, 0_ns));
+ projectile.addSecondary(std::make_tuple(
+ Code::Photon, 3_MeV, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
+ projectile.addSecondary(std::make_tuple(Code::Electron, 3_MeV - Electron::mass,
+ DirectionVector(rootCS, {1, 0, 0}), point0,
+ 0_ns));
+ projectile.addSecondary(std::make_tuple(Code::PiPlus, 4_GeV - PiPlus::mass,
+ DirectionVector(rootCS, {1, 0, 0}), point0,
+ 0_ns));
+
unsigned short A = 18;
unsigned short Z = 8;
projectile.addSecondary(std::make_tuple(Code::Nucleus, 4_GeV * A,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- point0, 0_ns, A, Z));
+ DirectionVector(rootCS, {1, 0, 0}), point0,
+ 0_ns, A, Z));
projectile.addSecondary(std::make_tuple(Code::Nucleus, 6_GeV * A,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- point0, 0_ns, A, Z));
+ DirectionVector(rootCS, {1, 0, 0}), point0,
+ 0_ns, A, Z));
cut.doSecondaries(view);
@@ -187,9 +186,8 @@ TEST_CASE("ParticleCut", "processes") {
const TimeType too_late = 1_s;
// add primary particle to stack
- auto particle = stack.addParticle(
- std::make_tuple(Code::Proton, Eabove,
- MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 1_ns));
+ auto particle = stack.addParticle(std::make_tuple(
+ Code::Proton, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, 1_ns));
// view on secondary particles
setup::StackView view(particle);
// ref. to primary particle through the secondary view.
@@ -198,9 +196,9 @@ TEST_CASE("ParticleCut", "processes") {
// add secondaries, all with energies above the threshold
// only cut is by time
for (auto proType : particleList) {
- projectile.addSecondary(
- std::make_tuple(proType, Eabove, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- point0, too_late));
+ projectile.addSecondary(std::make_tuple(proType, Eabove - get_mass(proType),
+ DirectionVector(rootCS, {1, 0, 0}), point0,
+ too_late));
}
cut.doSecondaries(view);
@@ -223,8 +221,9 @@ TEST_CASE("ParticleCut", "processes") {
// add particles, all with energies above the threshold
// only cut is by species
for (auto proType : particleList) {
- auto particle = stack.addParticle(std::make_tuple(
- proType, Eabove, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}), point0, 0_ns));
+ auto particle = stack.addParticle(
+ std::make_tuple(proType, Eabove - get_mass(proType),
+ DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
if (cut.doContinuous(particle, track) == ProcessReturn::ParticleAbsorbed) {
particle.erase();
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index a4a93126a..b2d683a88 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -102,9 +102,6 @@ TEST_CASE("PythiaInterface", "[processes]") {
SECTION("pythia decay") {
HEPEnergyType const P0 = 10_GeV;
- // HEPMomentumType const E0 = sqrt(P0*P0 + PiPlus::mass*PiPlus::mass);
-
- // feenableexcept(FE_INVALID); \todo how does this work nowadays...???
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::PiPlus, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
auto& stack = *stackPtr;
@@ -120,9 +117,13 @@ TEST_CASE("PythiaInterface", "[processes]") {
corsika::pythia8::Decay model(particleList);
+ CORSIKA_LOG_INFO("stack: {} {}", stack.asString(), particle.asString());
[[maybe_unused]] const TimeType time = model.getLifetime(particle);
+ double const gamma = particle.getEnergy() / get_mass(Code::PiPlus);
+ CHECK(time == get_lifetime(Code::PiPlus) * gamma);
model.doDecay(*secViewPtr);
- CHECK(stack.getEntries() == 3);
+ CORSIKA_LOG_INFO("piplus->{}", stack.asString());
+ CHECK(stack.getEntries() == 3); // piplus, muplu, numu
auto const pSum = sumMomentum(view, cs);
CHECK((pSum - plab).getNorm() / 1_GeV == Approx(0).margin(1e-4));
CHECK((pSum.getNorm() - plab.getNorm()) / 1_GeV == Approx(0).margin(1e-4));
@@ -153,8 +154,6 @@ TEST_CASE("PythiaInterface", "[processes]") {
SECTION("pythia interaction") {
- //! feenableexcept(FE_INVALID); \todo how does this work nowadays
-
// this will be a p-p collision at sqrts=3.5TeV -> no problem for pythia
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Proton, 0, 0, 7_TeV, (setup::Environment::BaseNodeType* const)nodePtr,
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
index 0b2a3ab7b..9905583ed 100644
--- a/tests/modules/testStackInspector.cpp
+++ b/tests/modules/testStackInspector.cpp
@@ -21,7 +21,7 @@
using namespace corsika;
-TEST_CASE("StackInspector", "[processes]") {
+TEST_CASE("StackInspector", "modules") {
logging::set_level(logging::level::info);
@@ -34,7 +34,7 @@ TEST_CASE("StackInspector", "[processes]") {
TestCascadeStack stack;
stack.clear();
HEPEnergyType E0 = 100_GeV;
- stack.addParticle(std::make_tuple(Code::Electron, E0,
+ stack.addParticle(std::make_tuple(Code::Electron,
MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
diff --git a/tests/stack/testHistoryView.cpp b/tests/stack/testHistoryView.cpp
index 6ee9300c3..56e118209 100644
--- a/tests/stack/testHistoryView.cpp
+++ b/tests/stack/testHistoryView.cpp
@@ -78,9 +78,9 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
TestStack stack;
// add primary particle
- auto p0 = stack.addParticle(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto p0 = stack.addParticle(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CHECK(stack.getEntries() == 1);
corsika::history::EventPtr evt = p0.getEvent();
@@ -99,9 +99,9 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
// add 5 secondaries
for (int i = 0; i < 5; ++i) {
- auto sec = hview0.addSecondary(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = hview0.addSecondary(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CHECK(sec.getParentEventIndex() == i);
CHECK(sec.getEvent() != nullptr);
@@ -120,9 +120,9 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
// add second generation of secondaries
// add 10 secondaries
for (int i = 0; i < 10; ++i) {
- auto sec = hview1.addSecondary(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = hview1.addSecondary(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CHECK(sec.getParentEventIndex() == i);
CHECK(sec.getEvent()->parentEvent() == ev1);
@@ -146,9 +146,9 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
for (int i = 0; i < 15; ++i) {
CORSIKA_LOG_TRACE("loop, view: " + std::to_string(i));
- auto sec = hview2.addSecondary(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = hview2.addSecondary(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CORSIKA_LOG_TRACE("loop, ---- ");
CHECK(sec.getParentEventIndex() == i);
@@ -175,9 +175,9 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
// add 5 secondaries
for (int i = 0; i < 5; ++i) {
CORSIKA_LOG_TRACE("loop " + std::to_string(i));
- auto sec = proj0.addSecondary(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = proj0.addSecondary(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CHECK(sec.getParentEventIndex() == i);
CHECK(sec.getEvent() != nullptr);
@@ -195,9 +195,9 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
// add second generation of secondaries
// add 10 secondaries
for (unsigned int i = 0; i < 10; ++i) {
- auto sec = proj1.addSecondary(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = proj1.addSecondary(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CHECK(sec.getParentEventIndex() == int(i));
CHECK(sec.getEvent()->parentEvent() == ev1);
diff --git a/tests/stack/testNuclearStackExtension.cpp b/tests/stack/testNuclearStackExtension.cpp
index f5a566813..fa9201023 100644
--- a/tests/stack/testNuclearStackExtension.cpp
+++ b/tests/stack/testNuclearStackExtension.cpp
@@ -17,18 +17,9 @@ using namespace corsika;
#include <iostream>
using namespace std;
-template <typename TParticle>
-HEPEnergyType kineticEnergy(TParticle const p) {
- if (p.getPID() == Code::Nucleus)
- return p.getEnergy() - get_nucleus_mass(p.getNuclearA(), p.getNuclearZ());
- else
- return p.getEnergy() - get_mass(p.getPID());
-}
-
-TEST_CASE("NuclearStackExtension", "[stack]") {
+TEST_CASE("NuclearStackExtension", "stack") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
@@ -36,9 +27,9 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
nuclear_stack::NuclearStackExtension<VectorStack,
nuclear_stack::ExtendedParticleInterfaceType>
s;
- s.addParticle(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ s.addParticle(
+ std::make_tuple(Code::Electron, 1.5_GeV, DirectionVector(dummyCS, {0, 1, 0}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CHECK(s.getEntries() == 1);
}
@@ -48,26 +39,26 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
nuclear_stack::ExtendedParticleInterfaceType>
s;
s.addParticle(std::make_tuple(
- Code::Nucleus, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Code::Nucleus, MomentumVector(dummyCS, {10_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 10));
CHECK(s.getEntries() == 1);
}
SECTION("write invalid nucleus") {
nuclear_stack::ParticleDataStack s;
- CHECK_THROWS(s.addParticle(std::make_tuple(
- Code::Nucleus, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 0, 0)));
+ CHECK_THROWS(s.addParticle(
+ std::make_tuple(Code::Nucleus, 1.5_GeV, DirectionVector(dummyCS, {1, 0, 0}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 0, 0)));
}
SECTION("read non nucleus") {
nuclear_stack::ParticleDataStack s;
- s.addParticle(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ s.addParticle(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {0_GeV, 10_GeV, 0_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
const auto pout = s.getNextParticle();
CHECK(pout.getPID() == Code::Electron);
- CHECK(pout.getEnergy() == 1.5_GeV);
+ CHECK(pout.getEnergy() == sqrt(square(10_GeV) + square(Electron::mass)));
CHECK(pout.getTime() == 100_s);
}
@@ -75,14 +66,14 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
auto const A = 10;
auto const Z = 9;
nuclear_stack::ParticleDataStack s;
- s.addParticle(std::make_tuple(
- Code::Nucleus, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, A, Z));
+ s.addParticle(
+ std::make_tuple(Code::Nucleus, 1.5_GeV, DirectionVector(dummyCS, {1, 0, 0}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, A, Z));
const auto pout = s.getNextParticle();
CHECK(pout.getPID() == Code::Nucleus);
- CHECK(pout.getEnergy() == 1.5_GeV);
+ CHECK(pout.getEnergy() == 1.5_GeV + pout.getMass());
CHECK(pout.getMass() == get_nucleus_mass(A, Z));
- CHECK(pout.getKineticEnergy() == kineticEnergy(pout));
+ CHECK(pout.getKineticEnergy() == 1.5_GeV);
CHECK(pout.getTime() == 100_s);
CHECK(pout.getNuclearA() == 10);
CHECK(pout.getNuclearZ() == 9);
@@ -90,9 +81,9 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
SECTION("read invalid nucleus") {
nuclear_stack::ParticleDataStack s;
- s.addParticle(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ s.addParticle(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 10_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
const auto pout = s.getNextParticle();
CHECK_THROWS(pout.getNuclearA());
CHECK_THROWS(pout.getNuclearZ());
@@ -105,11 +96,11 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
for (int i = 0; i < 99; ++i) {
if ((i + 1) % 10 == 0) {
s.addParticle(std::make_tuple(
- Code::Nucleus, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Code::Nucleus, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i / 2));
} else {
s.addParticle(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
}
}
@@ -126,13 +117,18 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
// i=9, 19, 29, etc. are nuclei
for (int i = 0; i < 99; ++i) {
if ((i + 1) % 10 == 0) {
+ unsigned int const A = i;
+ unsigned int const Z = i / 2;
s.addParticle(std::make_tuple(
- Code::Nucleus, i * 15_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i / 2));
+
+ Code::Nucleus, i * 15_GeV - get_nucleus_mass(A, Z),
+ DirectionVector(dummyCS, {0, 0, 1}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, A, Z));
} else {
- s.addParticle(std::make_tuple(
- Code::Electron, i * 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ s.addParticle(std::make_tuple(Code::Electron, i * 1.5_GeV - Electron::mass,
+ DirectionVector(dummyCS, {1, 0, 0}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}),
+ 100_s));
}
}
@@ -233,23 +229,23 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
// not valid:
CHECK_THROWS(s.addParticle(std::make_tuple(
- Code::Oxygen, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Code::Oxygen, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 16, 8)));
// valid
- auto particle = s.addParticle(std::make_tuple(
- Code::Nucleus, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9));
+ auto particle = s.addParticle(
+ std::make_tuple(Code::Nucleus, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9));
// not valid
CHECK_THROWS(particle.addSecondary(std::make_tuple(
- Code::Oxygen, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Code::Oxygen, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 16, 8)));
// add a another nucleus, so there are two now
- s.addParticle(std::make_tuple(
- Code::Nucleus, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9));
+ s.addParticle(
+ std::make_tuple(Code::Nucleus, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9));
// not valid, since end() is not a valid entry
CHECK_THROWS(s.swap(s.begin(), s.end()));
diff --git a/tests/stack/testVectorStack.cpp b/tests/stack/testVectorStack.cpp
index 074c34dea..7fa0903b4 100644
--- a/tests/stack/testVectorStack.cpp
+++ b/tests/stack/testVectorStack.cpp
@@ -17,32 +17,26 @@
using namespace corsika;
using namespace std;
-template <typename TParticle>
-HEPEnergyType kineticEnergy(TParticle const p) {
- return p.getEnergy() - get_mass(p.getPID());
-}
-
-TEST_CASE("VectorStack", "[stack]") {
+TEST_CASE("VectorStack", "stack") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
SECTION("read+write") {
VectorStack s;
- s.addParticle(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ s.addParticle(
+ std::make_tuple(Code::Electron, 1.5_GeV, DirectionVector(dummyCS, {1, 0, 0}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
// read
CHECK(s.getEntries() == 1);
CHECK(s.getSize() == 1);
auto pout = s.getNextParticle();
CHECK(pout.getPID() == Code::Electron);
- CHECK(pout.getEnergy() == 1.5_GeV);
- CHECK(pout.getKineticEnergy() == kineticEnergy(pout));
+ CHECK(pout.getEnergy() == 1.5_GeV + Electron::mass);
+ CHECK(pout.getKineticEnergy() == 1.5_GeV);
CHECK(pout.getTime() == 100_s);
}
@@ -51,9 +45,9 @@ TEST_CASE("VectorStack", "[stack]") {
VectorStack s;
for (int i = 0; i < 99; ++i)
- s.addParticle(std::make_tuple(
- Code::Electron, 1.5_GeV, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ s.addParticle(
+ std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CHECK(s.getSize() == 99);
--
GitLab