diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index 40edbeb439b616278990eb82b0613ef86e97b0bd..73c2e6f89cc1df311bc312a386b87945e5a1d6e5 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -14,7 +14,7 @@
namespace corsika {
- inline HEPEnergyType constexpr get_kinetic_energy_threshold(Code const code) {
+ inline HEPEnergyType constexpr calculate_kinetic_energy_threshold(Code const code) {
if (is_nucleus(code)) return particle::detail::threshold_nuclei;
return particle::detail::thresholds[static_cast<CodeIntType>(code)];
}
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index 44bcdf4244df288e1c3baa386fedc56840969997..0358bbb5a158ccf76765fac902280c3832618e7a 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -108,9 +108,9 @@ namespace corsika {
if (is_nucleus(pid)) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / get_nucleus_A(pid);
- return (ElabNuc < get_kinetic_energy_threshold(pid));
+ return (ElabNuc < calculate_kinetic_energy_threshold(pid));
} else {
- return (energyLab < get_kinetic_energy_threshold(pid));
+ return (energyLab < calculate_kinetic_energy_threshold(pid));
}
}
@@ -183,7 +183,7 @@ namespace corsika {
inline void ParticleCut::printThresholds() {
for (auto p : get_all_particles()) {
- auto const Eth = get_kinetic_energy_threshold(p);
+ auto const Eth = calculate_kinetic_energy_threshold(p);
CORSIKA_LOG_INFO("kinetic energy threshold for particle {} is {} GeV", p,
Eth / 1_GeV);
}
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 537a3d85fbf347303d00a1f519fa92ed8317e4d7..c840de798880c3c94a939aff13ffc7b02a32f28a 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -181,7 +181,7 @@ namespace corsika {
auto const energy = vParticle.getKineticEnergy();
auto const energy_lim =
std::max(energy * 0.9, // either 10% relative loss max., or
- get_kinetic_energy_threshold(
+ calculate_kinetic_energy_threshold(
vParticle.getPID()) // energy thresholds globally defined
// for individual particles
* 0.99999 // need to go slightly below global e-cut to assure
diff --git a/corsika/detail/modules/epos/InteractionModel.inl b/corsika/detail/modules/epos/InteractionModel.inl
index 35ed15a884c6e608056fcb5b562ddfd1d0c5e50e..f9e8925b7b2da6d94744a6d2afcad1c4fefdfa6a 100644
--- a/corsika/detail/modules/epos/InteractionModel.inl
+++ b/corsika/detail/modules/epos/InteractionModel.inl
@@ -430,11 +430,6 @@ namespace corsika::epos {
MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType E_final = 0_GeV;
- // position and time of interaction, not used in QgsjetII
- auto const& projectile = view.getProjectile();
- Point const& pOrig = projectile.getPosition();
- TimeType const tOrig = projectile.getTime();
-
// secondaries
EposStack es;
CORSIKA_LOGGER_DEBUG(logger_, "number of entries on Epos stack: {}", es.getSize());
@@ -449,12 +444,14 @@ namespace corsika::epos {
EposCode const eposId = psec.getPID();
Code const pid = epos::convertFromEpos(eposId);
+ HEPEnergyType const mass = get_mass(pid);
+ HEPEnergyType const Ekin = sqrt(p3output.getSquaredNorm() + mass * mass) - mass;
CORSIKA_LOGGER_TRACE(logger_,
" id= {}"
" p= {}",
pid, p3output.getComponents() / 1_GeV);
- auto pnew = view.addSecondary(std::make_tuple(pid, p3output, pOrig, tOrig));
+ auto pnew = view.addSecondary(std::make_tuple(pid, Ekin, p3output.normalized()));
P_final += pnew.getMomentum();
E_final += pnew.getEnergy();
}
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 30f0b10a1925d78fcdd7ad178b12accf31ffb358..818398ef889402c84f1b065686a0c2bea07870e4 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -29,7 +29,7 @@ namespace corsika::proposal {
// take some minutes if you have to build the tables and cannot read the
// from disk
auto const emCut =
- get_kinetic_energy_threshold(code) +
+ calculate_kinetic_energy_threshold(code) +
get_mass(code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
@@ -129,7 +129,7 @@ namespace corsika::proposal {
auto const energy = vP.getEnergy();
auto const energy_lim =
std::max(energy * 0.9, // either 10% relative loss max., or
- get_kinetic_energy_threshold(
+ calculate_kinetic_energy_threshold(
code) // energy thresholds globally defined for individual particles
* 0.9999 // need to go slightly below global e-cut to assure removal
// in ParticleCut. This does not matter since at cut-time
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index e81f9984b7b880ddd208556b087df0d0302b3d59..af6d603f6874113c4430c157085557956ca19581 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -33,7 +33,7 @@ namespace corsika::proposal {
// take some minutes if you have to build the tables and cannot read the
// from disk
auto const emCut =
- get_kinetic_energy_threshold(code) +
+ calculate_kinetic_energy_threshold(code) +
get_mass(code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
@@ -93,12 +93,10 @@ namespace corsika::proposal {
for (auto& s : sec) {
auto E = s.energy * 1_MeV;
auto vecProposal = s.direction;
- auto vec = QuantityVector(vecProposal.GetX() * E, vecProposal.GetY() * E,
- vecProposal.GetZ() * E);
- auto p = MomentumVector(labCS, vec);
+ auto dir = DirectionVector(
+ labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
- view.addSecondary(
- std::make_tuple(sec_code, p, projectile.getPosition(), projectile.getTime()));
+ view.addSecondary(std::make_tuple(sec_code, E - get_mass(sec_code), dir));
}
}
return ProcessReturn::Ok;
diff --git a/corsika/detail/modules/pythia8/Decay.inl b/corsika/detail/modules/pythia8/Decay.inl
index 81306df4509492448ba4425b1fd8e874f307ad41..3da0ebdbf3869072b97b40198cb53a10f0ec16b8 100644
--- a/corsika/detail/modules/pythia8/Decay.inl
+++ b/corsika/detail/modules/pythia8/Decay.inl
@@ -170,9 +170,6 @@ namespace corsika::pythia8 {
auto projectile = view.getProjectile();
- auto const& decayPoint = projectile.getPosition();
- auto const t0 = projectile.getTime();
-
auto const& labMomentum = projectile.getMomentum();
[[maybe_unused]] CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
@@ -230,14 +227,17 @@ namespace corsika::pythia8 {
{event[i].px() * 1_GeV, event[i].py() * 1_GeV, event[i].pz() * 1_GeV});
FourVector const fourMomRest{Erest, pRest};
auto const fourMomLab = boost.fromCoM(fourMomRest);
+ auto const p3 = fourMomLab.getSpaceLikeComponents();
+
+ HEPEnergyType const mass = get_mass(pyId);
+ HEPEnergyType const Ekin = sqrt(p3.getSquaredNorm() + mass * mass) - mass;
CORSIKA_LOG_TRACE(
"particle: id={} momentum={} energy={} ", pyId,
fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV,
fourMomLab.getTimeLikeComponent());
- view.addSecondary(
- std::make_tuple(pyId, fourMomLab.getSpaceLikeComponents(), decayPoint, t0));
+ view.addSecondary(std::make_tuple(pyId, Ekin, p3.normalized()));
}
// set particle stable
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index fe55d587412537f86dff4ca4f0ca0f9dd5ef99ab..2144380e7480f3325375a45052c0822880c9e036 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -247,8 +247,12 @@ namespace corsika::pythia8 {
MomentumVector const pyPlab(labCS,
{p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+ HEPEnergyType const mass = get_mass(pyId);
+ HEPEnergyType const Ekin = sqrt(pyPlab.getSquaredNorm() + mass * mass) - mass;
+
// add to corsika stack
- auto pnew = projectile.addSecondary(std::make_tuple(pyId, pyPlab, pOrig, tOrig));
+ auto pnew =
+ projectile.addSecondary(std::make_tuple(pyId, Ekin, pyPlab.normalized()));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
diff --git a/corsika/detail/modules/qgsjetII/InteractionModel.inl b/corsika/detail/modules/qgsjetII/InteractionModel.inl
index 8203faa49d1c9c6c3fe16de453f01dcde138e3bb..13ad2d0c814ff994a47c4bd00ec8d9632d85442d 100644
--- a/corsika/detail/modules/qgsjetII/InteractionModel.inl
+++ b/corsika/detail/modules/qgsjetII/InteractionModel.inl
@@ -7,14 +7,14 @@
*/
#include <corsika/modules/qgsjetII/InteractionModel.hpp>
-
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/geometry/Point.hpp>
-
#include <corsika/modules/qgsjetII/ParticleConversion.hpp>
#include <corsika/modules/qgsjetII/QGSJetIIFragmentsStack.hpp>
#include <corsika/modules/qgsjetII/QGSJetIIStack.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <sstream>
@@ -73,9 +73,8 @@ namespace corsika::qgsjetII {
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
auto const sqrtS = sqrt(sqrtS2);
if (!isValid(projectileId, targetId, sqrtS)) { return CrossSectionType::zero(); }
- HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
- static_pow<2>(get_mass(targetId))) /
- (2 * get_mass(targetId));
+ HEPEnergyType const Elab =
+ calculate_lab_energy(sqrtS2, get_mass(projectileId), get_mass(targetId));
int const iBeam = static_cast<QgsjetIIXSClassIntType>(
corsika::qgsjetII::getQgsjetIIXSCode(projectileId));
@@ -111,9 +110,8 @@ namespace corsika::qgsjetII {
!isValid(projectileId, targetId, sqrtS)) {
throw std::runtime_error("invalid target/projectile/energy combination.");
}
- HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
- static_pow<2>(get_mass(targetId))) /
- (2 * get_mass(targetId));
+ HEPEnergyType const Elab =
+ calculate_lab_energy(sqrtS2, get_mass(projectileId), get_mass(targetId));
int beamA = 0;
if (is_nucleus(projectileId)) { beamA = get_nucleus_A(projectileId); }
@@ -178,11 +176,6 @@ namespace corsika::qgsjetII {
// internal QGSJetII system
COMBoost boostInternal({Elab, pLab}, get_mass(targetId));
- // position and time of interaction, not used in QgsjetII
- auto const projectile = view.getProjectile();
- Point const pOrig = projectile.getPosition();
- TimeType const tOrig = projectile.getTime();
-
// fragments
QGSJetIIFragmentsStack qfs;
for (auto& fragm : qfs) {
@@ -208,11 +201,16 @@ namespace corsika::qgsjetII {
auto p3output = P4output.getSpaceLikeComponents();
p3output.rebase(originalCS); // transform back into standard lab frame
+ HEPEnergyType const mass = get_mass(idFragm);
+ HEPEnergyType const Ekin = sqrt(p3output.getSquaredNorm() + mass * mass) - mass;
+
CORSIKA_LOG_DEBUG(
"secondary fragment> id= {}"
" p={}",
idFragm, p3output.getComponents());
- auto pnew = view.addSecondary(std::make_tuple(idFragm, p3output, pOrig, tOrig));
+
+ auto pnew =
+ view.addSecondary(std::make_tuple(idFragm, Ekin, p3output.normalized()));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
@@ -241,8 +239,7 @@ namespace corsika::qgsjetII {
MomentumVector momentum{
csPrime,
{0.0_GeV, 0.0_GeV,
- sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
- (projectileEnergyLabPerNucleon * A - nucleusMass))}};
+ calculate_momentum(projectileEnergyLabPerNucleon * A, nucleusMass)}};
// this is not "CoM" here, but rather the system defined by projectile+target,
// which in Cascade-mode is already lab
@@ -252,24 +249,27 @@ namespace corsika::qgsjetII {
auto p3output = P4output.getSpaceLikeComponents();
p3output.rebase(originalCS); // transform back into standard lab frame
+ Code const idFragm = get_nucleus_code(A, Z);
+ HEPEnergyType const mass = get_mass(idFragm);
+ HEPEnergyType const Ekin = calculate_kinetic_energy(p3output.getNorm(), mass);
+
CORSIKA_LOG_DEBUG(
"secondary fragment> id={}"
" p={}"
" A={}"
" Z={}",
- get_nucleus_code(A, Z), p3output.getComponents(), A, Z);
+ idFragm, p3output.getComponents(), A, Z);
- auto pnew = view.addSecondary(
- std::make_tuple(get_nucleus_code(A, Z), p3output, pOrig, tOrig));
+ auto pnew =
+ view.addSecondary(std::make_tuple(idFragm, Ekin, p3output.normalized()));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
- }
+ } // namespace corsika::qgsjetII
// secondaries
QGSJetIIStack qs;
for (auto& psec : qs) {
-
auto momentum = psec.getMomentum(csPrime);
// this is not "CoM" here, but rather the system defined by projectile+target,
// which in Cascade-mode is already lab
@@ -278,11 +278,12 @@ namespace corsika::qgsjetII {
auto p3output = P4output.getSpaceLikeComponents();
p3output.rebase(originalCS); // transform back into standard lab frame
- CORSIKA_LOG_DEBUG("secondary> id= {}, p= {}",
- corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
- p3output.getComponents());
- auto pnew = view.addSecondary(std::make_tuple(
- corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), p3output, pOrig, tOrig));
+ Code const pid = corsika::qgsjetII::convertFromQgsjetII(psec.getPID());
+ HEPEnergyType const mass = get_mass(pid);
+ HEPEnergyType const Ekin = calculate_kinetic_energy(p3output.getNorm(), mass);
+
+ CORSIKA_LOG_DEBUG("secondary> id= {}, p= {}", pid, p3output.getComponents());
+ auto pnew = view.addSecondary(std::make_tuple(pid, Ekin, p3output.normalized()));
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 ab0114182b41f5db321f19ec6d91458cc9b5000d..b804feb534aedea966eddd1c5f378c7c4ba0329d 100644
--- a/corsika/detail/modules/sibyll/Decay.inl
+++ b/corsika/detail/modules/sibyll/Decay.inl
@@ -173,9 +173,7 @@ namespace corsika::sibyll {
throw std::runtime_error("STOP! Sibyll not configured to execute this decay!");
count_++;
- // remember position
- Point const& decayPoint = projectile.getPosition();
- TimeType const t0 = projectile.getTime();
+
// switch on decay for this particle
setUnstable(pCode);
printDecayConfig(pCode);
@@ -222,8 +220,11 @@ namespace corsika::sibyll {
CORSIKA_LOG_TRACE("Sibyll::Decay: i={} id={} p={} GeV", i, pid, components / 1_GeV);
- projectile.addSecondary(
- std::make_tuple(pid, MomentumVector(rootCS, components), decayPoint, t0));
+ auto const p3 = MomentumVector(rootCS, components);
+ HEPEnergyType const mass = get_mass(pid);
+ HEPEnergyType const Ekin = sqrt(p3.getSquaredNorm() + mass * mass) - mass;
+
+ projectile.addSecondary(std::make_tuple(pid, Ekin, p3.normalized()));
}
}
diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
index 99f078673c2d80df3c6251bd6d22eed25a077c4f..dd00e7606aabf439c677bf2dbaf74b34e1ddab82 100644
--- a/corsika/detail/modules/sibyll/InteractionModel.inl
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -138,11 +138,6 @@ namespace corsika::sibyll {
// add particles from sibyll to stack
- // position and time of interaction, not used in Sibyll
- auto const& projectile = secondaries.parent();
- Point const& pOrig = projectile.getPosition();
- TimeType const tOrig = projectile.getTime(); // no time in sibyll
-
// link to sibyll stack
SibStack ss;
@@ -163,9 +158,13 @@ namespace corsika::sibyll {
auto const P4lab = boost.fromCoM(FourVector{eCoM, pCoM});
auto const p3lab = P4lab.getSpaceLikeComponents();
+ Code const pid = corsika::sibyll::convertFromSibyll(psib.getPID());
+ HEPEnergyType const mass = get_mass(pid);
+ HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+
// add to corsika stack
- auto pnew = secondaries.addSecondary(std::make_tuple(
- corsika::sibyll::convertFromSibyll(psib.getPID()), p3lab, pOrig, tOrig));
+ auto pnew =
+ secondaries.addSecondary(std::make_tuple(pid, Ekin, p3lab.normalized()));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
diff --git a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
index 4a564c7e96bcc64e059dfec58b43eb3f5b8ed7ac..cc12434f2da112145a9bbe16960157391a21e05b 100644
--- a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
@@ -292,14 +292,6 @@ namespace corsika::sibyll {
}
// (LCOV_EXCL_STOP)
- // position and time of interaction, not used in NUCLIB
- auto const& projectile = view.parent();
- // position and time of interaction, not used in NUCLI
- Point const& pOrig = projectile.getPosition();
- TimeType const delay = projectile.getTime();
-
- CORSIKA_LOG_DEBUG("Interaction: position of interaction: {}, {} ns",
- pOrig.getCoordinates(), delay / 1_ns);
CORSIKA_LOG_DEBUG("number of fragments: {}", nFragments);
CORSIKA_LOG_DEBUG("adding nuclear fragments to particle stack..");
// put nuclear fragments on corsika stack
@@ -324,8 +316,11 @@ namespace corsika::sibyll {
// CORSIKA 7 way
// spectators inherit momentum from original projectile
auto const p3lab = p3NucleonLab * nuclA;
+
+ HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+
CORSIKA_LOG_DEBUG("fragment momentum {}", p3lab.getComponents() / 1_GeV);
- view.addSecondary(std::make_tuple(specCode, p3lab, pOrig, delay));
+ view.addSecondary(std::make_tuple(specCode, Ekin, p3lab.normalized()));
}
// add elastic nucleons to corsika stack
@@ -340,7 +335,11 @@ namespace corsika::sibyll {
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
auto const p3lab = p3NucleonLab;
- view.addSecondary(std::make_tuple(elaNucCode, p3lab, pOrig, delay));
+
+ HEPEnergyType const mass = get_mass(elaNucCode);
+ HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+
+ view.addSecondary(std::make_tuple(elaNucCode, Ekin, p3lab.normalized()));
}
// add inelastic interactions
@@ -348,12 +347,18 @@ namespace corsika::sibyll {
for (int j = 0; j < nInelNucleons; ++j) {
// TODO: sample neutron or proton
auto const pCode = Code::Proton;
+ HEPEnergyType const mass = get_mass(pCode);
+ HEPEnergyType const Ekin = sqrt(p3NucleonLab.getSquaredNorm() + mass * mass) - mass;
+
// temporarily add to stack, will be removed after interaction in DoInteraction
CORSIKA_LOG_DEBUG("inelastic interaction no. {}", j);
typename TSecondaryView::inner_stack_value_type nucleonStack;
- auto inelasticNucleon =
- nucleonStack.addParticle(std::make_tuple(pCode, p3NucleonLab, pOrig, delay));
+ Point const pDummy(boost.getOriginalCS(), {0_m, 0_m, 0_m});
+ TimeType const tDummy = 0_ns;
+ auto inelasticNucleon = nucleonStack.addParticle(
+ std::make_tuple(pCode, Ekin, p3NucleonLab.normalized(), pDummy, tDummy));
inelasticNucleon.setNode(view.getProjectile().getNode());
+
// create inelastic interaction for each nucleon
CORSIKA_LOG_TRACE("calling HadronicInteraction...");
// create new StackView for each of the nucleons
@@ -362,8 +367,12 @@ namespace corsika::sibyll {
hadronicInteraction_.doInteraction(nucleon_secondaries, pCode, targetId, nucleonP4,
targetP4);
for (const auto& pSec : nucleon_secondaries) {
- view.addSecondary(std::make_tuple(pSec.getPID(), pSec.getMomentum(),
- pSec.getPosition(), pSec.getTime()));
+
+ auto const p3lab = pSec.getMomentum();
+ Code const pid = pSec.getPID();
+ HEPEnergyType const mass = get_mass(pid);
+ HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+ view.addSecondary(std::make_tuple(pid, Ekin, p3lab.normalized()));
}
}
}
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index 1bf0f0d59a42428df231717f2df869f72c30bd45..9e070f1c62df8860f9bbb4d89df26f9c87728e09 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -259,10 +259,6 @@ namespace corsika::urqmd {
iflb_; // flag for retrying interaction in case of empty event, 0 means retry
::urqmd::urqmd_(iflb);
- auto projectile = view.getProjectile();
- auto const& projectilePosition = projectile.getPosition();
- auto const projectileTime = projectile.getTime();
-
// now retrieve secondaries from UrQMD
COMBoost const boost(projectileP4, targetP4);
auto const& originalCS = boost.getOriginalCS();
@@ -283,8 +279,15 @@ namespace corsika::urqmd {
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(" {} {} {} ", i, code, momentum.getComponents());
- view.addSecondary(
- std::make_tuple(code, momentum, projectilePosition, projectileTime));
+ HEPEnergyType const mass = get_mass(code);
+ HEPEnergyType Ekin = sqrt(momentum.getSquaredNorm() + mass * mass) - mass;
+ if (Ekin <= 0_GeV) {
+ CORSIKA_LOG_WARN("Negative kinetic energy {} {}. Skipping.", code, Ekin);
+ view.addSecondary(
+ std::make_tuple(code, 0_eV, DirectionVector{originalCS, {0, 0, 0}}));
+ } else {
+ view.addSecondary(std::make_tuple(code, Ekin, momentum.normalized()));
+ }
}
CORSIKA_LOG_DEBUG("UrQMD generated {} secondaries!", ::urqmd::sys_.npart);
}
diff --git a/corsika/detail/stack/VectorStack.inl b/corsika/detail/stack/VectorStack.inl
index 57bc094a438dc6f7383baa39215434099777690b..d777f24b3b9844dc6e4c78461fc8de37ee870fdf 100644
--- a/corsika/detail/stack/VectorStack.inl
+++ b/corsika/detail/stack/VectorStack.inl
@@ -24,58 +24,34 @@ namespace corsika {
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
- particle_data_momentum_type const& v) {
+ particle_data_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));
+ 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& parent,
- 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) + parent.getTime()); // parent time is added
- }
-
- template <typename StackIteratorInterface>
- inline void ParticleInterface<StackIteratorInterface>::setParticleData(
- particle_data_type const& v) {
+ secondary_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));
+ this->setPosition(parent.getPosition()); // position
+ this->setTime(parent.getTime()); // parent time
}
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
ParticleInterface<StackIteratorInterface> const& parent,
- particle_data_type const& v) {
+ secondary_extended_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) + parent.getTime()); // parent time is added
+ this->setPosition(parent.getPosition() + std::get<3>(v)); // + position
+ this->setTime(parent.getTime() + std::get<4>(v)); // + parent time
}
template <typename StackIteratorInterface>
diff --git a/corsika/framework/core/EnergyMomentumOperations.hpp b/corsika/framework/core/EnergyMomentumOperations.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c7abc940eeff4c5897deccb553ae4177e686eb53
--- /dev/null
+++ b/corsika/framework/core/EnergyMomentumOperations.hpp
@@ -0,0 +1,120 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+
+/**
+ * \file EnergyMomentumOperations.hpp
+ *
+ * Relativistic energy momentum calculations.
+ */
+
+namespace corsika {
+
+ namespace detail {
+ using HEPEnergyTypeSqr = decltype(1_GeV * 1_GeV);
+ }
+
+ /**
+ * \f[m^2=E_{tot}^2-p^2\f]
+ *
+ * @param E total energy.
+ * @param p particle momentum.
+ * @return HEPEnergyType-squared
+ */
+ auto constexpr calculate_mass_sqr(HEPEnergyType const E, HEPMomentumType const p) {
+ return (E - p) * (E + p);
+ }
+
+ /**
+ * \f[m=sqrt(E_{tot}^2-p^2) \f]
+ *
+ * @param E total energy.
+ * @param p particle momentum.
+ * @return HEPEnergyType
+ */
+ HEPEnergyType constexpr calculate_mass(HEPEnergyType const E, HEPMomentumType const p) {
+ return sqrt(calculate_mass_sqr(E, p));
+ }
+
+ /**
+ * \f[p^2=E_{tot}^2-m^2\f]
+ *
+ * @param E total energy.
+ * @param m particle mass.
+ * @return HEPEnergyType-squared
+ */
+ auto constexpr calculate_momentum_sqr(HEPEnergyType const E, HEPMassType const m) {
+ return (E - m) * (E + m);
+ }
+
+ /**
+ * \f[p=sqrt(E_{tot}^2-m^2) \f]
+ *
+ * @param E total energy.
+ * @param m particle mass.
+ * @return HEPEnergyType
+ */
+ HEPEnergyType constexpr calculate_momentum(HEPEnergyType const E, HEPMassType const m) {
+ return sqrt(calculate_momentum_sqr(E, m));
+ }
+
+ /**
+ * \f[E_{tot}^2=p^2+m^2\f]
+ *
+ * @param p momentum.
+ * @param m particle mass.
+ * @return HEPEnergyType-squared
+ */
+ auto constexpr calculate_total_energy_sqr(HEPMomentumType const p,
+ HEPMassType const m) {
+ return p * p + m * m;
+ }
+
+ /**
+ * \f[E_{tot}=sqrt(p^2+m^2)\f]
+ *
+ * @param p momentum.
+ * @param m particle mass.
+ * @return HEPEnergyType
+ */
+ HEPEnergyType constexpr calculate_total_energy(HEPMomentumType const p,
+ HEPMassType const m) {
+ return sqrt(calculate_total_energy_sqr(p, m));
+ }
+
+ /**
+ * \f[E_{kin}=sqrt(p^2+m^2) - m \f]
+ *
+ * @param p momentum.
+ * @param m particle mass.
+ * @return HEPEnergyType
+ */
+ HEPEnergyType constexpr calculate_kinetic_energy(HEPMomentumType const p,
+ HEPMassType const m) {
+ return calculate_total_energy(p, m) - m;
+ }
+
+ /**
+ * \f[E_{lab}=(\sqrt{s}^2 - m_{proj}^2 - m_{targ}^2) / (2 m_{targ}) \f]
+ *
+ * @param p momentum.
+ * @param m particle mass.
+ * @return HEPEnergyType
+ */
+ HEPEnergyType constexpr calculate_lab_energy(detail::HEPEnergyTypeSqr sqrtS_sqr,
+ HEPMassType const m_proj,
+ HEPMassType const m_targ) {
+ return (sqrtS_sqr - static_pow<2>(m_proj) - static_pow<2>(m_targ)) / (2 * m_targ);
+ }
+
+} // namespace corsika
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index 8f6391fee658264c0c830f823285cfc7efa02021..1b344954194c056005874eff5a28d0f15a1468ab 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -92,7 +92,7 @@ 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_kinetic_energy_threshold(
+ HEPEnergyType constexpr calculate_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(
diff --git a/corsika/framework/units/quantity.hpp b/corsika/framework/units/quantity.hpp
index 7b5bd35877b755fe1d751965975c148d42a8f20e..1ebd98015c757bfe451b838f89eea5995d22688e 100644
--- a/corsika/framework/units/quantity.hpp
+++ b/corsika/framework/units/quantity.hpp
@@ -454,7 +454,7 @@ namespace phys {
// powers and roots
template <int N, typename D, typename X>
- friend detail::Power<D, N, X> nth_power(quantity<D, X> const& x);
+ friend constexpr detail::Power<D, N, X> nth_power(quantity<D, X> const& x);
template <typename D, typename X>
friend constexpr detail::Power<D, 2, X> square(quantity<D, X> const& x);
@@ -463,13 +463,13 @@ namespace phys {
friend constexpr detail::Power<D, 3, X> cube(quantity<D, X> const& x);
template <int N, typename D, typename X>
- friend detail::Root<D, N, X> nth_root(quantity<D, X> const& x);
+ friend detail::Root<D, N, X> constexpr nth_root(quantity<D, X> const& x);
template <typename D, typename X>
- friend detail::Root<D, 2, X> sqrt(quantity<D, X> const& x);
+ friend detail::Root<D, 2, X> constexpr sqrt(quantity<D, X> const& x);
template <typename D, typename X>
- friend detail::Root<D, 3, X> cbrt(quantity<D, X> const& x);
+ friend detail::Root<D, 3, X> constexpr cbrt(quantity<D, X> const& x);
// comparison
@@ -629,7 +629,7 @@ namespace phys {
/// absolute value.
template <typename D, typename X>
- constexpr quantity<D, X> abs(quantity<D, X> const& x) {
+ quantity<D, X> constexpr abs(quantity<D, X> const& x) {
return quantity<D, X>(std::abs(x.m_value));
}
@@ -638,7 +638,7 @@ namespace phys {
/// N-th power.
template <int N, typename D, typename X>
- detail::Power<D, N, X> nth_power(quantity<D, X> const& x) {
+ detail::Power<D, N, X> constexpr nth_power(quantity<D, X> const& x) {
return detail::Power<D, N, X>(std::pow(x.m_value, X(N)));
}
@@ -663,7 +663,7 @@ namespace phys {
/// n-th root.
template <int N, typename D, typename X>
- detail::Root<D, N, X> nth_root(quantity<D, X> const& x) {
+ detail::Root<D, N, X> constexpr nth_root(quantity<D, X> const& x) {
static_assert(detail::root<D, N, X>::all_even_multiples,
"root result dimensions must be integral");
@@ -675,7 +675,7 @@ namespace phys {
/// square root.
template <typename D, typename X>
- detail::Root<D, 2, X> sqrt(quantity<D, X> const& x) {
+ detail::Root<D, 2, X> constexpr sqrt(quantity<D, X> const& x) {
static_assert(detail::root<D, 2, X>::all_even_multiples,
"root result dimensions must be integral");
@@ -685,7 +685,7 @@ namespace phys {
/// cubic root.
template <typename D, typename X>
- detail::Root<D, 3, X> cbrt(quantity<D, X> const& x) {
+ detail::Root<D, 3, X> constexpr cbrt(quantity<D, X> const& x) {
static_assert(detail::root<D, 3, X>::all_even_multiples,
"root result dimensions must be integral");
diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp
index bb26ddfa8004de4d593eb8c2376060bc677f62ca..159186e1df41de1ed27bd1c8a8d6b3d34ea1a0a1 100644
--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -67,16 +67,16 @@ namespace corsika {
void reset();
HEPEnergyType getElectronKineticECut() const {
- return get_kinetic_energy_threshold(Code::Electron);
+ return calculate_kinetic_energy_threshold(Code::Electron);
}
HEPEnergyType getPhotonKineticECut() const {
- return get_kinetic_energy_threshold(Code::Photon);
+ return calculate_kinetic_energy_threshold(Code::Photon);
}
HEPEnergyType getMuonKineticECut() const {
- return get_kinetic_energy_threshold(Code::MuPlus);
+ return calculate_kinetic_energy_threshold(Code::MuPlus);
}
HEPEnergyType getHadronKineticECut() const {
- return get_kinetic_energy_threshold(Code::Proton);
+ return calculate_kinetic_energy_threshold(Code::Proton);
}
//! returns total energy of particles that were removed by cut for invisible particles
HEPEnergyType getInvEnergy() const { return energy_invcut_; }
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index 22873790c22c1e478a0c6883008ed95e54a0562b..32dc0af1706197730dcf4d9fd911d80a2276b83c 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -33,19 +33,37 @@ namespace corsika {
typedef ParticleBase<TStackIterator> super_type;
public:
+ /**
+ * particle data information content.
+ *
+ * PID, Ekin, direction, position, time.
+ */
typedef std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType>
particle_data_type;
- typedef std::tuple<Code, MomentumVector, Point, TimeType> particle_data_momentum_type;
+ /**
+ * secondary particle data information content.
+ *
+ * PID, Ekin, direction.
+ */
+ typedef std::tuple<Code, HEPEnergyType, DirectionVector> secondary_data_type;
+
+ /**
+ * secondary particle data information content with position and time update.
+ *
+ * PID, Ekin, direction, delta-Position, delta-Time.
+ */
+ typedef std::tuple<Code, HEPEnergyType, DirectionVector, Vector<length_d>, TimeType>
+ secondary_extended_data_type;
std::string asString() const;
/**
* Set data of new particle.
*
- * @param v tuple containing: PID, Momentum Vector, Position, Time
+ * @param v tuple containing of type particle_data_type
*
- * MomentumVector is only used to determine the DirectionVector, the normalization
+ * MomentumVector is only used to determine the DirectionVector, the normalization
* is lost.
*/
void setParticleData(particle_data_type const& v);
@@ -54,29 +72,19 @@ namespace corsika {
* Set data of new particle.
*
* @param parent parent particle
- * @param v tuple containing: PID, Momentum Vector, Position, Time
- *
- * MomentumVector is only used to determine the DirectionVector, the normalization
- * is lost.
+ * @param v tuple containing of type secondary_data_type.
*/
void setParticleData(ParticleInterface<TStackIterator> const& parent,
- 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);
+ secondary_data_type const& v);
/**
* Set data of new particle.
*
* @param parent parent particle
- * @param v tuple containing: PID, kinetic Energy, Direction Vector, Position, Time
+ * @param v tuple containing of type secondary_extended_data_type.
*/
void setParticleData(ParticleInterface<TStackIterator> const& parent,
- particle_data_momentum_type const& v);
+ secondary_extended_data_type const& v);
///! Set particle corsika::Code
void setPID(Code const id) {
diff --git a/documentation/media.rst b/documentation/media.rst
index 6a24954c02ee6158c54015bd1be12592c4625d1f..7d7589788efd86b1f6d5070169f88ac7d037da24 100644
--- a/documentation/media.rst
+++ b/documentation/media.rst
@@ -1,5 +1,5 @@
-Media Properties
-================
+Medium Properties
+=================
.. toctree::
media_classes
diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp
index 8b28769038504483db2cc3b516d62f162d2dbd10..a11da70c6559c8d593a75dcaad2fe006b0c67927 100644
--- a/examples/boundary_example.cpp
+++ b/examples/boundary_example.cpp
@@ -6,13 +6,14 @@
* the license.
*/
-#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/output/OutputManager.hpp>
@@ -160,7 +161,9 @@ 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, plab, pos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), pos, 0_ns));
}
// define air shower object, run simulation
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index cf01b33ca499cf8d1df3d61bcc142f0aee921cdf..2f0d254472a477890d382ffe34ab56d6c6dabad8 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -6,14 +6,15 @@
* the license.
*/
-#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Logging.hpp>
#include <corsika/output/OutputManager.hpp>
@@ -125,7 +126,9 @@ 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, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
}
// setup processes, decays and interactions
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index fda0a3af6062b3011ee073c30e7fd2438349ae91..9934211b74786c9c9b0f7cdf3c70d5cf15c1005c 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -6,13 +6,15 @@
* the license.
*/
-#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Logging.hpp>
#include <corsika/output/OutputManager.hpp>
@@ -109,7 +111,9 @@ 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, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
}
// setup processes, decays and interactions
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index 7df6dcbde54c0d40190da9a377c166d9472c588b..6cf3051ad33a9ef51aa7bd8c91b8ed6a7e3dda81 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -12,18 +12,19 @@
// to include it first...
#include <corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/output/NoOutput.hpp>
@@ -395,7 +396,9 @@ int main(int argc, char** argv) {
stack.clear();
// add the desired particle to the stack
- stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
// if we want to fix the first location of the shower
if (force_interaction) {
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index b01d7b608b85d661d7fdce6004a1ef9efcd3cd4b..fd6191d84a5f24495b711b9ad794745368b25967 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -6,17 +6,20 @@
* the license.
*/
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
-#include <corsika/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/output/OutputManager.hpp>
@@ -130,7 +133,9 @@ int main(int argc, char** argv) {
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
- stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), 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 0ff10af164ff51fc3db6cde5bbb7324dea410d24..2de965757d4e0c3c0e0df83af5ab2cf47aa9055b 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -12,18 +12,19 @@
// to include it first...
#include <corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/output/OutputManager.hpp>
@@ -183,7 +184,9 @@ int main(int argc, char** argv) {
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
- stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ Code::Proton, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
<< std::endl;
diff --git a/examples/mars.cpp b/examples/mars.cpp
index afa8c71b70b8e2526785312acc38e1c7b684484c..12a4039bb89c983b50408a52369d4b0d5ef4af52 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -14,16 +14,19 @@
/* clang-format on */
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Cascade.hpp>
+
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/output/NoOutput.hpp>
@@ -418,7 +421,9 @@ int main(int argc, char** argv) {
stack.clear();
// add the desired particle to the stack
- stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
// run the shower
EAS.run();
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index e971a279d27edd940098b228d335efdca9d27594..8ca094f6e5c9f26b436236e21f71c12a7f5c2b7a 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -15,6 +15,7 @@
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <fstream>
@@ -75,7 +76,9 @@ int main() {
momentumComponents(theta / 180. * constants::pi, phi / 180. * constants::pi, P0);
auto plab = MomentumVector(rootCS, {px, py, pz});
- stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), 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 6bb97291aa2b2d7d3106fcf68a6b887aaf99c7d9..f38b1b2904d5eb777e0449571e892af58d82ccea 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -14,18 +14,19 @@
// to include it first...
#include <corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/output/NoOutput.hpp>
@@ -251,7 +252,9 @@ int main(int argc, char** argv) {
setup::Stack stack;
stack.clear();
- stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+ stack.addParticle(std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns));
BetheBlochPDG emContinuous(showerAxis);
diff --git a/tests/common/SetupTestStack.hpp b/tests/common/SetupTestStack.hpp
index a6141e582262b77e4062eb020bfdc7e3a19adece..476c99a54bcd2a5487b9d4cdee591e88607a3409 100644
--- a/tests/common/SetupTestStack.hpp
+++ b/tests/common/SetupTestStack.hpp
@@ -27,9 +27,6 @@ namespace corsika::setup::testing {
*
* standard stack setup for unit tests.
*
- *
- *
- *
* \return a tuple with element 0 being a Stack object filled with
* one particle, and element 1 the StackView on it.
*/
@@ -43,9 +40,11 @@ namespace corsika::setup::testing {
Point const origin(cs, {0_m, 0_m, 0_m});
MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
+ HEPMassType const mass = get_mass(vProjectileType);
+ HEPEnergyType const Ekin = sqrt(vMomentum * vMomentum + mass * mass) - mass;
- auto particle =
- stack->addParticle(std::make_tuple(vProjectileType, pLab, origin, 0_ns));
+ auto particle = stack->addParticle(
+ std::make_tuple(vProjectileType, Ekin, pLab.normalized(), 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 18732cb1f5faaa646fd9333a89faa1e9e2861973..624dc041575d1304c4e7bdcb012810e7ca3cc478 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -98,10 +98,10 @@ public:
++calls_;
auto vP = view.getProjectile();
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(), Ekin / 2, vP.getMomentum().normalized(),
- vP.getPosition(), vP.getTime()));
+ vP.addSecondary(
+ std::make_tuple(vP.getPID(), Ekin / 2, vP.getMomentum().normalized()));
+ vP.addSecondary(
+ std::make_tuple(vP.getPID(), Ekin / 2, vP.getMomentum().normalized()));
}
int getCalls() const { return calls_; }
@@ -162,11 +162,10 @@ TEST_CASE("Cascade", "[Cascade]") {
auto sequence = make_sequence(nullModel, stackInspect, split, cut);
TestCascadeStack stack;
stack.clear();
- stack.addParticle(std::make_tuple(
- 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));
+ stack.addParticle(std::make_tuple(Code::Electron,
+ E0 - get_mass(Code::Electron), // Ekin
+ DirectionVector(rootCS, {0, 0, -1}),
+ Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
DummyTracking tracking;
DummyOutputManager output;
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index f12faf8039a3baa7031416e49bfe5e9bc4c8d467..3ec511b627161dec1a09389f8fdad962e3bb13b6 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -87,11 +87,11 @@ TEST_CASE("ParticleProperties", "[Particles]") {
SECTION("Energy threshold") {
//! by default energy thresholds are set to zero
- CHECK(get_kinetic_energy_threshold(Electron::code) == 0_GeV);
+ CHECK(calculate_kinetic_energy_threshold(Electron::code) == 0_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);
+ CHECK_FALSE(calculate_kinetic_energy_threshold(Code::Electron) == 1_GeV);
+ CHECK(calculate_kinetic_energy_threshold(Code::Electron) == 10_GeV);
}
SECTION("Particle groups: electromagnetic") {
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index b293d5eb04a98c2dca8f375f094c22dbcb711ec2..ef2b6c6365056a51a2ea974e548353a2683f5baa 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -67,8 +67,8 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
// add secondaries, all with energies above the threshold
// only cut is by species
for (auto proType : particleList)
- projectile.addSecondary(std::make_tuple(
- proType, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
+ projectile.addSecondary(
+ std::make_tuple(proType, Eabove, DirectionVector(rootCS, {1, 0, 0})));
CHECK(view.getEntries() == 11);
CHECK(stack.getEntries() == 12);
@@ -94,8 +94,8 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
// add secondaries, all with energies above the threshold
// only cut is by species
for (auto proType : particleList) {
- projectile.addSecondary(std::make_tuple(
- proType, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
+ projectile.addSecondary(
+ std::make_tuple(proType, Eabove, DirectionVector(rootCS, {1, 0, 0})));
}
cut.doSecondaries(view);
@@ -118,16 +118,14 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
// add secondaries, all with energies below the threshold
// only cut is by species
for (auto proType : particleList)
- projectile.addSecondary(std::make_tuple(
- proType, Ebelow, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
+ projectile.addSecondary(
+ std::make_tuple(proType, Ebelow, DirectionVector(rootCS, {1, 0, 0})));
unsigned short A = 18;
unsigned short Z = 8;
projectile.addSecondary(std::make_tuple(get_nucleus_code(A, Z), Eabove * A,
- DirectionVector(rootCS, {1, 0, 0}), point0,
- 0_ns));
+ DirectionVector(rootCS, {1, 0, 0})));
projectile.addSecondary(std::make_tuple(get_nucleus_code(A, Z), Ebelow * A,
- DirectionVector(rootCS, {1, 0, 0}), point0,
- 0_ns));
+ DirectionVector(rootCS, {1, 0, 0})));
cut.doSecondaries(view);
@@ -148,21 +146,19 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
// only this way the secondary view is populated
auto projectile = view.getProjectile();
// add secondaries
- 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, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
- projectile.addSecondary(std::make_tuple(
- Code::PiPlus, 4_GeV, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
+ projectile.addSecondary(
+ std::make_tuple(Code::Photon, 3_MeV, DirectionVector(rootCS, {1, 0, 0})));
+ projectile.addSecondary(
+ std::make_tuple(Code::Electron, 3_MeV, DirectionVector(rootCS, {1, 0, 0})));
+ projectile.addSecondary(
+ std::make_tuple(Code::PiPlus, 4_GeV, DirectionVector(rootCS, {1, 0, 0})));
unsigned short A = 18;
unsigned short Z = 8;
projectile.addSecondary(std::make_tuple(get_nucleus_code(A, Z), 4_GeV * A,
- DirectionVector(rootCS, {1, 0, 0}), point0,
- 0_ns));
+ DirectionVector(rootCS, {1, 0, 0})));
projectile.addSecondary(std::make_tuple(get_nucleus_code(A, Z), 6_GeV * A,
- DirectionVector(rootCS, {1, 0, 0}), point0,
- 0_ns));
+ DirectionVector(rootCS, {1, 0, 0})));
cut.doSecondaries(view);
@@ -172,11 +168,11 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
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_kinetic_energy_threshold(Code::Electron) !=
- get_kinetic_energy_threshold(Code::Positron));
- CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == Electron::mass);
+ CHECK(calculate_kinetic_energy_threshold(Code::Electron) !=
+ calculate_kinetic_energy_threshold(Code::Positron));
+ CHECK_FALSE(calculate_kinetic_energy_threshold(Code::Electron) == Electron::mass);
// test default values still correct
- CHECK(get_kinetic_energy_threshold(Code::Proton) == 5_GeV);
+ CHECK(calculate_kinetic_energy_threshold(Code::Proton) == 5_GeV);
}
SECTION("cut on time") {
@@ -185,7 +181,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
// add primary particle to stack
auto particle = stack.addParticle(std::make_tuple(
- Code::Proton, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, 1_ns));
+ Code::Proton, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, too_late));
// view on secondary particles
setup::StackView view(particle);
// ref. to primary particle through the secondary view.
@@ -194,8 +190,8 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
// add secondaries, all with energies above the threshold
// only cut is by time
for (auto proType : particleList) {
- projectile.addSecondary(std::make_tuple(
- proType, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, too_late));
+ projectile.addSecondary(
+ std::make_tuple(proType, Eabove, DirectionVector(rootCS, {1, 0, 0})));
}
cut.doSecondaries(view);
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
index ec2157e3e79cff3a1680692b482bf6f5eac159f3..f76eab9a93f7177ee20cfc4bf6b531e153aba0e1 100644
--- a/tests/modules/testStackInspector.cpp
+++ b/tests/modules/testStackInspector.cpp
@@ -35,11 +35,11 @@ TEST_CASE("StackInspector", "modules") {
TestCascadeStack stack;
stack.clear();
HEPEnergyType E0 = 100_GeV;
- stack.addParticle(std::make_tuple(Code::Electron,
- MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
+ stack.addParticle(std::make_tuple(Code::Electron, 1_GeV,
+ DirectionVector(rootCS, {0, 0, -1}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
- stack.addParticle(std::make_tuple(get_nucleus_code(16, 8),
- MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
+ stack.addParticle(std::make_tuple(get_nucleus_code(16, 8), 1_GeV,
+ DirectionVector(rootCS, {0, 0, -1}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
SECTION("interface") {
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index f75d59b75b44a41d25376d013c9d9e0d813efa32..10f5954330887a57a09c7b4095311c99849feb31 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -121,8 +121,6 @@ TEST_CASE("UrQMD") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::PiPlus, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
- CHECK(stackPtr->getEntries() == 1);
- CHECK(secViewPtr->getEntries() == 0);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->getProjectile();
diff --git a/tests/stack/testVectorStack.cpp b/tests/stack/testVectorStack.cpp
index b44e3d55896e1c5d259642210f9b4dfa8de19a68..74d7d6f76c33e112fd576ac5427707300d53dff0 100644
--- a/tests/stack/testVectorStack.cpp
+++ b/tests/stack/testVectorStack.cpp
@@ -40,12 +40,6 @@ TEST_CASE("VectorStack", "stack") {
CHECK(pout.getTime() == 100_s);
CHECK(pout.getChargeNumber() == -1);
- // particle with no momentum, has no direction
- auto const p0 = s.addParticle(
- std::make_tuple(Code::Proton, MomentumVector(dummyCS, {0_GeV, 0_GeV, 0_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
- CHECK(p0.getDirection().getNorm() == 0);
-
s.clear();
CHECK(s.getEntries() == 0);
CHECK(s.getSize() == 0);
@@ -57,7 +51,7 @@ TEST_CASE("VectorStack", "stack") {
for (int i = 0; i < 99; ++i)
s.addParticle(
- std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ std::make_tuple(Code::Electron, 1_GeV, DirectionVector(dummyCS, {1, 0, 0}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
CHECK(s.getSize() == 99);
@@ -161,4 +155,36 @@ TEST_CASE("VectorStack", "stack") {
for (int i = 0; i < 99; ++i) s.last().erase();
CHECK(s.getEntries() == 0);
}
+
+ SECTION("secondary generation") {
+
+ // check inheritance of time and location
+
+ VectorStack s;
+ auto p1 = s.addParticle(std::make_tuple(
+ Code::Photon, 1.5_GeV - Electron::mass, DirectionVector(dummyCS, {1, 0, 0}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 10_s));
+ auto p2 = p1.addSecondary(
+ std::make_tuple(Code::Photon, 10_GeV, DirectionVector(dummyCS, {1, 0, 0})));
+ auto const delta_pos = Vector<length_d>(dummyCS, {1_m, 0_m, 0_m});
+ auto const delta_time = 1_s;
+ auto p3 = p1.addSecondary(std::make_tuple(Code::Photon, 10_GeV,
+ DirectionVector(dummyCS, {1, 0, 0}),
+ delta_pos, delta_time));
+
+ CHECK(p1.getTime() / 10_s == Approx(1));
+ CHECK(p1.getPosition().getX(dummyCS) / 1_m == Approx(1));
+ CHECK(p1.getPosition().getY(dummyCS) / 1_m == Approx(1));
+ CHECK(p1.getPosition().getZ(dummyCS) / 1_m == Approx(1));
+
+ CHECK(p2.getTime() / 10_s == Approx(1));
+ CHECK(p2.getPosition().getX(dummyCS) / 1_m == Approx(1));
+ CHECK(p2.getPosition().getY(dummyCS) / 1_m == Approx(1));
+ CHECK(p2.getPosition().getZ(dummyCS) / 1_m == Approx(1));
+
+ CHECK(p3.getTime() / 11_s == Approx(1));
+ CHECK(p3.getPosition().getX(dummyCS) / 1_m == Approx(2));
+ CHECK(p3.getPosition().getY(dummyCS) / 1_m == Approx(1));
+ CHECK(p3.getPosition().getZ(dummyCS) / 1_m == Approx(1));
+ }
}
\ No newline at end of file