diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl
index f91c4c3c472f564f2b6060c1a94b83fcdc1a8399..feb10651867f431b71afae8814353a11116cc4d8 100644
--- a/corsika/detail/framework/core/ParticleProperties.inl
+++ b/corsika/detail/framework/core/ParticleProperties.inl
@@ -10,35 +10,58 @@
#include <cmath>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
namespace corsika {
- inline HEPEnergyType constexpr get_kinetic_energy_threshold(Code const p) {
- return particle::detail::thresholds[static_cast<CodeIntType>(p)];
+ inline HEPEnergyType constexpr get_kinetic_energy_threshold(Code const code) {
+ if (is_nucleus(code)) return particle::detail::threshold_nuclei;
+ return particle::detail::thresholds[static_cast<CodeIntType>(code)];
}
- inline void constexpr set_kinetic_energy_threshold(Code const p,
+ inline void constexpr set_kinetic_energy_threshold(Code const code,
HEPEnergyType const val) {
- particle::detail::thresholds[static_cast<CodeIntType>(p)] = val;
+ if (is_nucleus(code))
+ particle::detail::threshold_nuclei = val;
+ else
+ particle::detail::thresholds[static_cast<CodeIntType>(code)] = val;
}
- inline HEPMassType constexpr get_mass(Code const p) {
- if (p == Code::Nucleus)
- throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
- return particle::detail::masses[static_cast<CodeIntType>(p)];
+ inline HEPMassType constexpr get_mass(Code const code) {
+ if (is_nucleus(code)) { return get_nucleus_mass(code); }
+ return particle::detail::masses[static_cast<CodeIntType>(code)];
}
- inline PDGCode constexpr get_PDG(Code const p) {
- return particle::detail::pdg_codes[static_cast<CodeIntType>(p)];
+ inline bool constexpr is_nucleus(Code const code) { return code >= Code::Nucleus; }
+
+ inline Code constexpr get_nucleus_code(unsigned int const A,
+ unsigned int const Z) { // 10LZZZAAAI
+ if (Z > A) { throw std::runtime_error("Z cannot be larger than A in nucleus."); }
+ return static_cast<Code>(static_cast<CodeIntType>(Code::Nucleus) + Z * 10000 +
+ A * 10);
+ }
+
+ inline unsigned int constexpr get_nucleus_Z(Code const code) {
+ return (static_cast<CodeIntType>(code) % static_cast<CodeIntType>(Code::Nucleus)) /
+ 10000;
+ }
+
+ inline unsigned int constexpr get_nucleus_A(Code const code) {
+ return (static_cast<CodeIntType>(code) % 10000) / 10;
}
- inline PDGCode constexpr get_PDG(unsigned int A, unsigned int Z) {
- return static_cast<PDGCode>(1000000000 + Z * 10000 + A * 10); // 10LZZZAAAI
+ inline PDGCode constexpr get_PDG(Code const code) {
+ if (code < Code::Nucleus) {
+ return particle::detail::pdg_codes[static_cast<CodeIntType>(code)];
+ }
+ unsigned int const Z = get_nucleus_Z(code);
+ unsigned int const A = get_nucleus_A(code);
+ return static_cast<PDGCode>(static_cast<CodeIntType>(Code::Nucleus) + Z * 10000 +
+ A * 10); // 10LZZZAAAI
}
inline int16_t constexpr get_charge_number(Code const code) {
- if (code == Code::Nucleus)
- throw std::runtime_error("charge of particle::Nucleus undefined");
+ if (is_nucleus(code)) return get_nucleus_Z(code);
return particle::detail::electric_charges[static_cast<CodeIntType>(code)];
}
@@ -47,6 +70,7 @@ namespace corsika {
}
inline std::string_view constexpr get_name(Code const code) {
+ if (is_nucleus(code)) { return "nucleus"; }
return particle::detail::names[static_cast<CodeIntType>(code)];
}
@@ -54,8 +78,9 @@ namespace corsika {
return particle::detail::lifetime[static_cast<CodeIntType>(p)] * second;
}
- inline bool constexpr is_hadron(Code const p) {
- return particle::detail::isHadron[static_cast<CodeIntType>(p)];
+ inline bool constexpr is_hadron(Code const code) {
+ if (is_nucleus(code)) return true;
+ return particle::detail::isHadron[static_cast<CodeIntType>(code)];
}
inline bool constexpr is_em(Code const c) {
@@ -71,24 +96,6 @@ namespace corsika {
c == Code::NuMuBar || c == Code::NuTauBar;
}
- inline int constexpr get_nucleus_A(Code const code) {
- if (code == Code::Nucleus) {
- throw std::runtime_error("get_nucleus_A(Code::Nucleus) is impossible!");
- }
- return particle::detail::nucleusA[static_cast<CodeIntType>(code)];
- }
-
- inline int constexpr get_nucleus_Z(Code const code) {
- if (code == Code::Nucleus) {
- throw std::runtime_error("get_nucleus_Z(Code::Nucleus) is impossible!");
- }
- return particle::detail::nucleusZ[static_cast<CodeIntType>(code)];
- }
-
- inline bool constexpr is_nucleus(Code const code) {
- return (code == Code::Nucleus) || (get_nucleus_A(code) != 0);
- }
-
inline std::ostream& operator<<(std::ostream& stream, corsika::Code const code) {
return stream << get_name(code);
}
@@ -97,7 +104,7 @@ namespace corsika {
static_assert(particle::detail::conversionArray.size() % 2 == 1);
// this will fail, for the strange case where the maxPDG is negative...
int constexpr maxPDG{(particle::detail::conversionArray.size() - 1) >> 1};
- auto const k = static_cast<PDGCodeType>(p);
+ auto const k = static_cast<PDGCodeIntType>(p);
if (std::abs(k) <= maxPDG) {
return particle::detail::conversionArray[k + maxPDG];
} else {
@@ -105,10 +112,18 @@ namespace corsika {
}
}
- inline HEPMassType get_nucleus_mass(unsigned int const A, unsigned int const Z) {
+ inline HEPMassType constexpr get_nucleus_mass(Code const code) {
+ unsigned int const A = get_nucleus_A(code);
+ unsigned int const Z = get_nucleus_Z(code);
return get_mass(Code::Proton) * Z + (A - Z) * get_mass(Code::Neutron);
}
+ inline std::string_view get_nucleus_name(Code const code) {
+ unsigned int const A = get_nucleus_A(code);
+ unsigned int const Z = get_nucleus_Z(code);
+ return fmt::format("Nucleus_A{}_Z{}", A, Z);
+ }
+
inline std::initializer_list<Code> constexpr get_all_particles() {
return particle::detail::all_particles;
}
diff --git a/corsika/detail/framework/process/InteractionCounter.inl b/corsika/detail/framework/process/InteractionCounter.inl
index d3572bf98ef63faf37f50d2a58cfdf4b610a4179..7a549d61da111e3a338b4b19d2b01cf7e5cc7ce2 100644
--- a/corsika/detail/framework/process/InteractionCounter.inl
+++ b/corsika/detail/framework/process/InteractionCounter.inl
@@ -25,14 +25,7 @@ namespace corsika {
.getNuclearComposition()
.getAverageMassNumber();
auto const massTarget = massNumber * constants::nucleonMass;
-
- if (auto const projectile_id = projectile.getPID(); projectile_id == Code::Nucleus) {
- auto const A = projectile.getNuclearA();
- auto const Z = projectile.getNuclearZ();
- histogram_.fill(projectile_id, projectile.getEnergy(), massTarget, A, Z);
- } else {
- histogram_.fill(projectile_id, projectile.getEnergy(), massTarget);
- }
+ histogram_.fill(projectile.getPID(), projectile.getEnergy(), massTarget);
process_.doInteraction(view);
}
diff --git a/corsika/detail/framework/process/InteractionHistogram.inl b/corsika/detail/framework/process/InteractionHistogram.inl
index efeb29c0027e98d52d49722334b9ebb51bf1460a..4f493a363bcef4a6f11a820f889c665ab0c1961b 100644
--- a/corsika/detail/framework/process/InteractionHistogram.inl
+++ b/corsika/detail/framework/process/InteractionHistogram.inl
@@ -21,32 +21,27 @@ namespace corsika {
, inthist_lab_{detail::hist_factory(num_bins_lab, lower_edge_lab, upper_edge_lab)} {
}
- inline void InteractionHistogram::fill(Code projectile_id, HEPEnergyType lab_energy,
- HEPEnergyType mass_target, int A, int Z) {
+ inline void InteractionHistogram::fill(Code const projectile_id,
+ HEPEnergyType const lab_energy,
+ HEPEnergyType const mass_target) {
auto constexpr inv_eV = 1 / 1_eV;
- if (projectile_id == Code::Nucleus) {
- auto const sqrtS = sqrt(A * A * (constants::nucleonMass * constants::nucleonMass) +
- mass_target * mass_target + 2 * lab_energy * mass_target);
-
- int32_t const pdg = 1'000'000'000l + Z * 10'000l + A * 10l;
-
- inthist_lab_(pdg, lab_energy * inv_eV);
- inthist_cms_(pdg, sqrtS * inv_eV);
- } else {
- auto const projectile_mass = get_mass(projectile_id);
- auto const sqrtS = sqrt(projectile_mass * projectile_mass +
- mass_target * mass_target + 2 * lab_energy * mass_target);
-
- inthist_cms_(static_cast<int>(get_PDG(projectile_id)), sqrtS * inv_eV);
- inthist_lab_(static_cast<int>(get_PDG(projectile_id)), lab_energy * inv_eV);
- }
+ auto const projectile_mass = get_mass(projectile_id);
+ auto const sqrtS = sqrt(projectile_mass * projectile_mass +
+ mass_target * mass_target + 2 * lab_energy * mass_target);
+
+ CORSIKA_LOG_INFO("pM={}, tM={}, pid={}, Elab={}, sqrtS={}, pdg={} a={} z={}",
+ projectile_mass / 1_GeV, mass_target / 1_GeV, projectile_id,
+ lab_energy / 1_GeV, sqrtS / 1_GeV, get_PDG(projectile_id),
+ get_nucleus_A(projectile_id), get_nucleus_Z(projectile_id));
+
+ inthist_cms_(static_cast<int>(get_PDG(projectile_id)), sqrtS * inv_eV);
+ inthist_lab_(static_cast<int>(get_PDG(projectile_id)), lab_energy * inv_eV);
}
inline InteractionHistogram& InteractionHistogram::operator+=(
InteractionHistogram const& other) {
inthist_lab_ += other.inthist_lab_;
inthist_cms_ += other.inthist_cms_;
-
return *this;
}
@@ -54,7 +49,6 @@ namespace corsika {
InteractionHistogram other) const {
other.inthist_lab_ += inthist_lab_;
other.inthist_cms_ += inthist_cms_;
-
return other;
}
} // namespace corsika
diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl
index af19f42f09e043c3d5a248c4d2c67f6cd70e262c..0e6163af59d4f5b16a282695d4c27bbe0f460a61 100644
--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -25,7 +25,7 @@ namespace corsika {
, em_count_(0)
, inv_count_(0)
, energy_count_() {
- for (auto p : get_all_particles())
+ for (auto p : get_all_particles()) {
if (is_hadron(p)) // nuclei are also hadrons
set_kinetic_energy_threshold(p, eHadCut);
else if (is_muon(p))
@@ -34,6 +34,8 @@ namespace corsika {
set_kinetic_energy_threshold(p, eEleCut);
else if (p == Code::Photon)
set_kinetic_energy_threshold(p, ePhoCut);
+ }
+ set_kinetic_energy_threshold(Code::Nucleus, eHadCut);
CORSIKA_LOG_DEBUG(
"setting kinetic energy thresholds: electrons = {} GeV, photons = {} GeV, "
"hadrons = {} GeV, "
@@ -53,11 +55,13 @@ namespace corsika {
, inv_count_(0)
, energy_count_(0) {
- for (auto p : get_all_particles())
+ for (auto p : get_all_particles()) {
if (is_hadron(p))
set_kinetic_energy_threshold(p, eHadCut);
else if (is_muon(p))
set_kinetic_energy_threshold(p, eMuCut);
+ }
+ set_kinetic_energy_threshold(Code::Nucleus, eHadCut);
CORSIKA_LOG_DEBUG(
"setting thresholds: hadrons = {} GeV, "
"muons = {} GeV",
@@ -75,6 +79,7 @@ namespace corsika {
, inv_count_(0)
, energy_count_(0) {
for (auto p : get_all_particles()) set_kinetic_energy_threshold(p, eCut);
+ set_kinetic_energy_threshold(Code::Nucleus, eCut);
CORSIKA_LOG_DEBUG("setting kinetic energy threshold for all particles to {} GeV",
eCut / 1_GeV);
}
@@ -100,9 +105,9 @@ namespace corsika {
auto const energyLab = vP.getKineticEnergy();
auto const pid = vP.getPID();
// nuclei
- if (pid == Code::Nucleus) {
+ if (is_nucleus(pid)) {
// calculate energy per nucleon
- auto const ElabNuc = energyLab / vP.getNuclearA();
+ auto const ElabNuc = energyLab / get_nucleus_A(pid);
return (ElabNuc < get_kinetic_energy_threshold(pid));
} else {
return (energyLab < get_kinetic_energy_threshold(pid));
diff --git a/corsika/detail/modules/StackInspector.inl b/corsika/detail/modules/StackInspector.inl
index 5ab9fb6a8a258e9f1747fdb02abb8ffe6d27869c..24a356fcb749f90796db5920bef351dee2d10ac9 100644
--- a/corsika/detail/modules/StackInspector.inl
+++ b/corsika/detail/modules/StackInspector.inl
@@ -51,9 +51,6 @@ namespace corsika {
" pos= {}"
" node = {}",
(i++), iterP.getPID(), (E / 1_GeV), pos, fmt::ptr(iterP.getNode()));
-
- if (iterP.getPID() == Code::Nucleus)
- CORSIKA_LOG_INFO("nuc_ref= {}", iterP.getNucleusRef());
}
}
diff --git a/corsika/detail/modules/epos/Interaction.inl b/corsika/detail/modules/epos/Interaction.inl
index 058dd2c7f21c3961c1de0030d3a0c86942123a28..dde502eb0b3284a83b1fd5af01708800d6e49573 100644
--- a/corsika/detail/modules/epos/Interaction.inl
+++ b/corsika/detail/modules/epos/Interaction.inl
@@ -63,21 +63,7 @@ namespace corsika::epos {
}
inline bool Interaction::isValidTarget(Code const TargetId) const {
- if (is_nucleus(TargetId)) {
- if (TargetId == Code::Nucleus) {
- // nuclearExtension for projectiles only
- CORSIKA_LOGGER_WARN(logger_,
- "Invalid target!"
- " Code::Nucleus only allowed for "
- "projectiles! "
- "This should not happen!");
- return false;
- } else {
- return (get_nucleus_Z(TargetId) < maxTargetMassNumber_ ? true : false);
- }
- } else {
- return false;
- }
+ return is_nucleus(TargetId) && (get_nucleus_A(TargetId) < maxTargetMassNumber_);
}
inline void Interaction::initialize() const {
@@ -418,8 +404,8 @@ namespace corsika::epos {
int beamA = 1;
int beamZ = 1;
if (is_nucleus(corsikaBeamId)) {
- beamA = projectile.getNuclearA();
- beamZ = projectile.getNuclearZ();
+ beamA = get_nucleus_A(corsikaBeamId);
+ beamZ = get_nucleus_Z(corsikaBeamId);
}
// get target from environment
@@ -485,8 +471,8 @@ namespace corsika::epos {
int beamA = 1;
int beamZ = 1;
if (is_nucleus(corsikaBeamId)) {
- beamA = projectile.getNuclearA();
- beamZ = projectile.getNuclearZ();
+ beamA = get_nucleus_A(corsikaBeamId);
+ beamZ = get_nucleus_Z(corsikaBeamId);
CORSIKA_LOGGER_DEBUG(logger_, "A={}, Z={} ", beamA, beamZ);
}
@@ -579,13 +565,11 @@ namespace corsika::epos {
A = 4;
Z = 2;
} else {
- // 100ZZZAAA0 -> std. pdg code
- EposCodeIntType const eposPdg = static_cast<EposCodeIntType>(eposId);
- Z = int(abs(eposPdg) / 10000) % 1000;
- A = int(abs(eposPdg) / 10) % 1000;
+ Z = get_nucleus_Z(eposId);
+ A = get_nucleus_Z(eposId);
}
auto pnew = view.addSecondary(
- std::make_tuple(Code::Nucleus, momentum, pOrig, tOrig, A, Z));
+ std::make_tuple(get_nucleus_code(A, Z), momentum, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
diff --git a/corsika/detail/modules/epos/ParticleConversion.inl b/corsika/detail/modules/epos/ParticleConversion.inl
index 38c52adbe9e0e89b19dd574dcc8f6c490394153a..b3f34dbfef4654e164666740cc9458c219f53526 100644
--- a/corsika/detail/modules/epos/ParticleConversion.inl
+++ b/corsika/detail/modules/epos/ParticleConversion.inl
@@ -14,8 +14,7 @@
namespace corsika::epos {
inline HEPMassType getEposMass(Code const pCode) {
- if (pCode == Code::Nucleus)
- throw std::runtime_error("Cannot getMass() of particle::Nucleus -> unspecified");
+ if (is_nucleus(pCode)) throw std::runtime_error("Not suited for Nuclei.");
auto sCode = convertToEposRaw(pCode);
if (sCode == 0)
throw std::runtime_error("getEposMass: unknown particle!");
diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl
index 1ee0f978d772159780968bab6625729606acfec0..c3147dbb83b00aff24e6f909220ff88586360ed4 100644
--- a/corsika/detail/modules/qgsjetII/Interaction.inl
+++ b/corsika/detail/modules/qgsjetII/Interaction.inl
@@ -56,7 +56,7 @@ namespace corsika::qgsjetII {
int iTarget = 1;
if (is_nucleus(targetId)) {
iTarget = targetA;
- if (iTarget > maxMassNumber_ || iTarget <= 0) {
+ if (iTarget > int(maxMassNumber_) || iTarget <= 0) {
std::ostringstream txt;
txt << "QgsjetII target outside range. Atarget=" << iTarget;
throw std::runtime_error(txt.str().c_str());
@@ -65,7 +65,7 @@ namespace corsika::qgsjetII {
int iProjectile = 1;
if (is_nucleus(beamId)) {
iProjectile = Abeam;
- if (iProjectile > maxMassNumber_ || iProjectile <= 0) {
+ if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) {
std::ostringstream txt;
txt << "QgsjetII projectile outside range. Aprojectile=" << iProjectile;
throw std::runtime_error(txt.str().c_str());
@@ -84,12 +84,12 @@ namespace corsika::qgsjetII {
}
template <typename TParticle>
- inline GrammageType Interaction::getInteractionLength(const TParticle& vP) const {
+ inline GrammageType Interaction::getInteractionLength(const TParticle& particle) const {
// coordinate system, get global frame of reference
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
- const Code corsikaBeamId = vP.getPID();
+ const Code corsikaBeamId = particle.getPID();
// beam particles for qgsjetII : 1, 2, 3 for p, pi, k
// read from cross section code table
@@ -99,19 +99,19 @@ namespace corsika::qgsjetII {
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
- HEPEnergyType Elab = vP.getEnergy();
+ HEPEnergyType const Elab = particle.getEnergy();
CORSIKA_LOG_DEBUG(
"Interaction: LambdaInt: \n"
" input energy: {} GeV"
" beam can interact: {}"
" beam pid: {}",
- vP.getEnergy() / 1_GeV, kInteraction, vP.getPID());
+ particle.getEnergy() / 1_GeV, kInteraction, corsikaBeamId);
if (kInteraction) {
int Abeam = 0;
- if (is_nucleus(vP.getPID())) Abeam = vP.getNuclearA();
+ if (is_nucleus(corsikaBeamId)) Abeam = get_nucleus_A(corsikaBeamId);
// get target from environment
/*
@@ -120,7 +120,7 @@ namespace corsika::qgsjetII {
and the boosts can be defined..
*/
- auto const* currentNode = vP.getNode();
+ auto const* currentNode = particle.getNode();
const auto& mediumComposition =
currentNode->getModelProperties().getNuclearComposition();
@@ -159,184 +159,177 @@ namespace corsika::qgsjetII {
inline void Interaction::doInteraction(TView& view) {
auto const projectile = view.getProjectile();
-
- const auto corsikaBeamId = projectile.getPID();
+ auto const corsikaBeamId = projectile.getPID();
CORSIKA_LOG_DEBUG(
"ProcessQgsjetII: "
- "DoInteraction: {} interaction possible? {}",
+ "doInteraction: {} interaction possible? {}",
corsikaBeamId, corsika::qgsjetII::canInteract(corsikaBeamId));
- if (corsika::qgsjetII::canInteract(corsikaBeamId)) {
+ if (!corsika::qgsjetII::canInteract(corsikaBeamId)) return;
- CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
+ CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
- // position and time of interaction, not used in QgsjetII
- Point pOrig = projectile.getPosition();
- TimeType tOrig = projectile.getTime();
+ // position and time of interaction, not used in QgsjetII
+ Point const pOrig = projectile.getPosition();
+ TimeType const tOrig = projectile.getTime();
- // define target
- // for QgsjetII is always a single nucleon
- // FOR NOW: target is always at rest
- auto const targetEnergyLab = 0_GeV + constants::nucleonMass;
- auto const targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
- FourVector const PtargLab(targetEnergyLab, targetMomentumLab);
+ // define target
+ // for QgsjetII is always a single nucleon
+ // FOR NOW: target is always at rest
+ auto const targetEnergyLab = 0_GeV + constants::nucleonMass;
+ auto const targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
+ FourVector const PtargLab(targetEnergyLab, targetMomentumLab);
- // define projectile
- HEPEnergyType const projectileEnergyLab = projectile.getEnergy();
- auto const projectileMomentumLab = projectile.getMomentum();
+ // define projectile
+ HEPEnergyType const projectileEnergyLab = projectile.getEnergy();
+ auto const projectileMomentumLab = projectile.getMomentum();
- int beamA = 0;
- if (is_nucleus(corsikaBeamId)) beamA = projectile.getNuclearA();
+ int beamA = 0;
+ if (is_nucleus(corsikaBeamId)) beamA = get_nucleus_A(corsikaBeamId);
- HEPEnergyType const projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
+ HEPEnergyType const projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
- CORSIKA_LOG_DEBUG(
- "Interaction: ebeam lab: {} GeV"
- "Interaction: pbeam lab: {} GeV",
- projectileEnergyLab / 1_GeV, projectileMomentumLab.getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG(
- "Interaction: etarget lab: {} GeV"
- "Interaction: ptarget lab: {} GeV",
- targetEnergyLab / 1_GeV, targetMomentumLab.getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG("Interaction: position of interaction: {}",
- pOrig.getCoordinates());
- CORSIKA_LOG_DEBUG("Interaction: time: {} ", tOrig);
-
- // sample target mass number
- auto const* currentNode = projectile.getNode();
- auto const& mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
- // get cross sections for target materials
- /*
- Here we read the cross section from the interaction model again,
- should be passed from getInteractionLength if possible
- */
- auto const& compVec = mediumComposition.getComponents();
- std::vector<CrossSectionType> cross_section_of_components(compVec.size());
-
- for (size_t i = 0; i < compVec.size(); ++i) {
- auto const targetId = compVec[i];
- int targetA = 0;
- if (is_nucleus(targetId)) targetA = get_nucleus_A(targetId);
- const auto sigProd =
- getCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA);
- cross_section_of_components[i] = sigProd;
- }
+ CORSIKA_LOG_DEBUG(
+ "ebeam lab: {} GeV "
+ "pbeam lab: {} GeV ",
+ projectileEnergyLab / 1_GeV, projectileMomentumLab.getComponents() / 1_GeV);
+ CORSIKA_LOG_DEBUG(
+ "etarget lab: {} GeV "
+ "ptarget lab: {} GeV ",
+ targetEnergyLab / 1_GeV, targetMomentumLab.getComponents() / 1_GeV);
+ CORSIKA_LOG_DEBUG("position of interaction: {}", pOrig.getCoordinates());
+ CORSIKA_LOG_DEBUG("time: {} ", tOrig);
+
+ // sample target mass number
+ auto const* currentNode = projectile.getNode();
+ auto const& mediumComposition =
+ currentNode->getModelProperties().getNuclearComposition();
+ // get cross sections for target materials
+ /*
+ Here we read the cross section from the interaction model again,
+ should be passed from getInteractionLength if possible
+ */
+ auto const& compVec = mediumComposition.getComponents();
+ std::vector<CrossSectionType> cross_section_of_components(compVec.size());
+
+ for (size_t i = 0; i < compVec.size(); ++i) {
+ auto const targetId = compVec[i];
+ int targetA = 0;
+ if (is_nucleus(targetId)) targetA = get_nucleus_A(targetId);
+ const auto sigProd =
+ getCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA);
+ cross_section_of_components[i] = sigProd;
+ }
- const auto targetCode =
- mediumComposition.sampleTarget(cross_section_of_components, rng_);
- CORSIKA_LOG_DEBUG("Interaction: target selected: {}", targetCode);
-
- int targetMassNumber = 1; // proton
- if (is_nucleus(targetCode)) { // nucleus
- targetMassNumber = get_nucleus_A(targetCode);
- if (targetMassNumber > maxMassNumber_)
- throw std::runtime_error(
- "QgsjetII target mass outside range."); // LCOV_EXCL_LINE there is no
- // allowed path here
- } else {
- if (targetCode != Proton::code) // LCOV_EXCL_LINE there is no allowed path here
- throw std::runtime_error(
- "QgsjetII Taget not possible."); // LCOV_EXCL_LINE there is no allowed path
- // here
- }
- CORSIKA_LOG_DEBUG("Interaction: target qgsjetII code/A: {}", targetMassNumber);
-
- int projectileMassNumber = 1; // "1" means "hadron"
- QgsjetIIHadronType qgsjet_hadron_type =
- qgsjetII::getQgsjetIIHadronType(corsikaBeamId);
- if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
- projectileMassNumber = projectile.getNuclearA();
- if (projectileMassNumber > maxMassNumber_)
- throw std::runtime_error(
- "QgsjetII projectile mass outside range."); // LCOV_EXCL_LINE there is no
- // allowed path here
- std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
- QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
- std::uniform_int_distribution select(0, 1);
- qgsjet_hadron_type = nucleons[select(rng_)];
- } else {
- // from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
- if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
- qgsjet_hadron_type = alternate_;
- alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
- ? QgsjetIIHadronType::PiMinusType
- : QgsjetIIHadronType::PiPlusType);
- }
+ const auto targetCode =
+ mediumComposition.sampleTarget(cross_section_of_components, rng_);
+
+ int targetMassNumber = 1; // proton
+ if (is_nucleus(targetCode)) { // nucleus
+ targetMassNumber = get_nucleus_A(targetCode);
+ if (targetMassNumber > int(maxMassNumber_))
+ throw std::runtime_error(
+ "QgsjetII target mass outside range."); // LCOV_EXCL_LINE there is no
+ // allowed path here
+ } else {
+ if (targetCode != Proton::code) // LCOV_EXCL_LINE there is no allowed path here
+ throw std::runtime_error(
+ "QgsjetII Taget not possible."); // LCOV_EXCL_LINE there is no allowed path
+ // here
+ }
+ CORSIKA_LOG_DEBUG("target: {}, qgsjetII code/A: {}", targetCode, targetMassNumber);
+
+ int projectileMassNumber = 1; // "1" means "hadron"
+ QgsjetIIHadronType qgsjet_hadron_type =
+ qgsjetII::getQgsjetIIHadronType(corsikaBeamId);
+ if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
+ projectileMassNumber = get_nucleus_A(corsikaBeamId);
+ if (projectileMassNumber > int(maxMassNumber_))
+ throw std::runtime_error(
+ "QgsjetII projectile mass outside range."); // LCOV_EXCL_LINE there is no
+ // allowed path here
+ std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
+ QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
+ std::uniform_int_distribution select(0, 1);
+ qgsjet_hadron_type = nucleons[select(rng_)];
+ } else {
+ // from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
+ if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
+ qgsjet_hadron_type = alternate_;
+ alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
+ ? QgsjetIIHadronType::PiMinusType
+ : QgsjetIIHadronType::PiPlusType);
}
+ }
- // beam id for qgsjetII
- int kBeam = 2; // default: proton Shouldn't we randomize neutron/proton for nuclei?
- if (corsikaBeamId != Code::Nucleus) {
- kBeam = corsika::qgsjetII::convertToQgsjetIIRaw(corsikaBeamId);
- // from conex
- if (kBeam == 0) { // replace pi0 or rho0 with pi+/pi-
- static int select = 1;
- kBeam = select;
- select *= -1;
- }
- // replace lambda by neutron
- if (kBeam == 6)
- kBeam = 3;
- else if (kBeam == -6)
- kBeam = -3;
- // else if (abs(kBeam)>6) -> throw
+ // beam id for qgsjetII
+ int kBeam = 2; // default: proton Shouldn't we randomize neutron/proton for nuclei?
+ if (!is_nucleus(corsikaBeamId)) {
+ kBeam = corsika::qgsjetII::convertToQgsjetIIRaw(corsikaBeamId);
+ // from conex
+ if (kBeam == 0) { // replace pi0 or rho0 with pi+/pi-
+ static int select = 1;
+ kBeam = select;
+ select *= -1;
}
+ // replace lambda by neutron
+ if (kBeam == 6)
+ kBeam = 3;
+ else if (kBeam == -6)
+ kBeam = -3;
+ // else if (abs(kBeam)>6) -> throw
+ }
+
+ count_++;
+ int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
+ CORSIKA_LOG_DEBUG(
+ "qgsjet_hadron_type_int={} projectileMassNumber={} targetMassNumber={}",
+ qgsjet_hadron_type_int, projectileMassNumber, targetMassNumber);
+ qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
+ targetMassNumber);
+ qgconf_();
+
+ // bookkeeping
+ MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV;
+
+ // to read the secondaries
+ // define rotation to and from CoM frame
+ // CoM frame definition in QgsjetII projectile: +z
+ auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
+ CoordinateSystemPtr const zAxisFrame =
+ make_rotationToZ(originalCS, projectileMomentumLab);
+
+ // fragments
+ QGSJetIIFragmentsStack qfs;
+ for (auto& fragm : qfs) {
+ int const A = fragm.getFragmentSize();
+ if (A == 1) { // nucleon
+ std::uniform_real_distribution<double> select;
+ Code idFragm = Code::Proton;
+ if (select(rng_) > 0.5) { idFragm = Code::Neutron; }
+
+ const HEPMassType nucleonMass = get_mass(idFragm);
+ // no pT, frgments just go forward
+ auto momentum =
+ Vector(zAxisFrame, corsika::QuantityVector<hepmomentum_d>{
+ 0.0_GeV, 0.0_GeV,
+ sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
+ (projectileEnergyLabPerNucleon - 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, momentum, pOrig, tOrig));
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
+
+ } else { // nucleus, A>1
- CORSIKA_LOG_DEBUG("Interaction: DoInteraction: E(GeV): {} GeV",
- projectileEnergyLab / 1_GeV);
- count_++;
- int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
- qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
- targetMassNumber);
- qgconf_();
-
- // bookkeeping
- MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV;
-
- // to read the secondaries
- // define rotation to and from CoM frame
- // CoM frame definition in QgsjetII projectile: +z
- auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
- CoordinateSystemPtr const zAxisFrame =
- make_rotationToZ(originalCS, projectileMomentumLab);
-
- // fragments
- QGSJetIIFragmentsStack qfs;
- for (auto& fragm : qfs) {
- Code idFragm = Code::Nucleus;
- int A = fragm.getFragmentSize();
int Z = 0;
switch (A) {
- case 1: { // proton/neutron
- std::uniform_real_distribution<double> select;
- if (select(rng_) > 0.5) {
- idFragm = Code::Proton;
- Z = 1;
- } else {
- idFragm = Code::Neutron;
- Z = 0;
- }
-
- const HEPMassType nucleonMass = get_mass(idFragm);
-
- auto momentum = Vector(
- zAxisFrame, corsika::QuantityVector<hepmomentum_d>{
- 0.0_GeV, 0.0_GeV,
- sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
- (projectileEnergyLabPerNucleon - 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, momentum, pOrig, tOrig));
- Plab_final += pnew.getMomentum();
- Elab_final += pnew.getEnergy();
- } break;
case 2: // deuterium
Z = 1;
break;
@@ -352,58 +345,53 @@ namespace corsika::qgsjetII {
}
}
- if (idFragm == Code::Nucleus) { // thus: not p or n
- const HEPMassType nucleusMass = Proton::mass * Z + Neutron::mass * (A - Z);
- auto momentum = Vector(
- zAxisFrame, QuantityVector<hepmomentum_d>{
- 0.0_GeV, 0.0_GeV,
- sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
- (projectileEnergyLabPerNucleon * A - nucleusMass))});
-
- momentum.rebase(originalCS); // transform back into standard lab frame
- CORSIKA_LOG_DEBUG(
- "secondary fragment> id={}"
- " p={}"
- " A={}"
- " Z= {}",
- idFragm, momentum.getComponents(), A, Z);
-
- auto pnew =
- view.addSecondary(std::make_tuple(idFragm, momentum, pOrig, tOrig, A, Z));
- Plab_final += pnew.getMomentum();
- Elab_final += pnew.getEnergy();
- }
- }
-
- // secondaries
- QGSJetIIStack qs;
- for (auto& psec : qs) {
-
- auto momentum = psec.getMomentum(zAxisFrame);
+ HEPMassType const nucleusMass = Proton::mass * Z + Neutron::mass * (A - Z);
+ // no pT, frgments just go forward
+ auto momentum = Vector(
+ zAxisFrame, QuantityVector<hepmomentum_d>{
+ 0.0_GeV, 0.0_GeV,
+ sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
+ (projectileEnergyLabPerNucleon * A - nucleusMass))});
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(
- "secondary fragment> id= {}"
- " p= {}",
- corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
- momentum.getComponents());
+ "secondary fragment> id={}"
+ " p={}"
+ " A={}"
+ " Z={}",
+ get_nucleus_code(A, Z), momentum.getComponents(), A, Z);
+
auto pnew = view.addSecondary(
- std::make_tuple(corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
- momentum, pOrig, tOrig));
+ std::make_tuple(get_nucleus_code(A, Z), momentum, pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
- CORSIKA_LOG_DEBUG(
- "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/
- "Elab_final="
- ", Plab_final={}"
- ", N_wounded,targ={}"
- ", N_wounded,proj={}"
- ", N_fragm,proj={}",
- Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents(),
- QGSJetIIFragmentsStackData::getWoundedNucleonsTarget(),
- QGSJetIIFragmentsStackData::getWoundedNucleonsProjectile(), qfs.getSize());
}
- }
+ // secondaries
+ QGSJetIIStack qs;
+ for (auto& psec : qs) {
+
+ auto momentum = psec.getMomentum(zAxisFrame);
+
+ momentum.rebase(originalCS); // transform back into standard lab frame
+ CORSIKA_LOG_DEBUG("secondary> id= {}, p= {}",
+ corsika::qgsjetII::convertFromQgsjetII(psec.getPID()),
+ momentum.getComponents());
+ auto pnew = view.addSecondary(std::make_tuple(
+ corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), momentum, pOrig, tOrig));
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
+ }
+ CORSIKA_LOG_DEBUG(
+ "conservation (all GeV): Ecm_final= n/a " /* << Ecm_final / 1_GeV*/
+ ", Elab_final={} "
+ ", Plab_final={}"
+ ", N_wounded,targ={}"
+ ", N_wounded,proj={}"
+ ", N_fragm,proj={}",
+ Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents(),
+ QGSJetIIFragmentsStackData::getWoundedNucleonsTarget(),
+ QGSJetIIFragmentsStackData::getWoundedNucleonsProjectile(), qfs.getSize());
+ }
} // namespace corsika::qgsjetII
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index c6adf5e105c3ca496658d6ea0c71a43a7a18e481..e4fdeec751684a0f143fbcc4420a51132bea5cf4 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -262,7 +262,7 @@ namespace corsika::sibyll {
if (is_nucleus(targetCode)) targetSibCode = get_nucleus_A(targetCode);
if (targetCode == Proton::code) targetSibCode = 1;
CORSIKA_LOG_DEBUG("Interaction: sibyll code: {}", targetSibCode);
- if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
+ if (targetSibCode > int(maxTargetMassNumber_) || targetSibCode < 1)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed.");
diff --git a/corsika/detail/modules/sibyll/NuclearInteraction.inl b/corsika/detail/modules/sibyll/NuclearInteraction.inl
index a75207f1c1372abfce16ff8f3a65f6e04140ccf0..db7788a117dd5a6cb5cb75422a10e37ebee0cffb 100644
--- a/corsika/detail/modules/sibyll/NuclearInteraction.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteraction.inl
@@ -160,11 +160,12 @@ namespace corsika::sibyll {
std::tuple<CrossSectionType, CrossSectionType> inline NuclearInteraction<
TEnvironment>::getCrossSection(TParticle const& projectile, Code const TargetId) {
- if (projectile.getPID() != Code::Nucleus)
+ if (!is_nucleus(projectile.getPID())) {
throw std::runtime_error(
"NuclearInteraction: getCrossSection: particle not a nucleus!");
+ }
- unsigned int const iBeamA = projectile.getNuclearA();
+ unsigned int const iBeamA = get_nucleus_A(projectile.getPID());
HEPEnergyType LabEnergyPerNuc = projectile.getEnergy() / iBeamA;
CORSIKA_LOG_DEBUG(
"NuclearInteraction: getCrossSection: called with: beamNuclA={} "
@@ -205,16 +206,9 @@ namespace corsika::sibyll {
const Code corsikaBeamId = projectile.getPID();
- if (corsikaBeamId != Code::Nucleus) {
- // check if target-style nucleus (enum), these are not allowed as projectile
- if (is_nucleus(corsikaBeamId)) {
- throw std::runtime_error(
- "NuclearInteraction: getInteractionLength: Wrong nucleus type. Nuclear "
- "projectiles should use NuclearStackExtension!");
- } else {
- // no nuclear interaction
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
+ if (!is_nucleus(corsikaBeamId)) {
+ // no nuclear interaction
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
// read from cross section code table
@@ -227,7 +221,7 @@ namespace corsika::sibyll {
// total momentum and energy
HEPEnergyType Elab = projectile.getEnergy() + constants::nucleonMass;
- int const nuclA = projectile.getNuclearA();
+ int const nuclA = get_nucleus_A(corsikaBeamId);
auto const ElabNuc = projectile.getEnergy() / nuclA;
MomentumVector pTotLab(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
@@ -320,14 +314,13 @@ namespace corsika::sibyll {
get_name(ProjId));
// check if target-style nucleus (enum)
- if (ProjId != Code::Nucleus)
+ if (!is_nucleus(ProjId)) {
throw std::runtime_error(
"NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
"should use NuclearStackExtension!");
+ }
- auto const ProjMass =
- projectile.getNuclearZ() * Proton::mass +
- (projectile.getNuclearA() - projectile.getNuclearZ()) * Neutron::mass;
+ auto const ProjMass = get_mass(ProjId);
CORSIKA_LOG_DEBUG("NuclearInteraction: projectile mass: {} ", ProjMass / 1_GeV);
count_++;
@@ -340,7 +333,7 @@ namespace corsika::sibyll {
CORSIKA_LOG_DEBUG("Interaction: time: {} ", tOrig / 1_s);
// projectile nucleon number
- const unsigned int kAProj = projectile.getNuclearA();
+ const unsigned int kAProj = get_nucleus_A(ProjId);
if (kAProj > getMaxNucleusAProjectile())
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
@@ -510,19 +503,14 @@ namespace corsika::sibyll {
for (int j = 0; j < nFragments; ++j) {
CORSIKA_LOG_DEBUG("fragment {}: A={} px={} py={} pz={}", j, AFragments[j],
fragments_.ppp[j][0], fragments_.ppp[j][1], fragments_.ppp[j][2]);
- Code specCode;
const auto nuclA = AFragments[j];
// get Z from stability line
const auto nuclZ = int(nuclA / 2.15 + 0.7);
// TODO: do we need to catch single nucleons??
- if (nuclA == 1)
- // TODO: sample neutron or proton
- specCode = Code::Proton;
- else
- specCode = Code::Nucleus;
-
- const HEPMassType mass = get_nucleus_mass(nuclA, nuclZ);
+ Code specCode = Code::Neutron; // sample neutron or proton ?
+ if (nuclA > 1) specCode = get_nucleus_code(nuclA, nuclZ);
+ HEPMassType const mass = get_mass(specCode);
CORSIKA_LOG_DEBUG("NuclearInteraction: adding fragment: {}", get_name(specCode));
CORSIKA_LOG_DEBUG("NuclearInteraction: A,Z: {}, {}", nuclA, nuclZ);
@@ -539,14 +527,8 @@ namespace corsika::sibyll {
CORSIKA_LOG_DEBUG("NuclearInteraction: fragment momentum: {}",
Plab.getSpaceLikeComponents().getComponents() / 1_GeV);
- if (nuclA == 1)
- // add nucleon
- projectile.addSecondary(
- std::make_tuple(specCode, Plab.getSpaceLikeComponents(), pOrig, tOrig));
- else
- // add nucleus
- projectile.addSecondary(std::make_tuple(specCode, Plab.getSpaceLikeComponents(),
- pOrig, tOrig, nuclA, nuclZ));
+ projectile.addSecondary(
+ std::make_tuple(specCode, Plab.getSpaceLikeComponents(), pOrig, tOrig));
}
// add elastic nucleons to corsika stack
diff --git a/corsika/detail/modules/sibyll/ParticleConversion.inl b/corsika/detail/modules/sibyll/ParticleConversion.inl
index d4459a7ba026d804ee3eb6aecc90cdeb81c2feb0..3bb07d5d711be5a92e376c1ac0553a535cd1630d 100644
--- a/corsika/detail/modules/sibyll/ParticleConversion.inl
+++ b/corsika/detail/modules/sibyll/ParticleConversion.inl
@@ -15,8 +15,7 @@
namespace corsika::sibyll {
inline HEPMassType getSibyllMass(Code const pCode) {
- if (pCode == corsika::Code::Nucleus)
- throw std::runtime_error("Cannot getMass() of particle::Nucleus -> unspecified");
+ if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
auto sCode = convertToSibyllRaw(pCode);
if (sCode == 0)
throw std::runtime_error("getSibyllMass: unknown particle!");
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index 12b9edbc372705ce20151b77fe855f8dc1e65024..80da719edae6312513999344fae7b819999980b9 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -129,7 +129,7 @@ namespace corsika::urqmd {
int vAProjectile = 1) {
// the following is a translation of ptsigtot() into C++
- if (vProjectileCode != Code::Nucleus &&
+ if (!is_nucleus(vProjectileCode) &&
!is_nucleus(vTargetCode)) { // both particles are "special"
auto const mProj = get_mass(vProjectileCode);
auto const mTar = get_mass(vTargetCode);
@@ -165,7 +165,7 @@ namespace corsika::urqmd {
Code targetCode) const {
auto const projectileCode = projectile.getPID();
- if (projectileCode == Code::Nucleus) {
+ if (is_nucleus(projectileCode)) {
/*
* unfortunately unavoidable at the moment until we have tools to get the actual
* inealstic cross-section from UrQMD
@@ -215,7 +215,7 @@ namespace corsika::urqmd {
auto projectile = view.getProjectile();
- auto projectileCode = projectile.getPID();
+ Code projectileCode = projectile.getPID();
auto const projectileEnergyLab = projectile.getEnergy();
auto const& projectileMomentumLab = projectile.getMomentum();
auto const& projectilePosition = projectile.getPosition();
@@ -246,14 +246,14 @@ namespace corsika::urqmd {
::urqmd::sys_.nsteps = 1;
// initialization regarding projectile
- if (Code::Nucleus == projectileCode) {
+ if (is_nucleus(projectileCode)) {
// is this everything?
::urqmd::inputs_.prspflg = 0;
- ::urqmd::sys_.Ap = projectile.getNuclearA();
- ::urqmd::sys_.Zp = projectile.getNuclearZ();
- ::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV) /
- projectile.getNuclearA();
+ ::urqmd::sys_.Ap = get_nucleus_A(projectileCode);
+ ::urqmd::sys_.Zp = get_nucleus_Z(projectileCode);
+ ::urqmd::rsys_.ebeam =
+ (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV) / ::urqmd::sys_.Ap;
::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) +
::urqmd::nucrad_(::urqmd::sys_.Ap) +
diff --git a/corsika/detail/setup/SetupStack.hpp b/corsika/detail/setup/SetupStack.hpp
index 1f62ab4baaa245cb47420aee20fa53dc3617ae5a..7a1327b2e54f6dfc81824e001d2a290a4f612e8b 100644
--- a/corsika/detail/setup/SetupStack.hpp
+++ b/corsika/detail/setup/SetupStack.hpp
@@ -10,7 +10,7 @@
#include <corsika/framework/stack/CombinedStack.hpp>
#include <corsika/stack/GeometryNodeStackExtension.hpp>
-#include <corsika/stack/NuclearStackExtension.hpp>
+#include <corsika/stack/VectorStack.hpp>
#include <corsika/stack/WeightStackExtension.hpp>
#include <corsika/stack/history/HistorySecondaryProducer.hpp>
#include <corsika/stack/history/HistoryStackExtension.hpp>
@@ -33,11 +33,11 @@ namespace corsika {
// combine particle data stack with geometry information for tracking
template <typename TStackIter>
using StackWithGeometryInterface =
- CombinedParticleInterface<nuclear_stack::ParticleDataStack::pi_type,
- SetupGeometryDataInterface, TStackIter>;
+ CombinedParticleInterface<VectorStack::pi_type, SetupGeometryDataInterface,
+ TStackIter>;
using StackWithGeometry =
- CombinedStack<typename nuclear_stack::ParticleDataStack::stack_data_type,
+ CombinedStack<typename VectorStack::stack_data_type,
node::GeometryData<setup::Environment>, StackWithGeometryInterface,
DefaultSecondaryProducer>;
diff --git a/corsika/detail/stack/NuclearStackExtension.inl b/corsika/detail/stack/NuclearStackExtension.inl
index fa0d5941fb0d7a8fc38a6371bb5475a3015e3689..cdd377d12ca9452e601847f9ad151e43c05fb4ad 100644
--- a/corsika/detail/stack/NuclearStackExtension.inl
+++ b/corsika/detail/stack/NuclearStackExtension.inl
@@ -126,7 +126,7 @@ namespace corsika::nuclear_stack {
setNuclearZ(Z);
HEPMassType m = 0_GeV;
if (PID == Code::Nucleus) {
- m = get_nucleus_mass(A, Z);
+ m = get_nucleus_mass(get_nucleus_code(A, Z));
} else {
m = get_mass(PID);
}
@@ -172,7 +172,7 @@ namespace corsika::nuclear_stack {
setNuclearZ(Z);
HEPMassType m = 0_GeV;
if (PID == Code::Nucleus) {
- m = get_nucleus_mass(A, Z);
+ m = get_nucleus_mass(get_nucleus_code(A, Z));
} else {
m = get_mass(PID);
}
@@ -240,7 +240,7 @@ namespace corsika::nuclear_stack {
inline HEPMassType NuclearParticleInterface<InnerParticleInterface,
StackIteratorInterface>::getMass() const {
if (super_type::getPID() == Code::Nucleus)
- return get_nucleus_mass(getNuclearA(), getNuclearZ());
+ return get_nucleus_mass(get_nucleus_code(getNuclearA(), getNuclearZ()));
return super_type::getMass();
}
diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp
index 85a65346b48b343e586516e2c8cbb1828db4ec25..8b9af126837c5d9c54942cb34b00c898283acd88 100644
--- a/corsika/framework/core/ParticleProperties.hpp
+++ b/corsika/framework/core/ParticleProperties.hpp
@@ -27,45 +27,49 @@
namespace corsika {
/**
- @defgroup Particles Particle Properties
-
- The properties of all particles are saved in static and flat
- arrays. There is a enum corsika::Code to identify each
- particles, and each individual particles has its own static class,
- which can be used to retrieve its physical properties.
-
- The properties of all elementary particles are accessible here. The data
- are taken from the Pythia ParticleData.xml file.
-
- Particle data can be accessed via global function in namespace corsika, or via
- static classes for each particle type. These classes all have the interface (example
- for the class corsika::Electron):
-
- @code{.cpp}
- static constexpr Code code{Code::Electron};
- static constexpr Code anti_code{Code::Positron};
- static constexpr HEPMassType mass{corsika::get_mass(code)};
- static constexpr ElectricChargeType charge{corsika::get_charge(code)};
- static constexpr int charge_number{corsika::get_charge_number(code)};
- static constexpr std::string_view name{corsika::get_name(code)};
- static constexpr bool is_nucleus{corsika::is_nucleus(code)};
- @endcode
-
- The names, relations and properties of all particles known to CORSIKA 8 are listed
- below.
-
- @addtogroup Particles
- @{
+ * @defgroup Particles Particle Properties
+ *
+ * The properties of all particles are saved in static and flat
+ * arrays. There is a enum corsika::Code to identify each
+ * particles, and each individual particles has its own static class,
+ * which can be used to retrieve its physical properties.
+ *
+ * The properties of all elementary particles are accessible here. The data
+ * are taken from the Pythia ParticleData.xml file.
+ *
+ * Particle data can be accessed via global function in namespace corsika, or via
+ * static classes for each particle type. These classes all have the interface (example
+ * for the class corsika::Electron):
+ *
+ * @code{.cpp}
+ * static constexpr Code code{Code::Electron};
+ * static constexpr Code anti_code{Code::Positron};
+ * static constexpr HEPMassType mass{corsika::get_mass(code)};
+ * static constexpr ElectricChargeType charge{corsika::get_charge(code)};
+ * static constexpr int charge_number{corsika::get_charge_number(code)};
+ * static constexpr std::string_view name{corsika::get_name(code)};
+ * static constexpr bool is_nucleus{corsika::is_nucleus(code)};
+ * @endcode
+ *
+ * The names, relations and properties of all particles known to CORSIKA 8 are listed
+ * below.
+ *
+ * @addtogroup Particles
+ * @{
*/
- /** The Code enum is the actual place to define CORSIKA 8 particle codes. */
- enum class Code : int16_t;
+ /**
+ * @brief The Code enum is the actual place to define CORSIKA 8 particle codes.
+ */
+ enum class Code : int32_t;
- /** Specifically for PDG ids */
+ /**
+ * @brief Specifically for PDG ids.
+ */
enum class PDGCode : int32_t;
using CodeIntType = std::underlying_type<Code>::type;
- using PDGCodeType = std::underlying_type<PDGCode>::type;
+ using PDGCodeIntType = std::underlying_type<PDGCode>::type;
// forward declarations to be used in GeneratedParticleProperties
int16_t constexpr get_charge_number(Code const); //!< electric charge in units of e
@@ -88,31 +92,81 @@ namespace corsika {
//! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme"
PDGCode constexpr get_PDG(Code const);
- PDGCode constexpr get_PDG(unsigned int A, unsigned int Z);
std::string_view constexpr get_name(Code const); //!< name of the particle as string
TimeType constexpr get_lifetime(Code const); //!< lifetime
- //! true iff the particle is a hard-coded nucleus or Code::Nucleus
- bool constexpr is_nucleus(Code const);
bool constexpr is_hadron(Code const); //!< true iff particle is hadron
bool constexpr is_em(Code const); //!< true iff particle is electron, positron or photon
bool constexpr is_muon(Code const); //!< true iff particle is mu+ or mu-
bool constexpr is_neutrino(Code const); //!< true iff particle is (anti-) neutrino
- int constexpr get_nucleus_A(
+
+ /**
+ * @brief Creates the Code for a nucleus of type 10LZZZAAAI.
+ *
+ * @return internal nucleus Code
+ */
+ Code constexpr get_nucleus_code(unsigned int const A, unsigned int const Z);
+ /**
+ * @brief Checks if Code corresponds to a nucleus.
+ *
+ * @return true if nucleus.
+ * @return false if not nucleus.
+ */
+ bool constexpr is_nucleus(Code const);
+
+ /**
+ * @brief Get the mass number A for nucleus.
+ *
+ * @return int size of nucleus.
+ */
+ unsigned int constexpr get_nucleus_A(
Code const); //!< returns A for hard-coded nucleus, otherwise 0
- int constexpr get_nucleus_Z(
+
+ /**
+ * @brief Get the charge number Z for nucleus.
+ *
+ * @return int charge of nucleus.
+ */
+ unsigned int constexpr get_nucleus_Z(
Code const); //!< returns Z for hard-coded nucleus, otherwise 0
- //! returns mass of (A,Z) nucleus, disregarding binding energy
- HEPMassType get_nucleus_mass(unsigned int const, unsigned int const);
+ /**
+ * @brief Calculates the mass of nucleus.
+ *
+ * @return HEPMassType the mass of (A,Z) nucleus, disregarding binding energy.
+ */
+ HEPMassType constexpr get_nucleus_mass(Code const code);
+
+ /**
+ * @brief Get the nucleus name.
+ *
+ * @param code
+ * @return std::string
+ */
+ inline std::string_view get_nucleus_name(Code const code);
- //! convert PDG code to CORSIKA 8 internal code
+ /**
+ * @brief convert PDG code to CORSIKA 8 internal code.
+ *
+ * @return Code internal code.
+ */
Code convert_from_PDG(PDGCode const);
+ /**
+ * @brief Returns list of all non-nuclei particles.
+ *
+ * @return std::initializer_list<Code> constexpr
+ */
std::initializer_list<Code> constexpr get_all_particles();
- //! the output stream operator for human-readable particle codes
+ /**
+ * @brief Code output operator.
+ *
+ * The output stream operator for human-readable particle codes.
+ *
+ * @return std::ostream&
+ */
std::ostream& operator<<(std::ostream&, corsika::Code);
/** @}*/
diff --git a/corsika/framework/process/InteractionHistogram.hpp b/corsika/framework/process/InteractionHistogram.hpp
index 78796f609e3a83316a52d818cf045425af7c507a..b687153b58e96bf4fb4235cb431318424c9136af 100644
--- a/corsika/framework/process/InteractionHistogram.hpp
+++ b/corsika/framework/process/InteractionHistogram.hpp
@@ -57,11 +57,11 @@ namespace corsika {
@param projectile_id corsika::Code of particle
@param lab_energy Energy in lab. frame
@param mass_target Mass of target particle
- @param A if projectile_id is corsika::Nucleus : Mass of nucleus
- @param Z if projectile_id is corsika::Nucleus : Charge of nucleus
+ @param A if projectile_id is Nucleus : Mass of nucleus
+ @param Z if projectile_id is Nucleus : Charge of nucleus
*/
- void fill(Code projectile_id, HEPEnergyType lab_energy, HEPEnergyType mass_target,
- int A = 0, int Z = 0);
+ void fill(Code const projectile_id, HEPEnergyType const lab_energy,
+ HEPEnergyType const mass_target);
//! return histogram in c.m.s. frame
hist_type const& CMSHist() const { return inthist_cms_; }
diff --git a/corsika/modules/epos/Interaction.hpp b/corsika/modules/epos/Interaction.hpp
index efaac0705f33174fa1910c780791e59651eaf1c1..56900ca01cdf99f41f6dc38b7cc9578bc6e7afad 100644
--- a/corsika/modules/epos/Interaction.hpp
+++ b/corsika/modules/epos/Interaction.hpp
@@ -77,8 +77,8 @@ namespace corsika::epos {
std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_epos_Interaction");
HEPEnergyType const minEnergyCoM_ = 6 * 1e9 * electronvolt;
HEPEnergyType const maxEnergyCoM_ = 2.e6 * 1e9 * electronvolt;
- int const maxTargetMassNumber_ = 20;
- int const minNuclearTargetA_ = 4;
+ static unsigned int constexpr maxTargetMassNumber_ = 20;
+ static unsigned int constexpr minNuclearTargetA_ = 4;
};
} // namespace corsika::epos
diff --git a/corsika/modules/epos/ParticleConversion.hpp b/corsika/modules/epos/ParticleConversion.hpp
index 3ee82a72168acf9479ae7a4e2e4d6a479def7fc1..736b93074d3e50d42c110a249283f7a31f779443 100644
--- a/corsika/modules/epos/ParticleConversion.hpp
+++ b/corsika/modules/epos/ParticleConversion.hpp
@@ -21,7 +21,7 @@ namespace corsika::epos {
using EposCodeIntType = std::underlying_type<EposCode>::type;
/**
- These are the possible projectile for which Epos knows the cross section
+ * These are the possible projectile for which Epos knows the cross section.
*/
enum class EposXSClass : int8_t {
CannotInteract = 0,
@@ -34,14 +34,27 @@ namespace corsika::epos {
#include <corsika/modules/epos/Generated.inc>
- EposCode constexpr convertToEpos(Code pCode) {
- return corsika2epos[static_cast<CodeIntType>(pCode)];
+ unsigned int constexpr get_nucleus_A(EposCode const eposId) {
+ // 100ZZZAAA0 -> std. pdg code
+ EposCodeIntType const eposPdg = static_cast<EposCodeIntType>(eposId);
+ return int(abs(eposPdg) / 10) % 1000;
+ }
+ unsigned int constexpr get_nucleus_Z(EposCode const eposId) {
+ // 100ZZZAAA0 -> std. pdg code
+ EposCodeIntType const eposPdg = static_cast<EposCodeIntType>(eposId);
+ return int(abs(eposPdg) / 10000) % 1000;
+ }
+
+ EposCode constexpr convertToEpos(Code const code) {
+ return corsika2epos[static_cast<CodeIntType>(code)];
}
- Code constexpr convertFromEpos(EposCode pCode) {
- EposCodeIntType const s = static_cast<EposCodeIntType>(pCode);
+ Code constexpr convertFromEpos(EposCode const eposId) {
+ EposCodeIntType const s = static_cast<EposCodeIntType>(eposId);
// if nucleus (pdg-id)
- if (s >= 1000000000) { return Code::Nucleus; }
+ if (s >= 1000000000) {
+ return get_nucleus_code(get_nucleus_A(eposId), get_nucleus_Z(eposId));
+ }
auto const corsikaCode = epos2corsika[s - minEpos];
if (corsikaCode == Code::Unknown) {
throw std::runtime_error(std::string("EPOS/CORSIKA conversion of ")
@@ -51,16 +64,19 @@ namespace corsika::epos {
return corsikaCode;
}
- int constexpr convertToEposRaw(Code pCode) {
- return static_cast<int>(convertToEpos(pCode));
+ int constexpr convertToEposRaw(Code const code) {
+ return static_cast<int>(convertToEpos(code));
}
- int constexpr getEposXSCode(Code pCode) {
+ int constexpr getEposXSCode(Code const code) {
+ if (is_nucleus(code)) { return static_cast<EposXSClassIntType>(EposXSClass::Baryon); }
return static_cast<EposXSClassIntType>(
- corsika2eposXStype[static_cast<CodeIntType>(pCode)]);
+ corsika2eposXStype[static_cast<CodeIntType>(code)]);
}
- bool constexpr canInteract(Code pCode) { return getEposXSCode(pCode) > 0; }
+ bool constexpr canInteract(Code const code) {
+ return is_nucleus(code) || getEposXSCode(code) > 0;
+ }
HEPMassType getEposMass(Code const);
diff --git a/corsika/modules/qgsjetII/Interaction.hpp b/corsika/modules/qgsjetII/Interaction.hpp
index 7af63723f5670fa0de59aea92c98f232b145d8d4..0030de2fccc3438377c11a801e0206e62d55a0c0 100644
--- a/corsika/modules/qgsjetII/Interaction.hpp
+++ b/corsika/modules/qgsjetII/Interaction.hpp
@@ -29,7 +29,7 @@ namespace corsika::qgsjetII {
~Interaction();
bool wasInitialized() { return initialized_; }
- int getMaxTargetMassNumber() const { return maxMassNumber_; }
+ unsigned int getMaxTargetMassNumber() const { return maxMassNumber_; }
bool isValidTarget(corsika::Code TargetId) const {
return is_nucleus(TargetId) && (get_nucleus_A(TargetId) < maxMassNumber_);
}
@@ -57,7 +57,7 @@ namespace corsika::qgsjetII {
corsika::default_prng_type& rng_ =
corsika::RNGManager<>::getInstance().getRandomStream("qgsjet");
- const int maxMassNumber_ = 208;
+ static unsigned int constexpr maxMassNumber_ = 208;
};
} // namespace corsika::qgsjetII
diff --git a/corsika/modules/qgsjetII/ParticleConversion.hpp b/corsika/modules/qgsjetII/ParticleConversion.hpp
index 473365f126a54e2249198d6ad1efc912363a7e96..3b4f689b8b7b1af73367f5f6a9355c273206d657 100644
--- a/corsika/modules/qgsjetII/ParticleConversion.hpp
+++ b/corsika/modules/qgsjetII/ParticleConversion.hpp
@@ -62,37 +62,37 @@ namespace corsika::qgsjetII {
return corsika2qgsjetII[static_cast<CodeIntType>(pCode)];
}
- Code constexpr convertFromQgsjetII(QgsjetIICode pCode) {
- auto const pCodeInt = static_cast<QgsjetIICodeIntType>(pCode);
- auto const corsikaCode = qgsjetII2corsika[pCodeInt - minQgsjetII];
+ Code constexpr convertFromQgsjetII(QgsjetIICode const code) {
+ auto const codeInt = static_cast<QgsjetIICodeIntType>(code);
+ auto const corsikaCode = qgsjetII2corsika[codeInt - minQgsjetII];
if (corsikaCode == Code::Unknown) {
throw std::runtime_error(std::string("QGSJETII/CORSIKA conversion of pCodeInt=")
- .append(std::to_string(pCodeInt))
+ .append(std::to_string(codeInt))
.append(" impossible"));
}
return corsikaCode;
}
- QgsjetIICodeIntType constexpr convertToQgsjetIIRaw(Code pCode) {
- return static_cast<QgsjetIICodeIntType>(convertToQgsjetII(pCode));
+ QgsjetIICodeIntType constexpr convertToQgsjetIIRaw(Code const code) {
+ return static_cast<QgsjetIICodeIntType>(convertToQgsjetII(code));
}
- QgsjetIIXSClass constexpr getQgsjetIIXSCode(Code pCode) {
- // if (pCode == corsika::particles::Code::Nucleus)
- // static_cast(QgsjetIIXSClassIntType>();
- return corsika2qgsjetIIXStype[static_cast<CodeIntType>(pCode)];
+ QgsjetIIXSClass constexpr getQgsjetIIXSCode(Code const code) {
+ return corsika2qgsjetIIXStype[static_cast<CodeIntType>(code)];
}
- QgsjetIIXSClassIntType constexpr getQgsjetIIXSCodeRaw(Code pCode) {
- return static_cast<QgsjetIIXSClassIntType>(getQgsjetIIXSCode(pCode));
+ QgsjetIIXSClassIntType constexpr getQgsjetIIXSCodeRaw(Code const code) {
+ return static_cast<QgsjetIIXSClassIntType>(getQgsjetIIXSCode(code));
}
- bool constexpr canInteract(Code pCode) {
- return getQgsjetIIXSCode(pCode) != QgsjetIIXSClass::CannotInteract;
+ bool constexpr canInteract(Code const code) {
+ if (is_nucleus(code)) return true;
+ return getQgsjetIIXSCode(code) != QgsjetIIXSClass::CannotInteract;
}
- QgsjetIIHadronType constexpr getQgsjetIIHadronType(Code pCode) {
- return corsika2qgsjetIIHadronType[static_cast<CodeIntType>(pCode)];
+ QgsjetIIHadronType constexpr getQgsjetIIHadronType(Code const code) {
+ if (is_nucleus(code)) return QgsjetIIHadronType::NucleusType;
+ return corsika2qgsjetIIHadronType[static_cast<CodeIntType>(code)];
}
} // namespace corsika::qgsjetII
diff --git a/corsika/modules/sibyll/Interaction.hpp b/corsika/modules/sibyll/Interaction.hpp
index 423d9319b9c9458265d65f1a28f0a0f1a17240f2..b53fb63080b4a0a44b751a9abf9e73f60b940401 100644
--- a/corsika/modules/sibyll/Interaction.hpp
+++ b/corsika/modules/sibyll/Interaction.hpp
@@ -52,15 +52,15 @@ namespace corsika::sibyll {
void doInteraction(TSecondaries&);
private:
- int getMaxTargetMassNumber() const { return maxTargetMassNumber_; }
+ unsigned int constexpr getMaxTargetMassNumber() const { return maxTargetMassNumber_; }
HEPEnergyType getMinEnergyCoM() const { return minEnergyCoM_; }
HEPEnergyType getMaxEnergyCoM() const { return maxEnergyCoM_; }
default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sibyll");
const HEPEnergyType minEnergyCoM_ = 10. * 1e9 * electronvolt;
const HEPEnergyType maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
- const int maxTargetMassNumber_ = 18;
- const int minNuclearTargetA_ = 4;
+ static unsigned int constexpr maxTargetMassNumber_ = 18;
+ static unsigned int constexpr minNuclearTargetA_ = 4;
// data members
int count_ = 0;
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index 09d997c5d7d343165874964418063fb6ce211d8a..f7fc083e86337ce34c94539fdb41b65219707dc1 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -397,11 +397,7 @@ int main(int argc, char** argv) {
stack.clear();
// add the desired particle to the stack
- if (A > 1) {
- stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
- } else {
- stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
- }
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
// if we want to fix the first location of the shower
if (force_interaction) {
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 6800a84be8858d31ed01d36d213ef56754a265b1..4b10e0436a04891b727c0310fec7a4b4ca950147 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -131,7 +131,7 @@ int main(int argc, char** argv) {
50_uT, 0_T});
builder.setNuclearComposition(
{{Code::Nitrogen, Code::Oxygen},
- {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+ {0.7847, 1. - 0.7847}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
diff --git a/src/framework/core/CMakeLists.txt b/src/framework/core/CMakeLists.txt
index 097908cf269ed53aa4465c0f920455f8ae8215a2..4faa8829e543b17d6905b8dbcd4c957bfff4761f 100644
--- a/src/framework/core/CMakeLists.txt
+++ b/src/framework/core/CMakeLists.txt
@@ -7,10 +7,10 @@ add_custom_command (
OUTPUT ${output_dir}/GeneratedParticleProperties.inc
${output_dir}/GeneratedParticleClasses.inc
${output_dir}/particle_db.pkl
- COMMAND ${input_dir}/pdxml_reader.py ${input_dir}/ParticleData.xml
- ${input_dir}/NuclearData.xml
- ${input_dir}/ParticleClassNames.xml
- DEPENDS ${input_dir}/pdxml_reader.py
+ COMMAND ${input_dir}/code_generator.py ${input_dir}/ParticleData.xml
+ ${input_dir}/NuclearData.xml
+ ${input_dir}/ParticleClassNames.xml
+ DEPENDS ${input_dir}/code_generator.py
${input_dir}/ParticleData.xml
${input_dir}/NuclearData.xml
${input_dir}/ParticleClassNames.xml
diff --git a/src/framework/core/NuclearData.xml b/src/framework/core/NuclearData.xml
index 73eda0ef2350aa6849df60c61a34231dd9fff869..257b6f260a43ddad22bdd6be5e99a545f6dd1cb4 100644
--- a/src/framework/core/NuclearData.xml
+++ b/src/framework/core/NuclearData.xml
@@ -1,8 +1,5 @@
<chapter name="Nuclear Data">
-<particle id="1000000000" name="nucleus" A="0" Z="0" >
-</particle>
-
<particle id="1000010010" name="hydrogen" A="1" Z="1" >
</particle>
@@ -60,10 +57,22 @@
<particle id="1000541280" name="xenon" A="128" Z="54" >
</particle>
+<particle id="1000741840" name="tungsten" A="184" Z="74" >
+</particle>
+
+<particle id="1000791970" name="gold" A="197" Z="79" >
+</particle>
+
<particle id="1000822080" name="lead" A="208" Z="82" >
</particle>
<particle id="1000862220" name="radon" A="222" Z="86" >
</particle>
+<particle id="1000922380" name="uranium" A="238" Z="92" >
+</particle>
+
+<particle id="1001243080" name="ralfium" A="308" Z="124" >
+</particle>
+
</chapter>
diff --git a/src/framework/core/pdxml_reader.py b/src/framework/core/code_generator.py
similarity index 84%
rename from src/framework/core/pdxml_reader.py
rename to src/framework/core/code_generator.py
index 13a79f33482d77e2e037217bf9ebd0e05a0b45bb..e06488c7b4da21840a26986e4514741d17a6c197 100755
--- a/src/framework/core/pdxml_reader.py
+++ b/src/framework/core/code_generator.py
@@ -19,6 +19,11 @@ mproton = 0.9382720813 # GeV
namespace = "corsika"
+# IDs of Nuclei are 10LZZZAAAI
+nucleusIdStr = "10{L:01d}{Z:03d}{A:03d}{I:01d}" # used with .format(L=,Z=,A=,I=)
+nucleusIdOffset = int(nucleusIdStr.format(L=0,A=0,Z=0,I=0))
+
+
##############################################################
#
# reading xml input data, return line by line particle data
@@ -238,7 +243,7 @@ def read_nuclei_db(filename, particle_db, classnames):
"mass": mass, # in GeV
"electric_charge": electric_charge, # in e/3
"lifetime": lifetime,
- "ngc_code": next(counter),
+ "ngc_code": int(nucleusIdStr.format(L=0,A=A,Z=Z,I=0)),
"A": A,
"Z": Z,
"isNucleus": True,
@@ -297,22 +302,37 @@ def gen_conversion_PDG_ngc(particle_db):
#
# return string with enum of all internal particle codes
#
-def gen_internal_enum(particle_db):
+def gen_internal_enum(particle_db, nuclei_db):
string = ("//! @cond EXCLUDE_DOXY\n"
"enum class Code : CodeIntType {\n"
" FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops...
- for k in filter(lambda k: "ngc_code" in particle_db[k], particle_db):
+ # non-nuclei loop
+ for k in filter(lambda k: "ngc_code" in particle_db[k] and not particle_db[k]["isNucleus"],
+ particle_db):
last_ngc_id = particle_db[k]['ngc_code']
string += " {key:s} = {code:d},\n".format(key=k, code=last_ngc_id)
- string += (" LastParticle = {:d},\n" # identifier for eventual loops...
- "}}; //! @endcond").format(last_ngc_id + 1)
+ string += " LastParticle = {:d},\n".format(last_ngc_id + 1) # identifier for eventual loops...
- if last_ngc_id > 0x7fff: # does not fit into int16_t
+ if last_ngc_id > 0x7fffffff: # does not fit into int32_t, this will never happen....
raise Exception(
"Integer overflow in internal particle code definition prevented!")
+ if last_ngc_id + 1 >= nucleusIdOffset:
+ raise Exception(
+ "Too many particles. Fix conflict with Code::Nuclear == {:d} (just increase id...) !").format(nucleusIdOffset)
+
+ # marker to mark the beginning of generic nuclear IDs: 1aaazzz
+ string += " Nucleus = {:d},\n".format(nucleusIdOffset) # identifier for Nuclei
+
+ # nuclei loop
+ for k in filter(lambda k: "ngc_code" in nuclei_db[k] and nuclei_db[k]["isNucleus"],
+ nuclei_db):
+ last_ngc_id = nuclei_db[k]['ngc_code']
+ string += " {key:s} = {code:d},\n".format(key=k, code=last_ngc_id)
+
+ string += "}; //! @endcond"
return string
@@ -322,7 +342,7 @@ def gen_internal_enum(particle_db):
#
def gen_pdg_enum(particle_db):
string = ("//! @cond EXCLUDE_DOXY\n"
- "enum class PDGCode : PDGCodeType {\n")
+ "enum class PDGCode : PDGCodeIntType {\n")
for cId in particle_db:
pdgCode = particle_db[cId]['pdg']
@@ -363,6 +383,7 @@ def gen_properties(particle_db):
for k in particle_db:
string += " 0 * corsika::units::si::electronvolt, // {name:s}\n".format( name = k)
string += "};\n\n"
+ string += "static corsika::units::si::HEPEnergyType threshold_nuclei = 0_eV;\n"
# PDG code table
string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n"
@@ -377,7 +398,7 @@ def gen_properties(particle_db):
string += "};\n"
# electric charges table
- string += "static constexpr std::array<int16_t, size> electric_charges = {\n"
+ string += "static constexpr std::array<int32_t, size> electric_charges = {\n"
for p in particle_db.values():
string += " {charge:d},\n".format(charge=p['electric_charge'] // 3)
string += "};\n"
@@ -409,22 +430,6 @@ def gen_properties(particle_db):
string += " {val},\n".format(val=value)
string += "};\n"
- ### nuclear data ###
-
- # nucleus mass number A
- string += "static constexpr std::array<int16_t, size> nucleusA = {\n"
- for p in particle_db.values():
- A = p.get('A', 0)
- string += " {val},\n".format(val=A)
- string += "};\n"
-
- # nucleus charge number Z
- string += "static constexpr std::array<int16_t, size> nucleusZ = {\n"
- for p in particle_db.values():
- Z = p.get('Z', 0)
- string += " {val},\n".format(val=Z)
- string += "};\n"
-
return string
@@ -432,35 +437,36 @@ def gen_properties(particle_db):
#
# return string with a list of classes for all particles
#
-def gen_classes(particle_db):
+def gen_classes(particle_db, nuclei_db):
string = ("// list of C++ classes to access particle properties\n"
"/** @defgroup ParticleClasses \n"
" @{ */\n")
- for cname in particle_db:
+ common_db = OrderedDict(list(particle_db.items()) + list(nuclei_db.items()))
+ for cname in common_db:
if cname == "Nucleus":
string += "// skipping Nucleus"
continue
antiP = 'Unknown'
- for cname_anti in particle_db:
- if (particle_db[cname_anti]['name'] == particle_db[cname]['antiName']):
+ for cname_anti in common_db:
+ if (common_db[cname_anti]['name'] == common_db[cname]['antiName']):
antiP = cname_anti
break
string += "\n"
string += "/** @class " + cname + "\n\n"
string += " * Particle properties are taken from the PYTHIA8 ParticleData.xml file:<br>\n"
- string += " * - pdg=" + str(particle_db[cname]['pdg']) + "\n"
- string += " * - mass=" + str(particle_db[cname]['mass']) + " GeV \n"
+ string += " * - pdg=" + str(common_db[cname]['pdg']) + "\n"
+ string += " * - mass=" + str(common_db[cname]['mass']) + " GeV \n"
string += " * - charge= " + \
- str(particle_db[cname]['electric_charge'] // 3) + " \n"
+ str(common_db[cname]['electric_charge'] // 3) + " \n"
string += " * - name=" + str(cname) + "\n"
string += " * - anti=" + str(antiP) + "\n"
- if (particle_db[cname]['isNucleus']):
- string += " * - nuclear A=" + str(particle_db[cname]['A']) + "\n"
- string += " * - nuclear Z=" + str(particle_db[cname]['Z']) + "\n"
+ if (common_db[cname]['isNucleus']):
+ string += " * - nuclear A=" + str(common_db[cname]['A']) + "\n"
+ string += " * - nuclear Z=" + str(common_db[cname]['Z']) + "\n"
string += "*/\n\n"
string += "class " + cname + " {\n"
string += " /** @cond EXCLUDE_DOXY */ \n"
@@ -473,7 +479,7 @@ def gen_classes(particle_db):
string += " static constexpr int charge_number{corsika::get_charge_number(code)};\n"
string += " static constexpr std::string_view name{corsika::get_name(code)};\n"
string += " static constexpr bool is_nucleus{corsika::is_nucleus(code)};\n"
- if particle_db[cname]['isNucleus']:
+ if common_db[cname]['isNucleus']:
string += " static constexpr int nucleus_A{corsika::get_nucleus_A(code)};\n"
string += " static constexpr int nucleus_Z{corsika::get_nucleus_Z(code)};\n"
string += " private:\n"
@@ -490,7 +496,7 @@ def gen_classes(particle_db):
#
#
def inc_start():
- string = ('// generated by pdxml_reader.py\n'
+ string = ('// generated by framework/core/code_generator.py\n'
'// MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN. \n'
'\n'
'namespace corsika {\n'
@@ -531,8 +537,9 @@ def inc_end():
#
# Serialize particle_db into file
#
-def serialize_particle_db(particle_db, file):
- pickle.dump(particle_db, file)
+def serialize_particle_db(particle_db, nuclei_db, file):
+ common_db = OrderedDict(list(particle_db.items()) + list(nuclei_db.items()))
+ pickle.dump(common_db, file)
###################################################################
@@ -546,16 +553,17 @@ if __name__ == "__main__":
sys.argv[0]), file=sys.stderr)
sys.exit(1)
- print("\n pdxml_reader.py: automatically produce particle properties from input files\n")
+ print("\n code_generator.py: automatically produce particle properties from input files\n")
- names = read_class_names(sys.argv[3])
- particle_db = OrderedDict()
- read_pythia_db(sys.argv[1], particle_db, names)
- read_nuclei_db(sys.argv[2], particle_db, names)
+ names = read_class_names(sys.argv[3]) # re-names and conventions
+ particle_db = OrderedDict() # the DB for pythia8 pdg particles
+ read_pythia_db(sys.argv[1], particle_db, names) # pythia8 pdg DB
+ nuclei_db = OrderedDict() # the DB for specific nuclei
+ read_nuclei_db(sys.argv[2], nuclei_db, names) # list of nuclei
with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f)
- print(gen_internal_enum(particle_db), file=f)
+ print(gen_internal_enum(particle_db, nuclei_db), file=f)
print(gen_pdg_enum(particle_db), file=f)
print(detail_start(), file=f)
print(gen_properties(particle_db), file=f)
@@ -565,8 +573,8 @@ if __name__ == "__main__":
with open("GeneratedParticleClasses.inc", "w") as f:
print(inc_start(), file=f)
- print(gen_classes(particle_db), file=f)
+ print(gen_classes(particle_db, nuclei_db), file=f)
print(inc_end(), file=f)
with open("particle_db.pkl", "wb") as f:
- serialize_particle_db(particle_db, f)
+ serialize_particle_db(particle_db, nuclei_db, f)
diff --git a/src/modules/epos/code_generator.py b/src/modules/epos/code_generator.py
index f47efbb965fafc025e6b5a43c9e11e99178d87e3..c3de621e9bff0f089dd550069e60a8a584bcec5f 100755
--- a/src/modules/epos/code_generator.py
+++ b/src/modules/epos/code_generator.py
@@ -39,7 +39,7 @@ def read_epos_codes(filename, particle_db):
particle_db[identifier]["epos_code"] = int(epos_code)
particle_db[identifier]["epos_xsType"] = xsType
except KeyError as e:
- raise Exception("Identifier '{:s}' not found in particle_db".format(identifier))
+ raise Exception("Identifier '{:s}' not found in CORSIKA8 particle_db".format(identifier))
def set_default_epos_definition(particle_db):
@@ -49,7 +49,7 @@ def set_default_epos_definition(particle_db):
This is achieved here.
- The function return nothing, but modified the input particle_db by adding the
+ The function returns nothing, but modified the input particle_db by adding the
fields 'xsType' and 'hadronType'
'''
for identifier, pData in particle_db.items():
@@ -83,6 +83,7 @@ def generate_corsika2epos(particle_db):
'''
string = "std::array<EposCode, {:d}> constexpr corsika2epos = {{\n".format(len(particle_db))
for identifier, pData in particle_db.items():
+ if pData['isNucleus']: continue
if 'epos_code' in pData:
string += " EposCode::{:s}, \n".format(identifier)
else:
@@ -98,6 +99,7 @@ def generate_corsika2epos_xsType(particle_db):
'''
string = "std::array<EposXSClass, {:d}> constexpr corsika2eposXStype = {{\n".format(len(particle_db))
for identifier, pData in particle_db.items():
+ if pData['isNucleus']: continue
if 'epos_xsType' in pData:
string += " EposXSClass::{:s}, // {:s}\n".format(pData['epos_xsType'], identifier)
else:
@@ -116,7 +118,6 @@ def generate_epos2corsika(particle_db) :
for identifier, pData in particle_db.items() :
if 'epos_code' in pData:
minID = min(minID, pData['epos_code'])
-
string += "EposCodeIntType constexpr minEpos = {:d};\n\n".format(minID)
pDict = {}
diff --git a/src/modules/epos/epos_codes.dat b/src/modules/epos/epos_codes.dat
index 9c942914e92a10fbd3bbb8b42750b66d315d8ce7..70b6cb380bf1e79651a426e7e5eae36aebbf0104 100644
--- a/src/modules/epos/epos_codes.dat
+++ b/src/modules/epos/epos_codes.dat
@@ -7,7 +7,7 @@
Unknown 0 0 CannotInteract
# Here is the list of particles known to EPOS
-Photon 10 0 CannotInteract
+Photon 10 0 CannotInteract
Electron 12 0 CannotInteract
Positron -12 0 CannotInteract
MuMinus 14 0 CannotInteract
diff --git a/src/modules/qgsjetII/code_generator.py b/src/modules/qgsjetII/code_generator.py
index 8321aad19ab1ae30252b930800a62aa7f0fd7b6b..b867e91dfeeb86a3bb44737d55ce5872f223b8cb 100755
--- a/src/modules/qgsjetII/code_generator.py
+++ b/src/modules/qgsjetII/code_generator.py
@@ -120,6 +120,7 @@ def generate_corsika2qgsjetII(particle_db):
'''
string = "std::array<QgsjetIICode, {:d}> constexpr corsika2qgsjetII = {{\n".format(len(particle_db))
for identifier, pData in particle_db.items():
+ if pData['isNucleus']: continue
if 'qgsjetII_code' in pData:
string += " QgsjetIICode::{:s}, \n".format(identifier)
else:
@@ -134,6 +135,7 @@ def generate_corsika2qgsjetII_xsType(particle_db):
'''
string = "std::array<QgsjetIIXSClass, {:d}> constexpr corsika2qgsjetIIXStype = {{\n".format(len(particle_db))
for identifier, pData in particle_db.items():
+ if pData['isNucleus']: continue
modelCodeXS = pData.get("qgsjetII_xsType", "CannotInteract")
string += " QgsjetIIXSClass::{:s}, // {:s}\n".format(modelCodeXS, identifier if modelCodeXS else identifier + " (not implemented in QGSJETII)")
string += "};\n"
@@ -146,6 +148,7 @@ def generate_corsika2qgsjetII_hadronType(particle_db):
'''
string = "std::array<QgsjetIIHadronType, {:d}> constexpr corsika2qgsjetIIHadronType = {{\n".format(len(particle_db))
for identifier, pData in particle_db.items():
+ if pData['isNucleus']: continue
modelCode = pData.get("qgsjetII_hadronType", "UndefinedType")
string += " QgsjetIIHadronType::{:s}, // {:s}\n".format(modelCode, identifier if modelCode else identifier + " (not implemented in QGSJETII)")
string += "};\n"
diff --git a/src/modules/sibyll/code_generator.py b/src/modules/sibyll/code_generator.py
index 22cb11022741597b6807c844766dcf0158f24365..20407f2c7e7cfe1808ab2c92b9c9ddf081feee7b 100755
--- a/src/modules/sibyll/code_generator.py
+++ b/src/modules/sibyll/code_generator.py
@@ -63,6 +63,7 @@ def generate_corsika2sibyll(particle_db):
'''
string = "std::array<SibyllCode, {:d}> constexpr corsika2sibyll = {{\n".format(len(particle_db))
for identifier, pData in particle_db.items():
+ if pData['isNucleus']: continue
if 'sibyll_code' in pData:
string += " SibyllCode::{:s}, \n".format(identifier)
else:
@@ -78,6 +79,7 @@ def generate_corsika2sibyll_xsType(particle_db):
'''
string = "std::array<SibyllXSClass, {:d}> constexpr corsika2sibyllXStype = {{\n".format(len(particle_db))
for identifier, pData in particle_db.items():
+ if pData['isNucleus']: continue
if 'sibyll_xsType' in pData:
string += " SibyllXSClass::{:s}, // {:s}\n".format(pData['sibyll_xsType'], identifier)
else:
diff --git a/tests/common/SetupTestStack.hpp b/tests/common/SetupTestStack.hpp
index d41f1aa2726f65a7eeb7a31b3132e38ad51e1be8..a5bfb3d23622a3da7f0a1ea4102f514735005906 100644
--- a/tests/common/SetupTestStack.hpp
+++ b/tests/common/SetupTestStack.hpp
@@ -35,7 +35,7 @@ namespace corsika::setup::testing {
**/
inline std::tuple<std::unique_ptr<setup::Stack>, std::unique_ptr<setup::StackView>>
- setup_stack(Code vProjectileType, int vA, int vZ, HEPEnergyType vMomentum,
+ setup_stack(Code vProjectileType, HEPEnergyType vMomentum,
setup::Environment::BaseNodeType* const vNodePtr,
CoordinateSystemPtr const& cs) {
@@ -44,19 +44,11 @@ namespace corsika::setup::testing {
Point const origin(cs, {0_m, 0_m, 0_m});
MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
- if (vProjectileType == Code::Nucleus) {
- 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
- auto particle =
- 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));
- }
+ auto particle =
+ 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));
}
} // namespace corsika::setup::testing
diff --git a/tests/framework/testCascade.hpp b/tests/framework/testCascade.hpp
index 638424f893c9fa1fef7dd957b1d7fc0010285cc2..cb6e0b5a0c8e06ab7a51fbf6c5a2565760171e97 100644
--- a/tests/framework/testCascade.hpp
+++ b/tests/framework/testCascade.hpp
@@ -14,7 +14,7 @@
#include <corsika/framework/stack/CombinedStack.hpp>
#include <corsika/framework/stack/SecondaryView.hpp>
#include <corsika/stack/GeometryNodeStackExtension.hpp>
-#include <corsika/stack/NuclearStackExtension.hpp>
+#include <corsika/stack/VectorStack.hpp>
using TestEnvironmentInterface = corsika::IEmpty;
using TestEnvironmentType = corsika::Environment<TestEnvironmentInterface>;
@@ -26,12 +26,13 @@ using SetupGeometryDataInterface =
// combine particle data stack with geometry information for tracking
template <typename StackIter>
using StackWithGeometryInterface =
- corsika::CombinedParticleInterface<corsika::nuclear_stack::ParticleDataStack::pi_type,
+ corsika::CombinedParticleInterface<corsika::VectorStack::pi_type,
SetupGeometryDataInterface, StackIter>;
-using TestCascadeStack = corsika::CombinedStack<
- typename corsika::nuclear_stack::ParticleDataStack::stack_data_type,
- corsika::node::GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
+using TestCascadeStack =
+ corsika::CombinedStack<typename corsika::VectorStack::stack_data_type,
+ corsika::node::GeometryData<TestEnvironmentType>,
+ StackWithGeometryInterface>;
/*
See also Issue 161
diff --git a/tests/framework/testInteractionCounter.cpp b/tests/framework/testInteractionCounter.cpp
index 4b4b7ff25dca0ba6a5cf44018c9efde66a150e82..4ea5d49b04f462e43517ed96b3bf6f565a29fba3 100644
--- a/tests/framework/testInteractionCounter.cpp
+++ b/tests/framework/testInteractionCounter.cpp
@@ -33,12 +33,11 @@ const std::string refDataDir = std::string(REFDATADIR); // from cmake
struct DummyProcess {
template <typename TParticle>
- GrammageType getInteractionLength([[maybe_unused]] TParticle const& particle) {
+ GrammageType getInteractionLength(TParticle const&) {
return 100_g / 1_cm / 1_cm;
}
-
template <typename TParticle>
- void doInteraction([[maybe_unused]] TParticle& projectile) {}
+ void doInteraction(TParticle&) {}
};
TEST_CASE("InteractionCounter", "[process]") {
@@ -58,7 +57,7 @@ TEST_CASE("InteractionCounter", "[process]") {
SECTION("DoInteraction nucleus") {
unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Nucleus, A, Z, 105_TeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ get_nucleus_code(A, Z), 105_TeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
CHECK(stackPtr->getEntries() == 1);
CHECK(secViewPtr->getEntries() == 0);
@@ -103,7 +102,7 @@ TEST_CASE("InteractionCounter", "[process]") {
SECTION("DoInteraction Lambda") {
auto constexpr code = Code::Lambda0;
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- code, 0, 0, 105_TeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ code, 105_TeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
CHECK(stackPtr->getEntries() == 1);
CHECK(secViewPtr->getEntries() == 0);
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index 322c553c7a45358a61e3e389fafb888a15b7d59d..f12faf8039a3baa7031416e49bfe5e9bc4c8d467 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -47,7 +47,7 @@ TEST_CASE("ParticleProperties", "[Particles]") {
CHECK(Electron::name == "e-");
CHECK(get_name(Code::Electron) == "e-");
CHECK(PiMinus::name == "pi-");
- CHECK(Iron::name == "iron");
+ CHECK(Iron::name == "nucleus");
CHECK(Photon::name == "photon");
}
@@ -112,7 +112,6 @@ TEST_CASE("ParticleProperties", "[Particles]") {
CHECK(is_hadron(Code::Proton));
CHECK(is_hadron(Code::PiPlus));
CHECK(is_hadron(Code::Oxygen));
- CHECK(is_hadron(Code::Nucleus));
}
SECTION("Particle groups: muons") {
@@ -155,10 +154,24 @@ TEST_CASE("ParticleProperties", "[Particles]") {
CHECK(Hydrogen::nucleus_Z == 1);
CHECK(Tritium::nucleus_A == 3);
- // Nucleus is a generic object, it has no specific properties
- CHECK_THROWS(get_nucleus_Z(Code::Nucleus));
- CHECK_THROWS(get_nucleus_A(Code::Nucleus));
- CHECK_THROWS(get_mass(Code::Nucleus));
- CHECK_THROWS(get_charge(Code::Nucleus));
+ CHECK(is_nucleus(get_nucleus_code(1, 1)));
+ CHECK(is_nucleus(get_nucleus_code(100, 100)));
+ CHECK(get_nucleus_code(208, 82) == Code::Lead);
+ CHECK_FALSE(is_nucleus(Code::Electron));
+ CHECK(is_nucleus(Code::Lead));
+ CHECK(get_nucleus_Z(Code::Lead) == 82);
+ CHECK(get_nucleus_A(Code::Lead) == 208);
+
+ // impossible nucleus
+ CHECK_THROWS(get_nucleus_code(20, 40));
+
+ // getters
+ auto const testId = get_nucleus_code(40, 20);
+ CHECK(get_nucleus_A(testId) == 40);
+ CHECK(get_nucleus_Z(testId) == 20);
+ CHECK(is_nucleus(testId));
+ CHECK(get_nucleus_mass(testId) == 20 * Proton::mass + 20 * Neutron::mass);
+ CHECK(get_name(testId) == "nucleus");
+ CHECK(get_charge(testId) == 20 * constants::e);
}
}
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index aa2a006f29602efe6b66e841d304df29a357921c..3eab3eaafb7ddf1cfc169798dddd9cc7120d80ec 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -29,7 +29,7 @@
using namespace corsika;
using namespace corsika::epos;
-TEST_CASE("epos", "modules") {
+TEST_CASE("epos", "module,process") {
logging::set_level(logging::level::trace);
@@ -209,7 +209,7 @@ TEST_CASE("EposInterface", "modules") {
const HEPEnergyType P0 = 60_GeV;
auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Proton, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ Code::Proton, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
MomentumVector plab =
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr;
@@ -238,7 +238,7 @@ TEST_CASE("EposInterface", "modules") {
const HEPEnergyType P0 = 10_TeV;
auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Nucleus, 8, 4, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ get_nucleus_code(8, 4), P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
MomentumVector plab =
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr;
diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp
index 1ed62d52453cf5c235d359e4958c2f2dbaa3a177..e350c2fb2f0f5d24a8f69c6719658e5c1d6e6391 100644
--- a/tests/modules/testObservationPlane.cpp
+++ b/tests/modules/testObservationPlane.cpp
@@ -41,8 +41,7 @@ TEST_CASE("ObservationPlane", "interface") {
ObservationPlane has origin at 10,0,0 and a normal in x-direction
*/
- auto [stack, viewPtr] =
- setup::testing::setup_stack(Code::NuE, 0, 0, 1_GeV, nodePtr, cs);
+ auto [stack, viewPtr] = setup::testing::setup_stack(Code::NuE, 1_GeV, nodePtr, cs);
[[maybe_unused]] setup::StackView& view = *viewPtr;
auto particle = stack->getNextParticle();
diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp
index ea3a8bbe552d874ae5d26d9841caf7b556381e5c..453764066699e3c5f4683720787083e32cc668ce 100644
--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -24,7 +24,7 @@
using namespace corsika;
-TEST_CASE("ParticleCut", "processes") {
+TEST_CASE("ParticleCut", "process,continuous,secondary") {
logging::set_level(logging::level::info);
@@ -123,12 +123,12 @@ TEST_CASE("ParticleCut", "processes") {
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,
+ projectile.addSecondary(std::make_tuple(get_nucleus_code(A, Z), Eabove * A,
DirectionVector(rootCS, {1, 0, 0}), point0,
- 0_ns, A, Z));
- projectile.addSecondary(std::make_tuple(Code::Nucleus, Ebelow * A,
+ 0_ns));
+ projectile.addSecondary(std::make_tuple(get_nucleus_code(A, Z), Ebelow * A,
DirectionVector(rootCS, {1, 0, 0}), point0,
- 0_ns, A, Z));
+ 0_ns));
cut.doSecondaries(view);
@@ -160,12 +160,12 @@ TEST_CASE("ParticleCut", "processes") {
unsigned short A = 18;
unsigned short Z = 8;
- projectile.addSecondary(std::make_tuple(Code::Nucleus, 4_GeV * A,
+ projectile.addSecondary(std::make_tuple(get_nucleus_code(A, Z), 4_GeV * A,
DirectionVector(rootCS, {1, 0, 0}), point0,
- 0_ns, A, Z));
- projectile.addSecondary(std::make_tuple(Code::Nucleus, 6_GeV * A,
+ 0_ns));
+ projectile.addSecondary(std::make_tuple(get_nucleus_code(A, Z), 6_GeV * A,
DirectionVector(rootCS, {1, 0, 0}), point0,
- 0_ns, A, Z));
+ 0_ns));
cut.doSecondaries(view);
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index f6809d4a61b4e6015e625a5c302ab8e7fd8ad713..0e449f16b6a38eaeac30b95a3475c7b302532f76 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -104,7 +104,7 @@ TEST_CASE("Pythia8Interface", "modules") {
SECTION("pythia decay") {
HEPEnergyType const P0 = 10_GeV;
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::PiPlus, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ Code::PiPlus, P0, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
auto& stack = *stackPtr;
auto& view = *secViewPtr;
@@ -164,8 +164,7 @@ TEST_CASE("Pythia8Interface", "modules") {
// 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,
- *csPtr);
+ Code::Proton, 7_TeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
auto& view = *secViewPtr;
auto const particle = stackPtr->getNextParticle();
@@ -198,8 +197,7 @@ TEST_CASE("Pythia8Interface", "modules") {
// this is a projectile nucleus with very little energy
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Oxygen, 0, 0, 17_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::Oxygen, 17_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
auto& view = *secViewPtr;
auto particle = stackPtr->first();
@@ -215,8 +213,7 @@ TEST_CASE("Pythia8Interface", "modules") {
// this is a projectile neutron with very little energy
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Neutron, 0, 0, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::Neutron, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
auto& view = *secViewPtr;
auto particle = stackPtr->first();
@@ -240,7 +237,7 @@ TEST_CASE("Pythia8Interface", "modules") {
// resonable projectile, but tool low energy
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Proton, 0, 0, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr_Fe,
+ Code::Proton, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr_Fe,
*csPtr_Fe);
auto& view = *secViewPtr;
{ [[maybe_unused]] auto const& dummy_StackPtr = stackPtr; }
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 33ba73d46943e9a98dcda0ee17f0ed1ad9654594..ccad0cd7eb42fc2c0574b627fe849be4030d8b86 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -129,7 +129,7 @@ TEST_CASE("QgsjetII", "[processes]") {
#include <SetupTestEnvironment.hpp>
#include <SetupTestStack.hpp>
-TEST_CASE("QgsjetIIInterface", "[processes]") {
+TEST_CASE("QgsjetIIInterface", "interaction,processes") {
logging::set_level(logging::level::info);
@@ -142,8 +142,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
SECTION("InteractionInterface") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Proton, 0, 0, 110_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::Proton, 110_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
auto projectile = secViewPtr->getProjectile();
@@ -155,16 +154,16 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
CHECK(length / (1_g / square(1_cm)) == Approx(93.04).margin(0.1));
- /***********************************
- It as turned out already two times (#291 and #307) that the detailed output of
- QGSJetII event generation depends on the gfortran version used. This is not reliable
- and cannot be tested in a unit test here. One related problem was already found (#291)
- and is realted to undefined behaviour in the evaluation of functions in logical
- expressions. It is not clear if #307 is the same issue.
+ /* **********************************
+ As it turned out already two times (#291 and #307) that the detailed output of
+ QGSJetII event generation depends on the gfortran version used. This is not reliable
+ and cannot be tested in a unit test here. One related problem was already found
+ (#291) and is realted to undefined behaviour in the evaluation of functions in logical
+ expressions. It is not clear if #307 is the same issue.
CHECK(view.getSize() == 14);
CHECK(sumCharge(view) == 2);
- ************************************/
+ *********************************** */
auto const secMomSum = sumMomentum(view, projectileMomentum.getCoordinateSystem());
CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
Approx(0).margin(1e-2));
@@ -173,7 +172,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
SECTION("InteractionInterface Nuclei") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Nucleus, 60, 30, 20100_GeV,
+ get_nucleus_code(60, 30), 20100_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
@@ -182,9 +181,9 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
corsika::qgsjetII::Interaction model;
model.doInteraction(view); // this also should produce some fragments
- CHECK(view.getSize() == Approx(350).margin(100)); // this is not physics validation
+ CHECK(view.getSize() == Approx(300).margin(100)); // this is not physics validation
int countFragments = 0;
- for (auto const& sec : view) { countFragments += (sec.getPID() == Code::Nucleus); }
+ for (auto const& sec : view) { countFragments += (is_nucleus(sec.getPID())); }
CHECK(countFragments == Approx(4).margin(2)); // this is not physics validation
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
@@ -195,7 +194,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
SECTION("Heavy nuclei") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Nucleus, 1000, 1000, 1100_GeV,
+ get_nucleus_code(1000, 1000), 1100_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
@@ -215,7 +214,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
SECTION("Allowed Particles") {
{ // electron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Electron, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ Code::Electron, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
@@ -225,8 +224,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
}
{ // pi0 is internally converted into pi+/pi-
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Pi0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::Pi0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
@@ -235,8 +233,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
}
{ // rho0 is internally converted into pi-/pi+
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Rho0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::Rho0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
@@ -245,7 +242,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
}
{ // Lambda is internally converted into neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Lambda0, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ Code::Lambda0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
@@ -255,8 +252,8 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
}
{ // AntiLambda is internally converted into anti neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Lambda0Bar, 0, 0, 1000_GeV,
- (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ Code::Lambda0Bar, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
[[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
[[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index f544204c238709dc9c1464694f87a7556f63a42b..68e69fe57c2a8237efd9e6c5bd1f7a1853478d88 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -150,7 +150,7 @@ TEST_CASE("SibyllInterface", "modules") {
const HEPEnergyType P0 = 60_GeV;
auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Proton, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ Code::Proton, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
MomentumVector plab =
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr;
@@ -237,7 +237,7 @@ TEST_CASE("SibyllInterface", "modules") {
const HEPEnergyType P0 = 5_GeV;
auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Proton, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ Code::Proton, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
MomentumVector plab =
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr;
@@ -256,7 +256,7 @@ TEST_CASE("SibyllInterface", "modules") {
const HEPEnergyType P0 = 1000_EeV;
auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Proton, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ Code::Proton, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
{ [[maybe_unused]] auto const& dummy1 = stack; }
MomentumVector plab =
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
@@ -272,7 +272,7 @@ TEST_CASE("SibyllInterface", "modules") {
auto const& cs1 = *csPtr1;
const HEPEnergyType P0 = 150_GeV;
auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Electron, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr1, cs1);
+ Code::Electron, P0, (setup::Environment::BaseNodeType* const)nodePtr1, cs1);
{ [[maybe_unused]] auto const& dummy1 = stack; }
MomentumVector plab = MomentumVector(
cs1, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
@@ -285,7 +285,7 @@ TEST_CASE("SibyllInterface", "modules") {
SECTION("NuclearInteractionInterface") {
auto [stack, viewPtr] =
- setup::testing::setup_stack(Code::Nucleus, 8, 4, 900_GeV,
+ setup::testing::setup_stack(get_nucleus_code(8, 4), 900_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, cs);
setup::StackView& view = *viewPtr;
auto particle = stack->first();
@@ -306,9 +306,8 @@ TEST_CASE("SibyllInterface", "modules") {
SECTION("DecayInterface") {
- auto [stackPtr, viewPtr] =
- setup::testing::setup_stack(Code::Lambda0, 0, 0, 10_GeV,
- (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ auto [stackPtr, viewPtr] = setup::testing::setup_stack(
+ Code::Lambda0, 10_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
setup::StackView& view = *viewPtr;
auto& stack = *stackPtr;
auto particle = stack.first();
@@ -327,7 +326,7 @@ TEST_CASE("SibyllInterface", "modules") {
SECTION("DecayInterface - decay not handled") {
// sibyll does not know the higgs for example
auto [stackPtr, viewPtr] = setup::testing::setup_stack(
- Code::H0, 0, 0, 10_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ Code::H0, 10_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
setup::StackView& view = *viewPtr;
auto& stack = *stackPtr;
auto particle = stack.first();
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
index 9464e54a55508bf5948ff4650a834f5a04bd34de..ec2157e3e79cff3a1680692b482bf6f5eac159f3 100644
--- a/tests/modules/testStackInspector.cpp
+++ b/tests/modules/testStackInspector.cpp
@@ -38,9 +38,9 @@ TEST_CASE("StackInspector", "modules") {
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));
- stack.addParticle(std::make_tuple(Code::Nucleus,
+ stack.addParticle(std::make_tuple(get_nucleus_code(16, 8),
MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
- Point(rootCS, {0_m, 0_m, 10_km}), 0_ns, 16, 8));
+ Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
SECTION("interface") {
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index 243cbffa7df9bc47b2486b048be356381587d6db..d1d967d6f38e58a9b5cff769589939bf0042ed48 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -167,7 +167,7 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
targetPtr->addChild(std::move(target_2_behind));
targetPtr->addChild(std::move(target_2_partly_behind));
- auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs);
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, P0, targetPtr, cs);
{ [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
auto particle = stack->first();
// Note: momentum in X-direction
@@ -257,7 +257,7 @@ TEST_CASE("TrackingLeapFrogCurved") {
tracking_leapfrog_curved::Tracking tracking;
Point const center(cs, {0_m, 0_m, 0_m});
- auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, worldPtr, cs);
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, P0, worldPtr, cs);
{ [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
auto particle = stack->first();
// Note: momentum in X-direction
@@ -291,7 +291,7 @@ TEST_CASE("TrackingLeapFrogCurved") {
auto* targetPtr = target.get();
worldPtr->addChild(std::move(target));
- auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs);
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, P0, targetPtr, cs);
{ [[maybe_unused]] auto& viewPtr_dum = viewPtr; } // prevent warning
auto particle = stack->first();
// Note: momentum in X-direction
@@ -323,8 +323,7 @@ TEMPLATE_TEST_CASE("TrackingFail", "doesntwork", tracking_leapfrog_curved::Track
TestType tracking;
Point const center(cs, {0_m, 0_m, 0_m});
- auto [stack, viewPtr] =
- setup::testing::setup_stack(Code::Proton, 0, 0, P0, worldPtr, cs);
+ auto [stack, viewPtr] = setup::testing::setup_stack(Code::Proton, P0, worldPtr, cs);
{ [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
auto particle = stack->first();
NonExistingDummyObject const dummy;
@@ -376,7 +375,7 @@ TEMPLATE_TEST_CASE("TrackingPlane", "plane", tracking_leapfrog_curved::Tracking,
TestType tracking;
Point const center(cs, {0_m, 0_m, 0_m});
- auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, worldPtr, cs);
+ auto [stack, viewPtr] = setup::testing::setup_stack(PID, P0, worldPtr, cs);
{ [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
auto particle = stack->first();
// Note: momentum in X-direction
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index b5c207a59e7b48487897404cdfee93a376aa1049..89bde00d2e1772b943f933311ce3f258e438502b 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -82,7 +82,7 @@ TEST_CASE("UrQMD") {
Code::K0Bar, Code::K0Long};
for (auto code : validProjectileCodes) {
- auto [stack, view] = setup::testing::setup_stack(code, 0, 0, 100_GeV, nodePtr, cs);
+ auto [stack, view] = setup::testing::setup_stack(code, 100_GeV, nodePtr, cs);
CHECK(stack->getEntries() == 1);
CHECK(view->getEntries() == 0);
@@ -96,8 +96,7 @@ TEST_CASE("UrQMD") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
auto const& cs = *csPtr;
{ [[maybe_unused]] auto const& env_dummy = env; }
- auto [stack, view] =
- setup::testing::setup_stack(Code::Proton, 0, 0, 100_GeV, nodePtr, cs);
+ auto [stack, view] = setup::testing::setup_stack(Code::Proton, 100_GeV, nodePtr, cs);
[[maybe_unused]] setup::StackView& viewRef = *(view.get());
CHECK(urqmd.getInteractionLength(stack->getNextParticle()) / 1_g * square(1_cm) ==
Approx(105).margin(5));
@@ -107,8 +106,7 @@ TEST_CASE("UrQMD") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Omega);
auto const& cs = *csPtr;
{ [[maybe_unused]] auto const& env_dummy = env; }
- auto [stack, view] =
- setup::testing::setup_stack(Code::Neutron, 0, 0, 100_GeV, nodePtr, cs);
+ auto [stack, view] = setup::testing::setup_stack(Code::Neutron, 100_GeV, nodePtr, cs);
[[maybe_unused]] setup::StackView& viewRef = *(view.get());
CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle()));
}
@@ -120,7 +118,7 @@ TEST_CASE("UrQMD") {
unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Nucleus, A, Z, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ get_nucleus_code(A, Z), 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
[[maybe_unused]] setup::StackView& viewRef = *(secViewPtr.get());
CHECK(stackPtr->getEntries() == 1);
@@ -145,8 +143,7 @@ TEST_CASE("UrQMD") {
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::PiPlus, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
CHECK(stackPtr->getEntries() == 1);
CHECK(secViewPtr->getEntries() == 0);
@@ -172,8 +169,7 @@ TEST_CASE("UrQMD") {
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::PiPlus, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
[[maybe_unused]] auto particle = stackPtr->first();
CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); // Code::Proton not a valid target
}
@@ -184,8 +180,7 @@ TEST_CASE("UrQMD") {
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::PiPlus, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
CHECK(stackPtr->getEntries() == 1);
CHECK(secViewPtr->getEntries() == 0);
@@ -211,8 +206,7 @@ TEST_CASE("UrQMD") {
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::K0Long, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
+ Code::K0Long, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
CHECK(stackPtr->getEntries() == 1);
CHECK(secViewPtr->getEntries() == 0);
diff --git a/tests/stack/CMakeLists.txt b/tests/stack/CMakeLists.txt
index 58aa2b3025ea50c492d7cba79e6d3f9749277e59..1314fd435a88474ce1df7cc8bf54cc918f30e5e9 100644
--- a/tests/stack/CMakeLists.txt
+++ b/tests/stack/CMakeLists.txt
@@ -1,12 +1,10 @@
set (test_stack_sources
TestMain.cpp
testHistoryStack.cpp
-# testHistoryView.cpp
testGeometryNodeStackExtension.cpp
testWeightStackExtension.cpp
testDummyStack.cpp
testVectorStack.cpp
- testNuclearStackExtension.cpp
)
CORSIKA_ADD_TEST (testStack SOURCES ${test_stack_sources})
diff --git a/tests/stack/testNuclearStackExtension.cpp b/tests/stack/testNuclearStackExtension.cpp
deleted file mode 100644
index 98727de42d47a31426877c35224782b675077107..0000000000000000000000000000000000000000
--- a/tests/stack/testNuclearStackExtension.cpp
+++ /dev/null
@@ -1,289 +0,0 @@
-/*
- * (c) Copyright 2018 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.
- */
-
-#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
-#include <corsika/stack/NuclearStackExtension.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-using namespace corsika;
-
-#include <catch2/catch.hpp>
-
-#include <iostream>
-using namespace std;
-
-TEST_CASE("NuclearStackExtension", "stack") {
-
- logging::set_level(logging::level::info);
-
- CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
-
- SECTION("write non nucleus") {
- nuclear_stack::NuclearStackExtension<VectorStack,
- nuclear_stack::ExtendedParticleInterfaceType>
- 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);
- }
-
- SECTION("write nucleus") {
-
- nuclear_stack::NuclearStackExtension<VectorStack,
- nuclear_stack::ExtendedParticleInterfaceType>
- s;
- s.addParticle(std::make_tuple(
- 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, 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, 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() == 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, 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 + 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);
-
- // -> set to rest
- auto pout_mod = s.begin();
- pout_mod.setMomentum(MomentumVector{dummyCS, {0_GeV, 0_GeV, 0_GeV}});
- CHECK(pout_mod.getKineticEnergy() / 1_GeV == Approx(0));
- }
-
- SECTION("read invalid nucleus") {
- nuclear_stack::ParticleDataStack 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());
- }
-
- SECTION("stack fill and cleanup") {
-
- nuclear_stack::ParticleDataStack s;
- // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
- for (int i = 0; i < 99; ++i) {
- if ((i + 1) % 10 == 0) {
- 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, i, i / 2));
- } else {
- 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.getEntries() == 99);
- for (int i = 0; i < 99; ++i) s.getNextParticle().erase();
- CHECK(s.getEntries() == 0);
- }
-
- SECTION("stack operations") {
-
- nuclear_stack::ParticleDataStack s;
- // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
- // 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 - 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 - Electron::mass,
- DirectionVector(dummyCS, {1, 0, 0}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}),
- 100_s));
- }
- }
-
- // copy
- {
- s.copy(s.begin() + 9, s.begin() + 10); // nuclei to non-nuclei
- const auto& p9 = s.cbegin() + 9;
- const auto& p10 = s.cbegin() + 10;
-
- CHECK(p9.getPID() == Code::Nucleus);
- CHECK(p9.getEnergy() == 9 * 15_GeV);
- CHECK(p9.getTime() == 100_s);
- CHECK(p9.getNuclearA() == 9);
- CHECK(p9.getNuclearZ() == 9 / 2);
-
- CHECK(p10.getPID() == Code::Nucleus);
- CHECK(p10.getEnergy() == 9 * 15_GeV);
- CHECK(p10.getTime() == 100_s);
- CHECK(p10.getNuclearA() == 9);
- CHECK(p10.getNuclearZ() == 9 / 2);
- }
-
- // copy
- {
- s.copy(s.begin() + 93, s.begin() + 9); // non-nuclei to nuclei
- const auto& p93 = s.cbegin() + 93;
- const auto& p9 = s.cbegin() + 9;
-
- CHECK(p9.getPID() == Code::Electron);
- CHECK(p9.getEnergy() == 93 * 1.5_GeV);
- CHECK(p9.getTime() == 100_s);
-
- CHECK(p93.getPID() == Code::Electron);
- CHECK(p93.getEnergy() == 93 * 1.5_GeV);
- CHECK(p93.getTime() == 100_s);
- }
-
- // copy
- {
- s.copy(s.begin() + 89, s.begin() + 79); // nuclei to nuclei
- const auto& p89 = s.cbegin() + 89;
- const auto& p79 = s.cbegin() + 79;
-
- CHECK(p89.getPID() == Code::Nucleus);
- CHECK(p89.getEnergy() == 89 * 15_GeV);
- CHECK(p89.getTime() == 100_s);
-
- CHECK(p79.getPID() == Code::Nucleus);
- CHECK(p79.getEnergy() == 89 * 15_GeV);
- CHECK(p79.getTime() == 100_s);
- }
-
- // invalid copy
- { CHECK_THROWS(s.copy(s.begin(), s.end() + 1000)); }
-
- // swap
- {
- s.swap(s.begin() + 11, s.begin() + 10);
- const auto& p11 = s.cbegin() + 11; // now: nucleus
- const auto& p10 = s.cbegin() + 10; // now: electron
-
- CHECK(p11.getPID() == Code::Nucleus);
- CHECK(p11.getEnergy() == 9 * 15_GeV);
- CHECK(p11.getTime() == 100_s);
- CHECK(p11.getNuclearA() == 9);
- CHECK(p11.getNuclearZ() == 9 / 2);
-
- CHECK(p10.getPID() == Code::Electron);
- CHECK(p10.getEnergy() == 11 * 1.5_GeV);
- CHECK(p10.getTime() == 100_s);
- }
-
- // swap two nuclei
- {
- s.swap(s.begin() + 29, s.begin() + 59);
- const auto& p29 = s.cbegin() + 29;
- const auto& p59 = s.cbegin() + 59;
-
- CHECK(p29.getPID() == Code::Nucleus);
- CHECK(p29.getEnergy() == 59 * 15_GeV);
- CHECK(p29.getTime() == 100_s);
- CHECK(p29.getNuclearA() == 59);
- CHECK(p29.getNuclearZ() == 59 / 2);
-
- CHECK(p59.getPID() == Code::Nucleus);
- CHECK(p59.getEnergy() == 29 * 15_GeV);
- CHECK(p59.getTime() == 100_s);
- CHECK(p59.getNuclearA() == 29);
- CHECK(p59.getNuclearZ() == 29 / 2);
- }
-
- for (int i = 0; i < 99; ++i) s.last().erase();
- CHECK(s.getEntries() == 0);
- }
-
- SECTION("not allowed") {
-
- nuclear_stack::NuclearStackExtension<VectorStack,
- nuclear_stack::ExtendedParticleInterfaceType>
- s;
-
- // not valid, no A,Z (implicit A=Z=0):
- CHECK_THROWS(s.addParticle(
- std::make_tuple(Code::Nucleus, 100_GeV, DirectionVector(dummyCS, {1, 0, 0}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s)));
- CHECK_THROWS(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)));
-
- // not valid, non-Code::Nuclei with A and Z
- CHECK_THROWS(s.addParticle(std::make_tuple(
- Code::Oxygen, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 20, 10)));
-
- // not valid, non-Code::Nuclei with A and Z
- CHECK_THROWS(s.addParticle(std::make_tuple(
- Code::Oxygen, 100_GeV, DirectionVector(dummyCS, {0, 0, 1}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 20, 10)));
-
- // valid, non-Nucleus with implicit A=Z=0
- s.addParticle(
- std::make_tuple(Code::Proton, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
-
- // valid, for further use
- 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, 5));
-
- // not valid, Oxygen, but A and Z
- CHECK_THROWS(particle.addSecondary(std::make_tuple(
- Code::Oxygen, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 20, 10)));
-
- // not valid, Nucleus, no A, Z
- CHECK_THROWS(particle.addSecondary(
- std::make_tuple(Code::Nucleus, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s)));
-
- // not valid, Nucleus, no A, Z, and Energy
- CHECK_THROWS(particle.addSecondary(
- std::make_tuple(Code::Nucleus, 100_GeV, DirectionVector(dummyCS, {1, 0, 0}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s)));
-
- // add a another nucleus, so there are two now
- 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()));
- }
-}