diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index 6c61177e3842f224ad31526fb7b32c31b7100412..d385ad5b93aa683b5270e8b3bf0414f7772df07f 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -13,11 +13,12 @@
namespace corsika {
- inline HEPEnergyType constexpr get_energy_threshold(Code const p) {
+ inline HEPEnergyType constexpr get_kinetic_energy_threshold(Code const p) {
return particle::detail::thresholds[static_cast<CodeIntType>(p)];
}
- inline void constexpr set_energy_threshold(Code const p, HEPEnergyType const val) {
+ inline void constexpr set_kinetic_energy_threshold(Code const p,
+ HEPEnergyType const val) {
particle::detail::thresholds[static_cast<CodeIntType>(p)] = val;
}
diff --git a/corsika/detail/framework/stack/Stack.inl b/corsika/detail/framework/stack/Stack.inl
index 4569a23a8ab49c16af4aa6a29a78759bcfc807a1..2929f93fe794f861e6ac6857fe3fc4c93305bc48 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/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
index 988bbe80b315d59344c8c71f8960ac5833e270aa..4562b39c0aa0fea4330fc686deab57d48683d98f 100644
--- a/corsika/detail/framework/utility/QuarticSolver.inl
+++ b/corsika/detail/framework/utility/QuarticSolver.inl
@@ -32,7 +32,7 @@ namespace corsika {
long double c3 = -b * b * e - d * d + 4. * c * e;
// cubic resolvent
- // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+ // y^3 − c*y^2 + (bd−4e)*y − b^2*e−d^2+4*c*e = 0
std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
long double y = x3[0]; // there is always at least one solution
@@ -52,6 +52,7 @@ namespace corsika {
q1 = q2 = y * 0.5;
// g1+g2 = b && g1+g2 = c-y <=> g^2 - b*g + c-y = 0 (p === g)
Det = b * b - 4 * (c - y);
+ CORSIKA_LOG_TRACE("Det={}", Det);
if (fabs(Det) < epsilon) { // in other words - D==0
p1 = p2 = b * 0.5;
} else {
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index 7eb690d64abf96e83af22fd255ce1814b315d2da..04aff2d69f0e237541c73389a10b5fc7823004b0 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -18,25 +18,27 @@ namespace corsika {
bool const inv)
: doCutEm_(false)
, doCutInv_(inv)
- , energy_(0_GeV)
- , em_energy_(0_GeV)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , energy_invcut_(0_GeV)
, em_count_(0)
- , inv_energy_(0_GeV)
, inv_count_(0) {
for (auto p : get_all_particles())
if (is_hadron(p))
- set_energy_threshold(p, eHadCut);
+ set_kinetic_energy_threshold(p, eHadCut);
else if (is_muon(p))
- set_energy_threshold(p, eMuCut);
+ set_kinetic_energy_threshold(p, eMuCut);
else if (p == Code::Electron || p == Code::Positron)
- set_energy_threshold(p, eEleCut);
+ set_kinetic_energy_threshold(p, eEleCut);
else if (p == Code::Photon)
- set_energy_threshold(p, ePhoCut);
+ set_kinetic_energy_threshold(p, ePhoCut);
else if (p == Code::Nucleus)
// nuclei have same threshold as hadrons on the nucleon level.
- set_energy_threshold(p, eHadCut);
+ set_kinetic_energy_threshold(p, eHadCut);
CORSIKA_LOG_DEBUG(
- "setting thresholds: electrons = {} GeV, photons = {} GeV, hadrons = {} GeV, "
+ "setting kinetic energy thresholds: electrons = {} GeV, photons = {} GeV, "
+ "hadrons = {} GeV, "
"muons = {} GeV",
eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV);
printThresholds();
@@ -46,17 +48,18 @@ namespace corsika {
bool const inv)
: doCutEm_(true)
, doCutInv_(inv)
- , energy_(0_GeV)
- , em_energy_(0_GeV)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , energy_invcut_(0_GeV)
, em_count_(0)
- , inv_energy_(0_GeV)
, inv_count_(0) {
for (auto p : get_all_particles())
if (is_hadron(p))
- set_energy_threshold(p, eHadCut);
+ set_kinetic_energy_threshold(p, eHadCut);
else if (is_muon(p))
- set_energy_threshold(p, eMuCut);
+ set_kinetic_energy_threshold(p, eMuCut);
CORSIKA_LOG_DEBUG(
"setting thresholds: hadrons = {} GeV, "
"muons = {} GeV",
@@ -67,13 +70,15 @@ namespace corsika {
inline ParticleCut::ParticleCut(HEPEnergyType const eCut, bool const em, bool const inv)
: doCutEm_(em)
, doCutInv_(inv)
- , energy_(0_GeV)
- , em_energy_(0_GeV)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , energy_invcut_(0_GeV)
, em_count_(0)
- , inv_energy_(0_GeV)
, inv_count_(0) {
- for (auto p : get_all_particles()) set_energy_threshold(p, eCut);
- CORSIKA_LOG_DEBUG("setting threshold for all particles to {} GeV", eCut / 1_GeV);
+ for (auto p : get_all_particles()) set_kinetic_energy_threshold(p, eCut);
+ CORSIKA_LOG_DEBUG("setting kinetic energy threshold for all particles to {} GeV",
+ eCut / 1_GeV);
printThresholds();
}
@@ -82,27 +87,28 @@ namespace corsika {
bool const inv)
: doCutEm_(em)
, doCutInv_(inv)
- , energy_(0_GeV)
- , em_energy_(0_GeV)
+ , energy_cut_(0_GeV)
+ , energy_timecut_(0_GeV)
+ , energy_emcut_(0_GeV)
+ , energy_invcut_(0_GeV)
, em_count_(0)
- , inv_energy_(0_GeV)
, inv_count_(0) {
- set_energy_thresholds(eCuts);
+ set_kinetic_energy_thresholds(eCuts);
CORSIKA_LOG_DEBUG("setting threshold particles individually");
printThresholds();
}
template <typename TParticle>
inline bool ParticleCut::isBelowEnergyCut(TParticle const& vP) const {
- auto const energyLab = vP.getEnergy();
+ auto const energyLab = vP.getKineticEnergy();
auto const pid = vP.getPID();
// nuclei
if (pid == Code::Nucleus) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / vP.getNuclearA();
- return (ElabNuc < get_energy_threshold(pid));
+ return (ElabNuc < get_kinetic_energy_threshold(pid));
} else {
- return (energyLab < get_energy_threshold(pid));
+ return (energyLab < get_kinetic_energy_threshold(pid));
}
}
@@ -114,26 +120,29 @@ namespace corsika {
inline bool ParticleCut::checkCutParticle(TParticle const& particle) {
Code const pid = particle.getPID();
- HEPEnergyType energy = particle.getEnergy();
- CORSIKA_LOG_DEBUG("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV", pid,
- energy / 1_GeV, (em_energy_ + inv_energy_ + energy_) / 1_GeV);
+ HEPEnergyType const kine_energy = particle.getKineticEnergy();
+ HEPEnergyType const energy = particle.getEnergy();
+ CORSIKA_LOG_DEBUG(
+ "ParticleCut: checking {}, E_kin= {} GeV, EcutTot={} GeV", pid,
+ kine_energy / 1_GeV,
+ (energy_emcut_ + energy_invcut_ + energy_cut_ + energy_timecut_) / 1_GeV);
if (doCutEm_ && is_em(pid)) {
CORSIKA_LOG_DEBUG("removing em. particle...");
- em_energy_ += energy;
+ energy_emcut_ += energy;
em_count_ += 1;
return true;
} else if (doCutInv_ && is_neutrino(pid)) {
CORSIKA_LOG_DEBUG("removing inv. particle...");
- inv_energy_ += energy;
+ energy_invcut_ += energy;
inv_count_ += 1;
return true;
} else if (isBelowEnergyCut(particle)) {
CORSIKA_LOG_DEBUG("removing low en. particle...");
- energy_ += energy;
+ energy_cut_ += energy;
return true;
} else if (particle.getTime() > 10_ms) {
CORSIKA_LOG_DEBUG("removing OLD particle...");
- energy_ += energy;
+ energy_timecut_ += energy;
return true;
}
return false; // this particle will not be removed/cut
@@ -161,29 +170,33 @@ namespace corsika {
inline void ParticleCut::printThresholds() {
for (auto p : get_all_particles()) {
- auto const Eth = get_energy_threshold(p);
- CORSIKA_LOG_INFO("energy threshold for particle {} is {} GeV", p, Eth / 1_GeV);
+ auto const Eth = get_kinetic_energy_threshold(p);
+ CORSIKA_LOG_INFO("kinetic energy threshold for particle {} is {} GeV", p,
+ Eth / 1_GeV);
}
}
inline void ParticleCut::showResults() {
CORSIKA_LOG_INFO(
" ******************************\n"
- " energy in em. component (GeV): {} \n "
- " no. of em. particles injected: {} \n "
- " energy in inv. component (GeV): {} \n "
- " no. of inv. particles injected: {} \n "
- " energy below particle cut (GeV): {} \n"
+ " energy removed by cut of electromagnetic (GeV): {} \n "
+ " no. of em. particles removed : {} \n "
+ " energy removed by cut of invisible (GeV): {} \n "
+ " no. of invisible particles removed : {} \n "
+ " energy removed by kinetic energy cut (GeV): {} \n"
+ " energy removed by time cut (GeV): {} \n"
" ******************************",
- em_energy_ / 1_GeV, em_count_, inv_energy_ / 1_GeV, inv_count_, energy_ / 1_GeV);
+ energy_emcut_ / 1_GeV, em_count_, energy_invcut_ / 1_GeV, inv_count_,
+ energy_cut_ / 1_GeV, energy_timecut_ / 1_GeV);
}
inline void ParticleCut::reset() {
- em_energy_ = 0_GeV;
+ energy_emcut_ = 0_GeV;
em_count_ = 0;
- inv_energy_ = 0_GeV;
+ energy_invcut_ = 0_GeV;
inv_count_ = 0;
- energy_ = 0_GeV;
+ energy_cut_ = 0_GeV;
+ energy_timecut_ = 0_GeV;
}
} // namespace corsika
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 310613d24d2d8fcbd75402822773891b8374dfca..31c34ce7b970dab276dad90a942a718422588ea1 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -179,16 +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_lim = std::max(
- energy * 0.9, // either 10% relative loss max., or
- get_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. This does not matter since at cut-time the entire
- // energy is removed.
- );
- auto const maxGrammage = (vParticle.getEnergy() - energy_lim) / dEdX;
+ 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()) // 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 = (energy - energy_lim) / dEdX;
return vParticle.getNode()->getModelProperties().getArclengthFromGrammage(
vTrack, maxGrammage);
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 975ff716650e137c485beee9c677e83ef0f46d3e..91308390ba4f9b5b8e87bb0d94e79646e9a6f101 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -33,8 +33,9 @@ namespace corsika::proposal {
// interpolate the crosssection for given media and energy cut. These may
// take some minutes if you have to build the tables and cannot read the
// from disk
- auto const emCut = get_energy_threshold(
- code); //! energy thresholds globally defined for individual particles
+ auto const emCut =
+ get_kinetic_energy_threshold(code) +
+ get_mass(code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Use higland multiple scattering and deactivate stochastic deflection by
@@ -129,11 +130,11 @@ namespace corsika::proposal {
auto const energy = vP.getEnergy();
auto const energy_lim =
std::max(energy * 0.9, // either 10% relative loss max., or
- get_energy_threshold(
+ get_kinetic_energy_threshold(
code) // energy thresholds globally defined for individual particles
- * 0.9999 // need to go 1% below global e-cut to assure removal in
- // ParticleCut. The 1% does not matter since at cut-time the
- // entire energy is removed.
+ * 0.9999 // need to go slightly below global e-cut to assure removal
+ // in ParticleCut. This does not matter since at cut-time
+ // the entire energy is removed.
);
// solving the track integral for giving energy lim
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index cd99ca9789ef51c371cbb5dc76e06bf3bdc25f7e..019d31270abeb1da6fe2a99121081363762dcac6 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -36,8 +36,9 @@ namespace corsika::proposal {
// interpolate the crosssection for given media and energy cut. These may
// take some minutes if you have to build the tables and cannot read the
// from disk
- auto const emCut = get_energy_threshold(
- code); //! energy thresholds globally defined for individual particles
+ auto const emCut =
+ get_kinetic_energy_threshold(code) +
+ get_mass(code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
@@ -98,8 +99,8 @@ namespace corsika::proposal {
vecProposal.GetZ() * E);
auto p = MomentumVector(labCS, vec);
auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
- view.addSecondary(std::make_tuple(sec_code, E, p, projectile.getPosition(),
- projectile.getTime()));
+ view.addSecondary(
+ std::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 9829283cfa01bebaf0aa1f032d35061adde11adb..163ec37a7864fdd2b2c85022aaafce73e73d6739 100644
--- a/corsika/detail/modules/pythia8/Decay.inl
+++ b/corsika/detail/modules/pythia8/Decay.inl
@@ -206,7 +206,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);
// LCOV_EXCL_START, we don't validate pythia8 internals
if (!Pythia8::Pythia::next())
@@ -238,9 +239,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 66afda1e36206be59616f55761ceb8224927ff24..74a910de769ff259aeafacc46e9b2d74d02e32e5 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -373,11 +373,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 9a86aa880e2333d9247bd2e18d347a83d9184541..f442b239cdefab4f237594d3659d30e91f362d1c 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 5ca2edb854bc24aaf32b04a9effeaae7c35b98cf..1474df5c33272409608134549b3c1ad7a981c9b2 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 1a06b8f95df9f9f62659ef68caeb84e56033940e..bbcacf4e5e7e90fca3fa5a21a0aa639848ccdc3c 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -355,9 +355,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 91b85c68126ffbdebb2417867fcdd2420361fcbc..6e212db022aae8071e3e83d1a665ec3e8936b70a 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 e6572b37bff6bfcd7e7115a4dbb9ce215f244e35..e806626865655616b9e00e9d4183dbd54c322a30 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 5dab28ab2506181580fd510e28164c507b70b878..fa0d5941fb0d7a8fc38a6371bb5475a3015e3689 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,9 +88,100 @@ namespace corsika::nuclear_stack {
super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
setNuclearA(A);
setNuclearZ(Z);
+ 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(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(v);
+ setNucleusRef(-1); // this is not a nucleus
+ }
+
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
+ 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());
+ }
+ setNucleusRef(
+ 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(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,
+ 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, v);
+ setNucleusRef(-1); // this is not a nucleus
+ }
+
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
+ setParticleData(super_type& parent, 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());
+ }
+
+ setNucleusRef(
+ 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_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,
@@ -104,6 +193,48 @@ namespace corsika::nuclear_stack {
(isNucleus() ? fmt::format("A={}, Z={}", getNuclearA(), getNuclearZ()) : "n/a"));
}
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline PDGCode NuclearParticleInterface<InnerParticleInterface,
+ StackIteratorInterface>::getPDG() const {
+ return (isNucleus() ? PDGCode(1000000000 + getNuclearZ() * 10000 + getNuclearA() * 10)
+ : super_type::getPDG());
+ }
+
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline void
+ NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::setMomentum(
+ MomentumVector const& v) {
+ HEPMomentumType const P = v.getNorm();
+ if (P == 0_eV) {
+ 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().setKineticEnergy(
+ super_type::getIndex(),
+ sqrt(square(this->getMass()) + square(P)) - this->getMass());
+ super_type::getStackData().setDirection(super_type::getIndex(), v / P);
+ }
+ }
+
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline MomentumVector NuclearParticleInterface<
+ InnerParticleInterface, StackIteratorInterface>::getMomentum() const {
+ 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,
@@ -121,6 +252,22 @@ namespace corsika::nuclear_stack {
return super_type::getCharge();
}
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline HEPEnergyType NuclearParticleInterface<
+ 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,
typename StackIteratorInterface>
inline int16_t NuclearParticleInterface<
diff --git a/corsika/detail/stack/VectorStack.inl b/corsika/detail/stack/VectorStack.inl
index a89c85ed951595ed585cee02649dc8e698733aac..40fb4e61c759e256af9d6fcac6283c6eee44fdb2 100644
--- a/corsika/detail/stack/VectorStack.inl
+++ b/corsika/detail/stack/VectorStack.inl
@@ -24,67 +24,106 @@ 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));
- this->setMomentum(std::get<2>(v));
+ 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 / sqrt(P2));
+ }
+ this->setPosition(std::get<2>(v));
+ this->setTime(std::get<3>(v));
+ }
+
+ template <typename StackIteratorInterface>
+ inline void ParticleInterface<StackIteratorInterface>::setParticleData(
+ ParticleInterface<StackIteratorInterface> const&,
+ particle_data_momentum_type const& v) {
+ this->setPID(std::get<0>(v));
+ 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 / sqrt(P2));
+ }
+ this->setPosition(std::get<2>(v));
+ this->setTime(std::get<3>(v));
+ }
+
+ template <typename StackIteratorInterface>
+ inline void ParticleInterface<StackIteratorInterface>::setParticleData(
+ particle_data_type const& v) {
+ this->setPID(std::get<0>(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 void ParticleInterface<StackIteratorInterface>::setParticleData(
- ParticleInterface<StackIteratorInterface> const&,
- std::tuple<Code, HEPEnergyType, MomentumVector, 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->setMomentum(std::get<2>(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();
- momentum_.clear();
+ dataEkin_.clear();
+ direction_.clear();
position_.clear();
time_.clear();
}
inline void VectorStackImpl::copy(size_t i1, size_t i2) {
dataPID_[i2] = dataPID_[i1];
- dataE_[i2] = dataE_[i1];
- momentum_[i2] = momentum_[i1];
+ dataEkin_[i2] = dataEkin_[i1];
+ direction_[i2] = direction_[i1];
position_[i2] = position_[i1];
time_[i2] = time_[i1];
}
inline void VectorStackImpl::swap(size_t i1, size_t i2) {
std::swap(dataPID_[i2], dataPID_[i1]);
- std::swap(dataE_[i2], dataE_[i1]);
- std::swap(momentum_[i2], momentum_[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]);
}
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();
- momentum_.push_back(
- MomentumVector(dummyCS, {0 * electronvolt, 0 * electronvolt, 0 * electronvolt}));
+ direction_.push_back(DirectionVector(dummyCS, {0, 0, 0}));
position_.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
time_.push_back(0 * second);
}
inline void VectorStackImpl::decrementSize() {
- if (dataE_.size() > 0) {
+ if (dataEkin_.size() > 0) {
dataPID_.pop_back();
- dataE_.pop_back();
- momentum_.pop_back();
+ dataEkin_.pop_back();
+ direction_.pop_back();
position_.pop_back();
time_.pop_back();
}
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index 7be6db2887dea865c53346caa83f3b6df5edfd84..6ec30d7487c847c7606aebb091e406cd809b8c6f 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -71,19 +71,19 @@ namespace corsika {
int16_t constexpr get_charge_number(Code const); //!< electric charge in units of e
ElectricChargeType constexpr get_charge(Code const); //!< electric charge
HEPMassType constexpr get_mass(Code const); //!< mass
- HEPEnergyType constexpr get_energy_threshold(
- Code const); //!< get energy threshold below which the particle is discarded, by
- //!< default set to particle mass
- void constexpr set_energy_threshold(
- Code const, HEPEnergyType const); //!< set energy threshold below which the particle
- //!< is discarded
-
- inline void set_energy_threshold(std::pair<Code const, HEPEnergyType const> p) {
- set_energy_threshold(p.first, p.second);
+ HEPEnergyType constexpr get_kinetic_energy_threshold(
+ Code const); //!< get kinetic energy threshold below which the particle is
+ //!< discarded, by default set to zero
+ void constexpr set_kinetic_energy_threshold(
+ Code const, HEPEnergyType const); //!< set kinetic energy threshold below which the
+ //!< particle is discarded
+
+ inline void set_kinetic_energy_threshold(std::pair<Code const, HEPEnergyType const> p) {
+ set_kinetic_energy_threshold(p.first, p.second);
}
- inline void set_energy_thresholds(
+ inline void set_kinetic_energy_thresholds(
std::unordered_map<Code const, HEPEnergyType const> const& eCuts) {
- for (auto v : eCuts) set_energy_threshold(v);
+ for (auto v : eCuts) set_kinetic_energy_threshold(v);
}
//! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
diff --git a/corsika/framework/stack/Stack.hpp b/corsika/framework/stack/Stack.hpp
index 2ff173fb3ccabcd8d5944ce4dc7452cd96314f80..caaf9d065b709fc0161c87f41bd2ac3661f8b2ac 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 76815746af5ab15ce5aa61fc48c2f2249f987b32..02ef7e5423c99e64248d4a5135f2d8262c221eb8 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/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index 11d8049d53e2b2ca817feaa16b7bfaea438ee70a..7857c691f53024263ab2dd8e284cab8777a9e317 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -21,15 +21,18 @@
namespace corsika {
/**
simple ParticleCut process. Goes through the secondaries of an interaction and
- removes particles according to their energy. Particles with a time delay of more than
- 10ms are removed as well. Invisible particles (neutrinos) can be removed if selected.
+ removes particles according to their kinetic energy. Particles with a time delay of
+ more than 10ms are removed as well. Invisible particles (neutrinos) can be removed if
+ selected. The threshold value is set to 0 by default but in principle can be configured
+ for each particle. Special constructors for cuts by the following groups are
+ implemented: (electrons,positrons), photons, hadrons and muons.
**/
class ParticleCut : public SecondariesProcess<ParticleCut>,
public ContinuousProcess<ParticleCut> {
public:
/**
- * particle cut with energy thresholds for electrons, photons,
+ * particle cut with kinetic energy thresholds for electrons, photons,
* hadrons (including nuclei with energy per nucleon) and muons
* invisible particles (neutrinos) can be cut or not
**/
@@ -63,14 +66,30 @@ namespace corsika {
void showResults(); // LCOV_EXCL_LINE
void reset();
- HEPEnergyType getElectronECut() const { return get_energy_threshold(Code::Electron); }
- HEPEnergyType getPhotonECut() const { return get_energy_threshold(Code::Photon); }
- HEPEnergyType getMuonECut() const { return get_energy_threshold(Code::MuPlus); }
- HEPEnergyType getHadronECut() const { return get_energy_threshold(Code::Proton); }
- HEPEnergyType getInvEnergy() const { return inv_energy_; }
- HEPEnergyType getCutEnergy() const { return energy_; }
- HEPEnergyType getEmEnergy() const { return em_energy_; }
+ HEPEnergyType getElectronKineticECut() const {
+ return get_kinetic_energy_threshold(Code::Electron);
+ }
+ HEPEnergyType getPhotonKineticECut() const {
+ return get_kinetic_energy_threshold(Code::Photon);
+ }
+ HEPEnergyType getMuonKineticECut() const {
+ return get_kinetic_energy_threshold(Code::MuPlus);
+ }
+ HEPEnergyType getHadronKineticECut() const {
+ return get_kinetic_energy_threshold(Code::Proton);
+ }
+ //! returns total energy of particles that were removed by cut for invisible particles
+ HEPEnergyType getInvEnergy() const { return energy_invcut_; }
+ //! returns total energy of particles that were removed by cut in time
+ HEPEnergyType getTimeCutEnergy() const { return energy_timecut_; }
+ //! returns total energy of particles that were removed by cut in kinetic energy
+ HEPEnergyType getCutEnergy() const { return energy_cut_; }
+ //! returns total energy of particles that were removed by cut for electromagnetic
+ //! particles
+ HEPEnergyType getEmEnergy() const { return energy_emcut_; }
+ //! returns number of electromagnetic particles
unsigned int getNumberEmParticles() const { return em_count_; }
+ //! returns number of invisible particles
unsigned int getNumberInvParticles() const { return inv_count_; }
private:
@@ -86,10 +105,11 @@ namespace corsika {
private:
bool doCutEm_;
bool doCutInv_;
- HEPEnergyType energy_ = 0 * electronvolt;
- HEPEnergyType em_energy_ = 0 * electronvolt;
+ HEPEnergyType energy_cut_ = 0 * electronvolt;
+ HEPEnergyType energy_timecut_ = 0 * electronvolt;
+ HEPEnergyType energy_emcut_ = 0 * electronvolt;
+ HEPEnergyType energy_invcut_ = 0 * electronvolt;
unsigned int em_count_ = 0;
- HEPEnergyType inv_energy_ = 0 * electronvolt;
unsigned int inv_count_ = 0;
};
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index 19e9a44c32aa401147d57585c89a7c1ab3ec6f4e..6c928656642959d93bf07a469a98e04637884bc8 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -47,20 +47,74 @@ namespace corsika::nuclear_stack {
typedef InnerParticleInterface<StackIteratorInterface> super_type;
public:
- typedef std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType>
- particle_data_type;
-
- typedef std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType,
+ typedef std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType,
unsigned short, unsigned short>
- altenative_particle_data_type;
+ nuclear_particle_data_type;
+
+ typedef std::tuple<Code, MomentumVector, Point, TimeType, unsigned short,
+ unsigned short>
+ nuclear_particle_data_momentum_type;
+
+ /**
+ *
+ * @param v which is a tuple containing: PID, kinetic Energy, DirectionVector,
+ * Position, Time
+ */
+ void setParticleData(typename super_type::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);
+
+ /**
+ *
+ * @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);
- void setParticleData(particle_data_type const& v);
+ /**
+ *
+ * @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(altenative_particle_data_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(super_type& p, particle_data_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, altenative_particle_data_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);
+
+ /**
+ *
+ * @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;
@@ -89,13 +143,45 @@ namespace corsika::nuclear_stack {
/// @}
/**
- * Overwrite normal getParticleMass function with nuclear version
+ * Overwrite normal getPDG function with nuclear version
+ */
+ PDGCode getPDG() const;
+
+ /**
+ * 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
*/
HEPMassType getMass() const;
+
/**
* Overwrite normal getParticleCharge function with nuclear version
*/
ElectricChargeType getCharge() const;
+
+ /**
+ * Overwrite normal getEnergy function with nuclear version
+ */
+ HEPEnergyType getEnergy() const;
+
/**
* Overwirte normal getChargeNumber function with nuclear version
**/
@@ -172,6 +258,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 e5a765781ae44f38b38da7d4135a6363c697f5d3..223d651f90a6a3d37c5d45e362fa8702e21e9c54 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -33,48 +33,118 @@ 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;
- void setParticleData(
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v);
+ std::string asString() const;
- void setParticleData(
- ParticleInterface<StackIteratorInterface> const&,
- std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v);
+ /**
+ * 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(particle_data_type const& v);
+
+ /**
+ * 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& p,
+ particle_data_type 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);
+
+ /**
+ * 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);
- /// individual setters
+ ///! 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.
+ */
void setMomentum(MomentumVector const& v) {
- super_type::getStackData().setMomentum(super_type::getIndex(), v);
+ HEPMomentumType const P = v.getNorm();
+ if (P == 0_eV) {
+ 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().setKineticEnergy(
+ super_type::getIndex(),
+ sqrt(square(getMass()) + square(P)) - this->getMass());
+ super_type::getStackData().setDirection(super_type::getIndex(), v / P);
+ }
+ }
+ //! Set direction
+ void setDirection(DirectionVector const& v) {
+ super_type::getStackData().setDirection(super_type::getIndex(), v);
}
+ //! Set position
void setPosition(Point const& v) {
super_type::getStackData().setPosition(super_type::getIndex(), v);
}
+ //! Set time
void setTime(TimeType const& v) {
super_type::getStackData().setTime(super_type::getIndex(), v);
}
- /// individual getters
+ //! Get corsika::Code
Code getPID() const {
return super_type::getStackData().getPID(super_type::getIndex());
}
- HEPEnergyType getEnergy() const {
- return super_type::getStackData().getEnergy(super_type::getIndex());
+ //! Get PDG code
+ PDGCode getPDG() const { return get_PDG(getPID()); }
+ //! Get kinetic energy
+ HEPEnergyType getKineticEnergy() const {
+ return super_type::getStackData().getKineticEnergy(super_type::getIndex());
}
- MomentumVector getMomentum() const {
- return super_type::getStackData().getMomentum(super_type::getIndex());
+ //! Get direction
+ DirectionVector getDirection() const {
+ return super_type::getStackData().getDirection(super_type::getIndex());
}
+ //! Get position
Point getPosition() const {
return super_type::getStackData().getPosition(super_type::getIndex());
}
+ //! Get time
TimeType getTime() const {
return super_type::getStackData().getTime(super_type::getIndex());
}
@@ -83,18 +153,25 @@ namespace corsika {
*
* @{
*/
- DirectionVector getDirection() const {
- return this->getMomentum() / this->getEnergy();
- }
-
+ //! 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()); }
+ //! Get kinetic energy
+ HEPEnergyType getEnergy() const { return this->getKineticEnergy() + this->getMass(); }
+
+ //! Get charge number
int16_t getChargeNumber() const { return get_charge_number(this->getPID()); }
///@}
};
@@ -102,16 +179,18 @@ namespace corsika {
/**
* Memory implementation of the most simple (stupid) particle stack object.
*
+ * @note if we ever want to have off-shell particles, we need to
+ * add momentum as HEPMomentumType, and a lot of care.
*/
class VectorStackImpl {
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<MomentumVector> momentum_vector_type;
+ typedef std::vector<DirectionVector> direction_vector_type;
VectorStackImpl() = default;
@@ -131,27 +210,17 @@ 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 setMomentum(size_t i, MomentumVector const& v) { momentum_[i] = v; }
+ 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]; }
-
- MomentumVector getMomentum(size_t i) const { return momentum_[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]; }
- HEPEnergyType getDataE(size_t i) const { return dataE_[i]; }
-
- void setDataE(size_t i, HEPEnergyType const& dataE) { dataE_[i] = dataE; }
-
- Code getDataPid(size_t i) const { return dataPID_[i]; }
-
- void setDataPid(size_t i, Code const& dataPid) { dataPID_[i] = dataPid; }
/**
* Function to copy particle at location i2 in stack to i1
*/
@@ -168,8 +237,8 @@ namespace corsika {
private:
/// the actual memory to store particle data
code_vector_type dataPID_;
- energy_vector_type dataE_;
- momentum_vector_type momentum_;
+ 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 ebb202597bc8f7d09e88b815e89e5f55cf5a3225..3ce8df868981a77af9d4df509f414264a2c3138d 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 f7758db0a3163521274ee5b0b33e21cc478c4405..2852428d2c074b0afaa22fcdc58ec5eeeaa0825c 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 a02c914e4ad1d554160d36745e5d00a2ee6d981e..71be6f4f9edc102cbc5d5d190955abe757d07191 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 606d72b41cf10445d62bb3d769e1465e90a2068b..12400723bf9c29550f8c424c9b55d117c8b7a82c 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 f1602f5ce24a7359d8e1de500369f70d32177b5f..25603ad8af4c3f797d5cde7ae47d0190cf960a29 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 30187c520a4d0badbef4126f5e192ffd5906702f..41beda4f9c5c6f220411f675fe8e3d8c76e7eda3 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 d419efbb86b9cbfb3bc5a2b8db3cc8adb3a6e2bb..4c9e43efc1883e933433aa83240896442491bb9b 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 b6f346546176ed5cd9e0b9a5ffd960d9893fb438..21354adcfc31d8fec2ec28381a7aeb4ceb39e39e 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -6,6 +6,8 @@
* the license.
*/
+#define TRACE
+
/* clang-format off */
// InteractionCounter used boost/histogram, which
// fails if boost/type_traits have been included before. Thus, we have
@@ -318,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));
}
}
@@ -392,7 +394,7 @@ int main(int argc, char** argv) {
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
// register the observation plane with the output
- output.add("obsplane", observationLevel);
+ output.add("particles", observationLevel);
auto sequence = make_sequence( // emCascadeCounted,
stackInspect, hadronSequence, reset_particle_mass, decaySequence,
diff --git a/src/framework/core/pdxml_reader.py b/src/framework/core/pdxml_reader.py
index d8809a004f5e21690aa37acbea1712da714a2c87..13a79f33482d77e2e037217bf9ebd0e05a0b45bb 100755
--- a/src/framework/core/pdxml_reader.py
+++ b/src/framework/core/pdxml_reader.py
@@ -358,11 +358,10 @@ def gen_properties(particle_db):
mass=p['mass'], name=p['name'])
string += "};\n\n"
- # particle threshold table, initially set to the particle mass
- string += "static std::array<corsika::units::si::HEPEnergyType, size> thresholds = {\n"
- for p in particle_db.values():
- string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(
- mass=p['mass'], name=p['name'])
+ # particle threshold table, initially set to 0
+ string += "static std::array<corsika::units::si::HEPEnergyType, size> thresholds = {\n"
+ for k in particle_db:
+ string += " 0 * corsika::units::si::electronvolt, // {name:s}\n".format( name = k)
string += "};\n\n"
# PDG code table
diff --git a/tests/common/SetupTestStack.hpp b/tests/common/SetupTestStack.hpp
index 3d5322d4a79f5a4be0d397348daff9bbf5530c87..d41f1aa2726f65a7eeb7a31b3132e38ad51e1be8 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 44d277bae1bf03a7fdfe9fa12eef3ae2c2197f95..1d5cb6be27ba186e7ebaccc7872d2dbbabfc9cd8 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/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp
index 809ae5e1391071b948a6307600cabb4bfc201136..79073518bd6e8ea302f8a1414b73cacc03ccfdd2 100644
--- a/tests/framework/testGeometry.cpp
+++ b/tests/framework/testGeometry.cpp
@@ -30,7 +30,6 @@ double constexpr absMargin = 1.0e-8;
TEST_CASE("Geometry CoordinateSystems") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
@@ -66,6 +65,9 @@ TEST_CASE("Geometry CoordinateSystems") {
Point p2(translatedCS, {0_m, 0_m, 0_m});
CHECK(((p2 - p1).getComponents() - translationVector).getNorm().magnitude() ==
Approx(0).margin(absMargin));
+ CHECK(p2.getX(rootCS) == 0_m);
+ CHECK(p2.getY(rootCS) == 4_m);
+ CHECK(p2.getZ(rootCS) == 0_m);
}
SECTION("multiple translations") {
@@ -332,6 +334,8 @@ TEST_CASE("Distance between points") {
CHECK(distance(p5, p6) / 1_m == Approx(1));
}
+TEST_CASE("Geometry Tree") {}
+
TEST_CASE("Path") {
// define a known CS
CoordinateSystemPtr root = get_root_CoordinateSystem();
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index c35f8bca97982e2d2388ab122cfad22269dd4960..322c553c7a45358a61e3e389fafb888a15b7d59d 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -86,12 +86,12 @@ TEST_CASE("ParticleProperties", "[Particles]") {
}
SECTION("Energy threshold") {
- //! by default energy thresholds are set to particle mass
- CHECK(get_energy_threshold(Electron::code) / Electron::mass == Approx(1));
+ //! by default energy thresholds are set to zero
+ CHECK(get_kinetic_energy_threshold(Electron::code) == 0_GeV);
- set_energy_threshold(Electron::code, 10_GeV);
- CHECK_FALSE(get_energy_threshold(Code::Electron) == 1_GeV);
- CHECK(get_energy_threshold(Code::Electron) == 10_GeV);
+ set_kinetic_energy_threshold(Electron::code, 10_GeV);
+ CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == 1_GeV);
+ CHECK(get_kinetic_energy_threshold(Code::Electron) == 10_GeV);
}
SECTION("Particle groups: electromagnetic") {
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index 02105eb9e343f8cf51fc034a3d708b6f88ec2b76..1335974cdd70104ced3ff639b8afb3610fffcec3 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -62,10 +62,11 @@ struct DummyView {
int globalCount = 0; // simple counter
-int checkDecay = 0; // use this as a bit field
-int checkInteract = 0; // use this as a bit field
-int checkSec = 0; // use this as a bit field
-int checkCont = 0; // use this as a bit field
+int checkDecay = 0; // use this as a bit field
+int checkInteract = 0; // use this as a bit field
+int checkSec = 0; // use this as a bit field
+int checkCont = 0; // use this as a bit field
+int checkSecondaries = 0; // use this as a bit field
class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
public:
@@ -332,6 +333,22 @@ private:
int count_ = 0;
};
+class Secondaries1 : public SecondariesProcess<Secondaries1> {
+public:
+ template <typename TView>
+ void doSecondaries(TView const&) {
+ checkSecondaries |= 1;
+ }
+};
+
+class Secondaries2 : public SecondariesProcess<Secondaries2> {
+public:
+ template <typename TView>
+ void doSecondaries(TView const&) {
+ checkSecondaries |= 2;
+ }
+};
+
class Boundary1 : public BoundaryCrossingProcess<Boundary1> {
public:
Boundary1(double const v = 1.0)
@@ -580,7 +597,6 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
/**
* In this example switching is done only by "data_[0]>0", where
@@ -596,8 +612,11 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
auto cp2 = ContinuousProcess2(0, 2_m);
auto cp3 = ContinuousProcess3(0, 3_m);
- auto sequence1 = make_sequence(Process1(0), cp2, Decay1(0), Boundary1(1.0));
- auto sequence2 = make_sequence(cp3, Process2(0), Boundary1(-1.0), Decay2(0));
+ auto sec1 = Secondaries1();
+ auto sec2 = Secondaries2();
+
+ auto sequence1 = make_sequence(Process1(0), cp2, Decay1(0), sec1, Boundary1(1.0));
+ auto sequence2 = make_sequence(cp3, Process2(0), Boundary1(-1.0), Decay2(0), sec2);
auto sequence3 = make_sequence(cp1, Process3(0),
SwitchProcessSequence(select1, sequence1, sequence2));
@@ -650,7 +669,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkInteract = 0;
checkSec = 0;
checkCont = 0;
- particle.data_[0] = 100; // data positive
+ particle.data_[0] = 100; // data positive --> sequence1
sequence3.doContinuous(particle, track, ContinuousProcessIndex(1));
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
@@ -661,7 +680,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkInteract = 0;
checkSec = 0;
checkCont = 0;
- particle.data_[0] = -100; // data negative
+ particle.data_[0] = -100; // data negative --> sequence2
sequence3.doContinuous(particle, track, ContinuousProcessIndex(1));
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
@@ -676,7 +695,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkInteract = 0;
checkSec = 0;
checkCont = 0;
- particle.data_[0] = 100; // data positive
+ particle.data_[0] = 100; // data positive --> sequence1
sequence3.selectInteraction(view, lambda_select);
sequence3.selectDecay(view, time_select);
CHECK(checkInteract == 0b100); // this is Process3
@@ -692,7 +711,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkInteract = 0;
checkSec = 0;
checkCont = 0;
- particle.data_[0] = -100; // data negative
+ particle.data_[0] = -100; // data negative --> sequence2
sequence3.selectInteraction(view, lambda_select);
sequence3.selectDecay(view, time_select);
CHECK(checkInteract == 0b010); // this is Process2
@@ -704,7 +723,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkInteract = 0;
checkSec = 0;
checkCont = 0;
- particle.data_[0] = -100; // data negative
+ particle.data_[0] = -100; // data negative --> sequence2
sequence3.doSecondaries(view);
Stack1 stack(0);
sequence3.doStack(stack);
@@ -720,7 +739,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkInteract = 0;
checkSec = 0;
checkCont = 0;
- particle.data_[0] = -100; // data positive
+ particle.data_[0] = -100; // data negative --> sequence1
sequence4.selectInteraction(view, lambda_select);
sequence4.doSecondaries(view);
sequence4.selectDecay(view, time_select);
@@ -741,6 +760,22 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
CHECK(checkDecay == 0);
}
+ SECTION("Check SecondariesProcesses in SwitchProcessSequence") {
+
+ DummyData particle;
+ DummyView view(particle);
+
+ checkSecondaries = 0;
+ particle.data_[0] = 100; // data positive --> sequence1
+ sequence3.doSecondaries(view);
+ CHECK(checkSecondaries == 1);
+
+ checkSecondaries = 0;
+ particle.data_[0] = -100; // data positive --> sequence1
+ sequence3.doSecondaries(view);
+ CHECK(checkSecondaries == 2);
+ }
+
SECTION("Check ContinuousProcesses in SwitchProcessSequence") {
DummyData particle;
diff --git a/tests/framework/testSecondaryView.cpp b/tests/framework/testSecondaryView.cpp
index fa113f88b6e111e9837d9ef96d28a432024974dd..9ca07a13d1e567398456207874cbd47093902956 100644
--- a/tests/framework/testSecondaryView.cpp
+++ b/tests/framework/testSecondaryView.cpp
@@ -76,6 +76,7 @@ TEST_CASE("SecondaryStack", "[stack]") {
CHECK(view.getSize() == 0);
CHECK(view.getEntries() == 0);
CHECK(view.isEmpty());
+ CHECK_THROWS(view.erase(view.begin())); // nothing to delete yet
{
auto proj = view.getProjectile();
@@ -112,6 +113,7 @@ TEST_CASE("SecondaryStack", "[stack]") {
auto pDel = view.getNextParticle();
view.erase(pDel);
+ CHECK_THROWS(view.erase(pDel)); // already erased
CHECK(view.getSize() == 2);
CHECK(stack.getSize() == 4);
diff --git a/tests/framework/testSolver.cpp b/tests/framework/testSolver.cpp
index 52f1fa13fc29f59d7b90e4471b76bb51997b0613..423b73b8b60bcb4c90c4f3e2ff9adbf732b2e724 100644
--- a/tests/framework/testSolver.cpp
+++ b/tests/framework/testSolver.cpp
@@ -53,7 +53,6 @@ TEST_CASE("Solver") {
{7e8, -1e-7}};
for (auto v : vals) {
-
{
double a = v.first;
double b = v.second;
@@ -74,7 +73,9 @@ TEST_CASE("Solver") {
}
CHECK(solve_linear_real(0, 55.).size() == 0);
- }
+ CHECK(solve_linear(0, 55.).size() == 0);
+
+ } // linear
SECTION("quadratic") {
@@ -129,7 +130,7 @@ TEST_CASE("Solver") {
}
}
}
- }
+ } // quadratic
SECTION("cubic") {
diff --git a/tests/framework/testStackInterface.cpp b/tests/framework/testStackInterface.cpp
index 4171647af5a547baeb55e6a2e3c3701eaeb77b29..1bb77d166ebd71e86efa9b6901144247f8b0eec8 100644
--- a/tests/framework/testStackInterface.cpp
+++ b/tests/framework/testStackInterface.cpp
@@ -71,6 +71,9 @@ TEST_CASE("Stack", "[Stack]") {
SECTION("delete from stack") {
StackTest stack;
+
+ CHECK_THROWS(stack.erase(stack.begin())); // nothing to delete
+
CHECK(stack.getSize() == 0);
StackTest::stack_iterator_type p =
stack.addParticle(std::tuple{0.}); // valid way to access particle data
@@ -95,7 +98,8 @@ TEST_CASE("Stack", "[Stack]") {
CHECK(stack.getEntries() == 3);
CHECK(!stack.isEmpty());
- p.erase(); // mark for deletion: size=3, entries=2
+ p.erase(); // mark for deletion: size=3, entries=2
+ CHECK_THROWS(p.erase()); // already deleted
CHECK(stack.getSize() == 3);
CHECK(stack.getEntries() == 2);
CHECK(!stack.isEmpty());
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index 3eab2cd6b9e5cdce3efc87ae48a9281145fbcb72..9a5cfab0e81d99f816abfc7da29c31e4e8ebe378 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -38,10 +38,43 @@ CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+TEST_CASE("VolumeTree") {
+ logging::set_level(logging::level::info);
+ Environment<IEmpty> env;
+ auto& universe = *(env.getUniverse());
+ auto world = Environment<IEmpty>::createNode<Sphere>(Point{gCS, 0_m, 0_m, 0_m}, 150_km);
+ // volume cut partly by "world"
+ auto vol1 =
+ Environment<IEmpty>::createNode<Sphere>(Point{gCS, 0_m, 0_m, 140_km}, 20_km);
+ // partly overlap with "vol1"
+ auto vol2 = Environment<IEmpty>::createNode<Sphere>(Point{gCS, 0_m, 0_m, 120_km}, 5_km);
+ vol1->excludeOverlapWith(vol2);
+ world->addChild(std::move(vol1));
+ world->addChild(std::move(vol2));
+
+ // in world
+ CHECK(dynamic_cast<Sphere const&>(
+ world->getContainingNode(Point(gCS, 0_m, 0_m, 0_m))->getVolume())
+ .getRadius() == 150_km);
+ // in vol1
+ CHECK(dynamic_cast<Sphere const&>(
+ world->getContainingNode(Point(gCS, 0_m, 0_m, 149_km))->getVolume())
+ .getRadius() == 20_km);
+ // outside world, in universe
+ CHECK(world->getContainingNode(Point(gCS, 0_m, 151_km, 0_m)) == nullptr);
+ // in vol2
+ CHECK(dynamic_cast<Sphere const&>(
+ world->getContainingNode(Point(gCS, 0_m, 0_km, 119_km))->getVolume())
+ .getRadius() == 5_km);
+ CHECK(dynamic_cast<Sphere const&>(
+ world->getContainingNode(Point(gCS, 0_m, 0_km, 121_km))->getVolume())
+ .getRadius() == 5_km);
+ universe.addChild(std::move(world));
+}
+
TEST_CASE("HomogeneousMedium") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
std::vector<float>{1.f});
@@ -54,7 +87,6 @@ TEST_CASE("HomogeneousMedium") {
TEST_CASE("FlatExponential") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
std::vector<float>{1.f});
@@ -117,7 +149,6 @@ TEST_CASE("FlatExponential") {
TEST_CASE("SlidingPlanarExponential") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
std::vector<float>{1.f});
@@ -178,7 +209,6 @@ struct Exponential {
TEST_CASE("InhomogeneousMedium") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
Vector direction(gCS, QuantityVector<dimensionless_d>(1, 0, 0));
@@ -227,7 +257,6 @@ TEST_CASE("InhomogeneousMedium") {
TEST_CASE("LayeredSphericalAtmosphereBuilder") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
LayeredSphericalAtmosphereBuilder builder =
make_layered_spherical_atmosphere_builder<>::create(gOrigin,
@@ -267,7 +296,6 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// setup our interface types
using ModelInterface = IMagneticFieldModel<IMediumModel>;
diff --git a/tests/media/testMedium.cpp b/tests/media/testMedium.cpp
index edcba45fc6fe01a3b5554878fcca735627b7a9d7..709185d18524b7d5536657ffbd005bbafde9cfb9 100644
--- a/tests/media/testMedium.cpp
+++ b/tests/media/testMedium.cpp
@@ -27,7 +27,6 @@ using namespace corsika;
TEST_CASE("MediumProperties") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// test access of medium properties via enum and class types
@@ -45,7 +44,6 @@ TEST_CASE("MediumProperties") {
TEST_CASE("MediumPropertyModel w/ Homogeneous") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CoordinateSystemPtr gCS = get_root_CoordinateSystem();
diff --git a/tests/media/testShowerAxis.cpp b/tests/media/testShowerAxis.cpp
index 19d7635aedebd4ff901396af86cc4651ba94b3e2..9378d9b41fcab19f9bf3a1b1547c3f9cdb59e6b4 100644
--- a/tests/media/testShowerAxis.cpp
+++ b/tests/media/testShowerAxis.cpp
@@ -48,7 +48,6 @@ auto setupEnvironment(Code vTargetCode) {
TEST_CASE("Homogeneous Density") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto [env, csPtr, nodePtr] = setupEnvironment(Code::Nitrogen);
auto const& cs = *csPtr;
@@ -62,7 +61,7 @@ TEST_CASE("Homogeneous Density") {
Point const injectionPos = showerCore + Vector<dimensionless_d>{cs, {0, 0, 1}} * t;
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos), *env,
- true, // -> do not throw exceptions
+ true, // -> throw exceptions
20}; // -> number of bins
CHECK(showerAxis.getSteplength() == 500_m);
@@ -83,4 +82,11 @@ TEST_CASE("Homogeneous Density") {
CHECK_THROWS(showerAxis.getX(-1_m));
CHECK_THROWS(showerAxis.getX((injectionPos - showerCore).getNorm() + 1_m));
+
+ ShowerAxis const showerAxisNoThrow{injectionPos, (showerCore - injectionPos), *env,
+ false, // -> do not throw exceptions
+ 20}; // -> number of bins
+ CHECK(showerAxisNoThrow.getX(-1_m) == showerAxis.getMinimumX());
+ CHECK(showerAxisNoThrow.getX((injectionPos - showerCore).getNorm() + 1_m) ==
+ showerAxis.getMaximumX());
}
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index d72de5b5314c8af20d58e8b04adead3994e5db95..5c0b9e860a75c4f673d6cabdb156293702afd4aa 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -8,7 +8,7 @@ set (test_modules_sources
testPythia8.cpp
testUrQMD.cpp
testCONEX.cpp
- testOnShellCheck.cpp
+ # testOnShellCheck.cpp
testParticleCut.cpp
testSibyll.cpp
)
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index ab47527d80a254975e2585d6a8a356ebf7a1a420..94ec3f27f194160547d35fbab2df13007b6c93b2 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -50,13 +50,14 @@ TEST_CASE("ParticleCut", "processes") {
SECTION("cut on particle type: inv") {
+ // particle cut with 20GeV threshold for all, also cut invisible
ParticleCut cut(20_GeV, false, true);
- CHECK(cut.getHadronECut() == 20_GeV);
+ 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.
@@ -66,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);
@@ -82,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.
@@ -94,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.
@@ -119,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);
@@ -139,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);
@@ -174,10 +174,11 @@ TEST_CASE("ParticleCut", "processes") {
SECTION("cut low energy: reset thresholds of arbitrary set of particles") {
ParticleCut cut({{Code::Electron, 5_MeV}, {Code::Positron, 50_MeV}}, false, true);
- CHECK(get_energy_threshold(Code::Electron) != get_energy_threshold(Code::Positron));
- CHECK_FALSE(get_energy_threshold(Code::Electron) == Electron::mass);
+ CHECK(get_kinetic_energy_threshold(Code::Electron) !=
+ get_kinetic_energy_threshold(Code::Positron));
+ CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == Electron::mass);
// test default values still correct
- CHECK(get_energy_threshold(Code::Proton) == 5_GeV);
+ CHECK(get_kinetic_energy_threshold(Code::Proton) == 5_GeV);
}
SECTION("cut on time") {
@@ -185,27 +186,27 @@ 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.
// only this way the secondary view is populated
auto projectile = view.getProjectile();
// add secondaries, all with energies above the threshold
- // only cut is by species
+ // 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);
CHECK(view.getEntries() == 0);
- CHECK(cut.getCutEnergy() / 1_GeV == 11000);
- cut.reset();
+ CHECK(cut.getTimeCutEnergy() == 11 * Eabove);
CHECK(cut.getCutEnergy() == 0_GeV);
+ cut.reset();
+ CHECK(cut.getTimeCutEnergy() == 0_GeV);
}
setup::Trajectory const track = setup::testing::make_track<setup::Trajectory>(
@@ -215,13 +216,14 @@ TEST_CASE("ParticleCut", "processes") {
SECTION("cut on DoContinous, just invisibles") {
ParticleCut cut(20_GeV, false, true);
- CHECK(cut.getHadronECut() == 20_GeV);
+ CHECK(cut.getHadronKineticECut() == 20_GeV);
// 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 fdbbc82e120ea4f675295143e1050acb76e94e73..846c734877badaaa8cbb508f30da1716ed24b78d 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -101,9 +101,6 @@ TEST_CASE("Pythia8Interface", "modules") {
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;
@@ -161,8 +158,6 @@ TEST_CASE("Pythia8Interface", "modules") {
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/testSibyll.cpp b/tests/modules/testSibyll.cpp
index dc5249a7218e2e2dea96b7a23f791833ef4e8c9a..85b512740c1471e9a83d40947330ddbb7694dcc0 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -203,8 +203,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
*/
CHECK(pSum.getComponents(cs).getX() / P0 == Approx(1).margin(0.05));
- CHECK(pSum.getComponents(cs).getY() / 1_GeV == Approx(0).margin(1e-4));
- CHECK(pSum.getComponents(cs).getZ() / 1_GeV == Approx(0).margin(1e-4));
+ CHECK(pSum.getComponents(cs).getY() / 1_GeV == Approx(0).margin(1e-3));
+ CHECK(pSum.getComponents(cs).getZ() / 1_GeV == Approx(0).margin(1e-3));
CHECK((pSum - plab).getNorm() / 1_GeV ==
Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
@@ -212,8 +212,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
CHECK(length / 1_g * 1_cm * 1_cm == Approx(88.7).margin(0.1));
// CHECK(view.getEntries() == 9); //! \todo: this was 20 before refactory-2020: check
- // also sibyll not stable wrt. to compiler
- // changes
+ // "also sibyll not stable wrt. to compiler
+ // changes"
}
SECTION("NuclearInteractionInterface") {
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
index 0b2a3ab7b9a489df955f00814ce371dfc4a026e9..9905583ed579e73cd059330805dd49f146e0c009 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 6ee9300c3ab4fff3c74b834bb4cbfca8df828c0f..56e11820979f3b7b0bf9ee2dd72a473118cf7d16 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 1e33f6c0570230bed7fc31baa13f0f2de1cb0c91..fa9201023ee13e6d4fcc00460ecd67175f0b950d 100644
--- a/tests/stack/testNuclearStackExtension.cpp
+++ b/tests/stack/testNuclearStackExtension.cpp
@@ -17,10 +17,9 @@ using namespace corsika;
#include <iostream>
using namespace std;
-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();
@@ -28,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);
}
@@ -40,37 +39,41 @@ 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);
}
SECTION("read nucleus") {
+ 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, 10, 9));
+ 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() == 1.5_GeV);
CHECK(pout.getTime() == 100_s);
CHECK(pout.getNuclearA() == 10);
CHECK(pout.getNuclearZ() == 9);
@@ -78,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());
@@ -93,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));
}
}
@@ -114,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));
}
}
@@ -221,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 9767b6c847beff9097e873aed8a390ccb62a45fd..7fa0903b456d6bcd26ebf72a728beeb88746eba3 100644
--- a/tests/stack/testVectorStack.cpp
+++ b/tests/stack/testVectorStack.cpp
@@ -17,26 +17,26 @@
using namespace corsika;
using namespace std;
-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.getEnergy() == 1.5_GeV + Electron::mass);
+ CHECK(pout.getKineticEnergy() == 1.5_GeV);
CHECK(pout.getTime() == 100_s);
}
@@ -45,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);