From e177228beacd26e82d411d5e7d0724e3814ac0fa Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 1 Dec 2021 08:12:18 +0100
Subject: [PATCH] fixed recent timing bug
---
.../detail/modules/epos/InteractionModel.inl | 9 ++--
corsika/detail/modules/pythia8/Decay.inl | 10 ++---
.../detail/modules/pythia8/Interaction.inl | 6 ++-
.../modules/qgsjetII/InteractionModel.inl | 33 ++++++++-------
corsika/detail/modules/sibyll/Decay.inl | 11 ++---
.../modules/sibyll/InteractionModel.inl | 13 +++---
.../sibyll/NuclearInteractionModel.inl | 38 ++++++++++-------
corsika/detail/modules/urqmd/UrQMD.inl | 15 ++++---
corsika/detail/stack/VectorStack.inl | 41 ++----------------
corsika/stack/VectorStack.hpp | 39 +++++++----------
tests/common/SetupTestStack.hpp | 9 ++--
tests/framework/testCascade.cpp | 17 ++++----
tests/modules/testParticleCut.cpp | 42 +++++++++----------
tests/modules/testStackInspector.cpp | 8 ++--
tests/modules/testUrQMD.cpp | 2 -
tests/stack/testVectorStack.cpp | 8 +---
16 files changed, 131 insertions(+), 170 deletions(-)
diff --git a/corsika/detail/modules/epos/InteractionModel.inl b/corsika/detail/modules/epos/InteractionModel.inl
index 35ed15a88..f9e8925b7 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/pythia8/Decay.inl b/corsika/detail/modules/pythia8/Decay.inl
index 81306df45..3da0ebdbf 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 fe55d5874..2144380e7 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 8203faa49..81b62894b 100644
--- a/corsika/detail/modules/qgsjetII/InteractionModel.inl
+++ b/corsika/detail/modules/qgsjetII/InteractionModel.inl
@@ -178,11 +178,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 +203,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();
@@ -252,15 +252,19 @@ 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 = sqrt(p3output.getSquaredNorm() + mass * mass) - 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();
}
@@ -278,11 +282,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 = sqrt(p3output.getSquaredNorm() + mass * mass) - 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 ab0114182..b804feb53 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 99f078673..dd00e7606 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 4a564c7e9..06de518c4 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,13 @@ 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 1bf0f0d59..9e070f1c6 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 57bc094a4..fc62767e1 100644
--- a/corsika/detail/stack/VectorStack.inl
+++ b/corsika/detail/stack/VectorStack.inl
@@ -22,41 +22,6 @@
namespace corsika {
- template <typename StackIteratorInterface>
- inline void ParticleInterface<StackIteratorInterface>::setParticleData(
- 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(
- 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) {
@@ -70,12 +35,12 @@ namespace corsika {
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
ParticleInterface<StackIteratorInterface> const& parent,
- 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) + parent.getTime()); // parent time is added
+ this->setPosition(parent.getPosition()); // position
+ this->setTime(parent.getTime()); // parent time
}
template <typename StackIteratorInterface>
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index 22873790c..029f75490 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -33,50 +33,41 @@ 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;
-
- std::string asString() const;
-
/**
- * Set data of new particle.
- *
- * @param v tuple containing: PID, Momentum Vector, Position, Time
+ * secondary particle data information content.
*
- * MomentumVector is only used to determine the DirectionVector, the normalization
- * is lost.
+ * PID, Ekin, direction.
*/
- void setParticleData(particle_data_type const& v);
+ typedef std::tuple<Code, HEPEnergyType, DirectionVector> secondary_data_type;
+
+ std::string asString() const;
/**
* Set data of new particle.
*
- * @param parent parent 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(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);
+ void setParticleData(particle_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_data_type.
*/
void setParticleData(ParticleInterface<TStackIterator> const& parent,
- particle_data_momentum_type const& v);
+ secondary_data_type const& v);
///! Set particle corsika::Code
void setPID(Code const id) {
diff --git a/tests/common/SetupTestStack.hpp b/tests/common/SetupTestStack.hpp
index a6141e582..476c99a54 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 18732cb1f..624dc0415 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/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index b293d5eb0..fc546ebf3 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);
@@ -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 ec2157e3e..f76eab9a9 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 f75d59b75..10f595433 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 b44e3d558..2adb1d49c 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);
--
GitLab