diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index f3e5938fe0573d96ddb61df91b996373f575d6a2..f11427148c96e927176168f6d9d8c1ee3dc9fc15 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -47,7 +47,7 @@ namespace corsika {
auto pNext = stack_.getNextParticle();
CORSIKA_LOG_TRACE(
- "============== next particle : count={}, pid={}, "
+ "============== next particle : count={}, pid={}"
", stack entries={}"
", stack deleted={}",
count_, pNext.getPID(), stack_.getEntries(), stack_.getErased());
@@ -134,11 +134,9 @@ namespace corsika {
ExponentialDistribution expDist(total_lambda);
GrammageType const next_interact = expDist(rng_);
- CORSIKA_LOG_DEBUG(
- "total_lambda={} g/cm2, "
- ", next_interact={} g/cm2",
- double(total_lambda / 1_g * 1_cm * 1_cm),
- double(next_interact / 1_g * 1_cm * 1_cm));
+ CORSIKA_LOG_DEBUG("total_lambda={} g/cm2, next_interact={} g/cm2",
+ double(total_lambda / 1_g * 1_cm * 1_cm),
+ double(next_interact / 1_g * 1_cm * 1_cm));
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime = sequence_.getInverseLifetime(particle);
@@ -147,10 +145,8 @@ namespace corsika {
ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
TimeType const next_decay = expDistDecay(rng_);
- CORSIKA_LOG_DEBUG(
- "total_lifetime={} s"
- ", next_decay={} s",
- (1 / total_inv_lifetime) / 1_s, next_decay / 1_s);
+ CORSIKA_LOG_DEBUG("total_lifetime={} ns, next_decay={} ns",
+ (1 / total_inv_lifetime) / 1_ns, next_decay / 1_ns);
// convert next_decay from time to length [m]
LengthType const distance_decay = next_decay * particle.getMomentum().getNorm() /
@@ -198,10 +194,6 @@ namespace corsika {
// move particle along the trajectory to new position
// also update momentum/direction/time
step.setLength(min_distance);
- particle.setPosition(step.getPosition(1));
- // assumption: tracking does not change absolute momentum (continuous physics can and
- // will):
- particle.setMomentum(step.getDirection(1) * particle.getMomentum().getNorm());
// apply all continuous processes on particle + track
if (sequence_.doContinuous(particle, step, limitingId) ==
@@ -215,9 +207,12 @@ namespace corsika {
} else {
particle.erase();
}
- return;
+ return; // particle is gone -> return
}
particle.setTime(particle.getTime() + step.getDuration());
+ particle.setPosition(step.getPosition(1));
+ particle.setMomentum(step.getDirection(1) * particle.getMomentum().getNorm());
+
if (isContinuous) {
return; // there is nothing further, step is finished
}
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index a716b156d5e5b052b62dcf0bc771d41a1ea2fb8e..fe55d587412537f86dff4ca4f0ca0f9dd5ef99ab 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -80,26 +80,20 @@ namespace corsika::pythia8 {
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
}
- inline void Interaction::isValid(Code const projectileId, Code const targetId,
+ inline bool Interaction::isValid(Code const projectileId, Code const targetId,
HEPEnergyType const sqrtS) const {
- if ((10_GeV > sqrtS) || (sqrtS > 1_PeV)) {
- throw std::runtime_error("energy out of bounds for PYTHIA");
- }
+ if ((10_GeV > sqrtS) || (sqrtS > 1_PeV)) { return false; }
if (targetId != Code::Hydrogen && targetId != Code::Neutron &&
targetId != Code::Proton) {
- throw std::runtime_error("wrong target for PYTHIA");
+ return false;
}
- if (is_nucleus(projectileId)) {
- // nuclei handled by different process, this should not happen
- throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
- }
+ if (is_nucleus(projectileId)) { return false; }
- if (!canInteract(projectileId)) {
- throw std::runtime_error("Projectile not supported by PYTHIA!");
- }
+ if (!canInteract(projectileId)) { return false; }
+ return true;
}
inline void Interaction::configureLabFrameCollision(Code const projectileId,
@@ -150,7 +144,9 @@ namespace corsika::pythia8 {
HEPEnergyType const CoMenergy = (projectileP4 + targetP4).getNorm();
- isValid(projectileId, targetId, CoMenergy); // throws
+ if (!isValid(projectileId, targetId, CoMenergy)) {
+ return {CrossSectionType::zero(), CrossSectionType::zero()};
+ }
// input particle PDG
auto const pdgCodeBeam = static_cast<int>(get_PDG(projectileId));
@@ -195,11 +191,13 @@ namespace corsika::pythia8 {
static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId));
- isValid(projectileId, targetId, sqrtS); // throws
+ if (!isValid(projectileId, targetId, sqrtS)) {
+ throw std::runtime_error("invalid target,projectile,energy combination.");
+ }
// position and time of interaction
- Point pOrig = projectile.getPosition();
- TimeType tOrig = projectile.getTime();
+ Point const& pOrig = projectile.getPosition();
+ TimeType const tOrig = projectile.getTime();
CORSIKA_LOG_DEBUG("Interaction: ebeam lab: {} GeV", eProjectileLab / 1_GeV);
diff --git a/corsika/detail/modules/qgsjetII/InteractionModel.inl b/corsika/detail/modules/qgsjetII/InteractionModel.inl
index 0b9bb370ac0cfb32721b9c9718cf7fe1c4c1d700..b65d5b3f72a64e576711aeee639514bfdc4d355f 100644
--- a/corsika/detail/modules/qgsjetII/InteractionModel.inl
+++ b/corsika/detail/modules/qgsjetII/InteractionModel.inl
@@ -41,44 +41,32 @@ namespace corsika::qgsjetII {
CORSIKA_LOG_DEBUG("QgsjetII::InteractionModel n= {}", count_);
}
- inline void InteractionModel::isValid(Code const projectileId,
+ inline bool InteractionModel::isValid(Code const projectileId,
Code const targetId) const {
if (is_nucleus(targetId)) {
size_t iTarget = get_nucleus_A(targetId);
- if (iTarget > int(maxMassNumber_) || iTarget <= 0) {
- std::ostringstream txt;
- txt << "QgsjetII target outside range. Atarget=" << iTarget;
- throw std::runtime_error(txt.str().c_str());
- }
+ if (iTarget > int(maxMassNumber_) || iTarget <= 0) { return false; }
} else if (targetId != Proton::code) {
- std::ostringstream txt;
- txt << "QgsjetII not valid target=" << targetId;
- throw std::runtime_error(txt.str().c_str());
+ return false;
}
if (is_nucleus(projectileId)) {
size_t iProjectile = get_nucleus_A(projectileId);
- if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) {
- std::ostringstream txt;
- txt << "QgsjetII projectile outside range. Aprojectile=" << iProjectile;
- throw std::runtime_error(txt.str().c_str());
- }
+ if (iProjectile > int(maxMassNumber_) || iProjectile <= 0) { return false; }
} else if (!is_hadron(projectileId)) {
- std::ostringstream txt;
- txt << "QgsjetII projectile can only be hadrons. Is=" << projectileId;
- throw std::runtime_error(txt.str().c_str());
+ return false;
}
+ return true;
}
inline CrossSectionType InteractionModel::getCrossSection(
Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
- double sigProd = std::numeric_limits<double>::infinity();
-
- if (!corsika::qgsjetII::canInteract(projectileId)) { return sigProd * 1_mb; }
-
- isValid(projectileId, targetId); // throws
+ if (!corsika::qgsjetII::canInteract(projectileId)) {
+ return CrossSectionType::zero();
+ }
+ if (!isValid(projectileId, targetId)) { return CrossSectionType::zero(); }
// define projectile, in lab frame
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
@@ -97,9 +85,8 @@ namespace corsika::qgsjetII {
"QgsjetII::getCrossSection Elab= {} GeV iBeam= {}"
" iProjectile= {} iTarget= {}",
Elab / 1_GeV, iBeam, iProjectile, iTarget);
- sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
+ double sigProd = qgsect_(Elab / 1_GeV, iBeam, iProjectile, iTarget);
CORSIKA_LOG_DEBUG("QgsjetII::getCrossSection sigProd= {} mb", sigProd);
-
return sigProd * 1_mb;
}
@@ -114,9 +101,10 @@ namespace corsika::qgsjetII {
"doInteraction: {} interaction possible? {}",
projectileId, corsika::qgsjetII::canInteract(projectileId));
- if (!corsika::qgsjetII::canInteract(projectileId)) return;
-
- isValid(projectileId, targetId); // throws
+ if (!corsika::qgsjetII::canInteract(projectileId) ||
+ !isValid(projectileId, targetId)) {
+ throw std::runtime_error("invalid target/projectile/energy combination.");
+ }
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
diff --git a/corsika/detail/modules/sibyll/Decay.inl b/corsika/detail/modules/sibyll/Decay.inl
index 6b8ae281d4a4be07408785e1a846c376c8b4bc0f..ab0114182b41f5db321f19ec6d91458cc9b5000d 100644
--- a/corsika/detail/modules/sibyll/Decay.inl
+++ b/corsika/detail/modules/sibyll/Decay.inl
@@ -174,7 +174,7 @@ namespace corsika::sibyll {
count_++;
// remember position
- Point const decayPoint = projectile.getPosition();
+ Point const& decayPoint = projectile.getPosition();
TimeType const t0 = projectile.getTime();
// switch on decay for this particle
setUnstable(pCode);
diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
index 38acb4b3d457bef02f81270761f7fc82716d5821..9ecb434bc62cb5a0bbccc6a13ce889032b4bdcda 100644
--- a/corsika/detail/modules/sibyll/InteractionModel.inl
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -36,26 +36,24 @@ namespace corsika::sibyll {
CORSIKA_LOG_DEBUG("Sibyll::Model n={}, Nnuc={}", count_, nucCount_);
}
- inline void constexpr InteractionModel::isValid(Code const projectileId,
+ inline bool constexpr InteractionModel::isValid(Code const projectileId,
Code const targetId,
HEPEnergyType const sqrtSnn) const {
- if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) {
- // i.e. nuclei handled by different process, this should not happen
- throw std::runtime_error("CoM energy out of bounds for SIBYLL");
- }
+ if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
if (is_nucleus(targetId)) {
- unsigned int const targA = get_nucleus_A(targetId);
+ size_t const targA = get_nucleus_A(targetId);
if (targA != 1 && (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) {
- throw std::runtime_error("Target outside of allowed range for SIBYLL");
+ return false;
}
} else if (targetId != Code::Proton && targetId != Code::Neutron &&
targetId != Code::Hydrogen) {
- throw std::runtime_error("Target cannot be handled by SIBYLL");
+ return false;
}
if (is_nucleus(projectileId) || !corsika::sibyll::canInteract(projectileId)) {
- throw std::runtime_error("Projectile cannot be handled by SIBYLL");
+ return false;
}
+ return true;
}
inline std::tuple<CrossSectionType, CrossSectionType>
@@ -68,6 +66,10 @@ namespace corsika::sibyll {
// sqrtS per target nucleon
HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
+ if (!isValid(projectileId, targetId, sqrtSnn)) {
+ return {CrossSectionType::zero(), CrossSectionType::zero()};
+ }
+
double dummy, dum1, dum3, dum4, dumdif[3]; // dummies needed for fortran call
int const iBeam = corsika::sibyll::getSibyllXSCode(
projectileId); // 0 (can not interact, 1: proton-like, 2: pion-like,
@@ -109,7 +111,9 @@ namespace corsika::sibyll {
HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
COMBoost const boost(projectileP4, targetP4 / targetSibCode);
- isValid(projectileId, targetId, sqrtSnn); // throws
+ if (!isValid(projectileId, targetId, sqrtSnn)) {
+ throw std::runtime_error("Invalid target/projectile/energy combination");
+ }
CORSIKA_LOG_DEBUG("pId={} tId={} sqrtSnn={}GeV", projectileId, targetId, sqrtSnn);
@@ -135,10 +139,9 @@ namespace corsika::sibyll {
// add particles from sibyll to stack
// position and time of interaction, not used in Sibyll
- Point const pOrig = Point(csPrime, {0_m, 0_m, 0_m});
- TimeType const tOrig = 0_s; // no time in sibyll
- CORSIKA_LOG_DEBUG("position of interaction: {}, time {} ", pOrig.getCoordinates(),
- tOrig);
+ auto const& projectile = secondaries.parent();
+ Point const& pOrig = projectile.getPosition();
+ TimeType const tOrig = projectile.getTime(); // no time in sibyll
// link to sibyll stack
SibStack ss;
diff --git a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
index e9ad447f33e0d25bd4d14fa735ed630d083662d9..bda2996b5a1c2eacd9d2f9a5e81c2f85e3e91987 100644
--- a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
@@ -40,20 +40,17 @@ namespace corsika::sibyll {
}
template <typename TEnvironment, typename TNucleonModel>
- inline void constexpr NuclearInteractionModel<TEnvironment, TNucleonModel>::isValid(
+ inline bool constexpr NuclearInteractionModel<TEnvironment, TNucleonModel>::isValid(
Code const projectileId, Code const targetId, HEPEnergyType const sqrtSnn) const {
// also depends on underlying model, for Proton/Neutron projectile
- hadronicInteraction_.isValid(Code::Proton, targetId, sqrtSnn); // throws
+ if (!hadronicInteraction_.isValid(Code::Proton, targetId, sqrtSnn)) { return false; }
// projectile limits:
- if (!is_nucleus(projectileId)) {
- throw std::runtime_error("can only handle nuclear projectile");
- }
+ if (!is_nucleus(projectileId)) { return false; }
unsigned int projectileA = get_nucleus_A(projectileId);
- if (projectileA > getMaxNucleusAProjectile() || projectileA < 2) {
- throw std::runtime_error("projectile mass A out of bounds");
- }
+ if (projectileA > getMaxNucleusAProjectile() || projectileA < 2) { return false; }
+ return true;
} // namespace corsika::sibyll
template <typename TEnvironment, typename TNucleonModel>
@@ -114,28 +111,35 @@ namespace corsika::sibyll {
for (Code const ptarg : allElementsInUniverse) {
if (ptarg == Code::Argon) continue; // NEED TO IGNORE Argon ....
++k;
- CORSIKA_LOG_DEBUG("init target component: {}", ptarg);
+ CORSIKA_LOG_DEBUG("init target component: {} A={}", ptarg, get_nucleus_A(ptarg));
int const ib = get_nucleus_A(ptarg);
- hadronicInteraction_.isValid(Code::Proton, ptarg, 100_GeV); // throws
+ if (!hadronicInteraction_.isValid(Code::Proton, ptarg, 100_GeV)) {
+ throw std::runtime_error("Invalid target type.");
+ }
targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
// loop over energies, fNEnBins log. energy bins
- for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
+ for (size_t i = 0; i < getNEnergyBins(); ++i) {
// hard coded energy grid, has to be aligned to definition in signuc2!!, no
// comment..
HEPEnergyType const Ecm = pow(10., 1. + 1. * i) * 1_GeV;
+ // head-on pp collision:
+ HEPEnergyType const EcmHalve = Ecm / 2;
+ HEPMomentumType const pcm =
+ sqrt(EcmHalve * EcmHalve - Proton::mass * Proton::mass);
CoordinateSystemPtr cs = get_root_CoordinateSystem();
- HEPMomentumType const pcm = sqrt(Ecm * Ecm - Proton::mass * Proton::mass);
- FourMomentum projectileP4(Ecm,
- {cs, pcm, 0_eV, 0_eV}); // this is ONLY needd for sqrtS
- FourMomentum targetP4(0_eV, {cs, 0_eV, 0_eV, 0_eV});
+ FourMomentum projectileP4(EcmHalve, {cs, pcm, 0_eV, 0_eV});
+ FourMomentum targetP4(EcmHalve, {cs, -pcm, 0_eV, 0_eV});
// get p-p cross sections
- auto const protonId = Code::Proton;
+ if (!hadronicInteraction_.isValid(Code::Proton, Code::Proton, Ecm)) {
+ throw std::runtime_error("invalid projectile,target,ecm combination");
+ }
auto const [siginel, sigela] = hadronicInteraction_.getCrossSectionInelEla(
- protonId, protonId, projectileP4, targetP4);
- const double dsig = siginel / 1_mb;
- const double dsigela = sigela / 1_mb;
+ Code::Proton, Code::Proton, projectileP4, targetP4);
+ double const dsig = siginel / 1_mb;
+ double const dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
- for (unsigned int j = 1; j < gMaxNucleusAProjectile_; ++j) {
+ CORSIKA_LOG_TRACE("Ecm={} siginel={} sigela={}", Ecm / 1_GeV, dsig, dsigela);
+ for (size_t j = 1; j < gMaxNucleusAProjectile_; ++j) {
const int jj = j + 1;
double sig_out, dsig_out, sigqe_out, dsigqe_out;
sigma_mc_(jj, ib, dsig, dsigela, gNSample_, sig_out, dsig_out, sigqe_out,
@@ -143,6 +147,7 @@ namespace corsika::sibyll {
// write to table
cnucsignuc_.sigma[j][k][i] = sig_out;
cnucsignuc_.sigqe[j][k][i] = sigqe_out;
+ CORSIKA_LOG_TRACE("nuc A={} sig={} qe={}", j, sig_out, sigqe_out);
}
}
}
@@ -177,7 +182,7 @@ namespace corsika::sibyll {
FourMomentum const& targetP4) const {
HEPEnergyType const sqrtSnn = (projectileP4 + targetP4).getNorm();
- isValid(projectileId, targetId, sqrtSnn); // throws
+ if (!isValid(projectileId, targetId, sqrtSnn)) { return CrossSectionType::zero(); }
HEPEnergyType const LabEnergyPerNuc =
static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
auto const sigProd =
@@ -199,7 +204,9 @@ namespace corsika::sibyll {
// this is center-of-mass for projectile_nucleon - target
FourMomentum const nucleonP4 = projectileP4 / projectileA;
HEPEnergyType const sqrtSnucleon = (nucleonP4 + targetP4).getNorm();
- isValid(projectileId, targetId, sqrtSnucleon); // throws
+ if (!isValid(projectileId, targetId, sqrtSnucleon)) {
+ throw std::runtime_error("Invalid projectile/target/energy combination.");
+ }
// projectile is always nucleus!
// Elab corresponding to sqrtSnucleon -> fixed target projectile
COMBoost const boost(nucleonP4, targetP4);
@@ -287,8 +294,8 @@ namespace corsika::sibyll {
Point const& pOrig = projectile.getPosition();
TimeType const delay = projectile.getTime();
- CORSIKA_LOG_DEBUG("Interaction: position of interaction: {} {}",
- pOrig.getCoordinates(), delay / 1_s);
+ CORSIKA_LOG_DEBUG("Interaction: position of interaction: {}, {} ns",
+ pOrig.getCoordinates(), delay / 1_ns);
CORSIKA_LOG_DEBUG("number of fragments: {}", nFragments);
CORSIKA_LOG_DEBUG("adding nuclear fragments to particle stack..");
// put nuclear fragments on corsika stack
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index b8393b0e91da8c63cae6ddf5015df605abee4762..1bf0f0d59a42428df231717f2df869f72c30bd45 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -37,14 +37,13 @@ namespace corsika::urqmd {
::urqmd::iniurqmdc8_();
}
- inline void UrQMD::isValid(Code const projectileId, Code const targetId) const {
+ inline bool UrQMD::isValid(Code const projectileId, Code const targetId) const {
if (!is_hadron(projectileId) || !corsika::urqmd::canInteract(projectileId)) {
- throw std::runtime_error("UrQMD projectile is not a compatible hadron.");
- }
- if (!is_nucleus(targetId)) {
- throw std::runtime_error("UrQMD target is not a nucleus .");
+ return false;
}
+ if (!is_nucleus(targetId)) { return false; }
+ return true;
}
inline CrossSectionType UrQMD::getTabulatedCrossSection(
@@ -141,7 +140,7 @@ namespace corsika::urqmd {
FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
- if (is_nucleus(projectileId)) {
+ if (!isValid(projectileId, targetId)) {
/*
* unfortunately unavoidable at the moment until we have tools to get the actual
* inealstic cross-section from UrQMD
@@ -149,8 +148,6 @@ namespace corsika::urqmd {
return CrossSectionType::zero();
}
- isValid(projectileId, targetId); // throws
-
// define projectile, in lab frame
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
@@ -200,7 +197,9 @@ namespace corsika::urqmd {
static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId));
- isValid(projectileId, targetId); // throws
+ if (!isValid(projectileId, targetId)) {
+ throw std::runtime_error("invalid target,projectile,energy combination");
+ }
size_t const targetA = get_nucleus_A(targetId);
size_t const targetZ = get_nucleus_Z(targetId);
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 74dad4f2ff1f71049efde8198ccbaaa94ff84910..e6d3600b8cff72e305a1b57f64a526a39772341f 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -273,7 +273,7 @@ namespace corsika {
FourMomentum const& targetP4) const;
template <typename TParticle>
- TimeType getLifetime(TParticle&& particle) {
+ TimeType getLifetime(TParticle& particle) {
return 1. / getInverseLifetime(particle);
}
diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp
index 02a5313b70c752a6680cf462ab4aee1e0c2be187..25ebf826b3f886065d5759fffa4eaa94f3e659ac 100644
--- a/corsika/modules/pythia8/Interaction.hpp
+++ b/corsika/modules/pythia8/Interaction.hpp
@@ -35,7 +35,7 @@ namespace corsika::pythia8 {
bool canInteract(Code const) const;
void configureLabFrameCollision(Code const, Code const, HEPEnergyType const);
- void isValid(Code const projectileId, Code const targetId,
+ bool isValid(Code const projectileId, Code const targetId,
HEPEnergyType const sqrtS) const;
/**
* Returns inelastic AND elastic cross sections.
diff --git a/corsika/modules/qgsjetII/InteractionModel.hpp b/corsika/modules/qgsjetII/InteractionModel.hpp
index 4ec1bf6a279a51e8e0c2f2fbd0aa4857d14be3eb..45e9f8ef083b5af9d8ee2c270c9629db2c010260 100644
--- a/corsika/modules/qgsjetII/InteractionModel.hpp
+++ b/corsika/modules/qgsjetII/InteractionModel.hpp
@@ -33,7 +33,7 @@ namespace corsika::qgsjetII {
* @param beamId
* @param targetId
*/
- void isValid(Code const beamId, Code const targetId) const;
+ bool isValid(Code const beamId, Code const targetId) const;
/**
* Return the QGSJETII inelastic/production cross section.
diff --git a/corsika/modules/sibyll/InteractionModel.hpp b/corsika/modules/sibyll/InteractionModel.hpp
index 0f65dcab7d8b38581e141807d5eff1a234f8818d..177e77857b305d385da162a3c0b705a8c5c276f9 100644
--- a/corsika/modules/sibyll/InteractionModel.hpp
+++ b/corsika/modules/sibyll/InteractionModel.hpp
@@ -45,7 +45,7 @@ namespace corsika::sibyll {
* sibyll only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
* neutrons (p,n == nucleon).
*/
- void constexpr isValid(Code const projectileId, Code const targetId,
+ bool constexpr isValid(Code const projectileId, Code const targetId,
HEPEnergyType const sqrtSnn) const;
/**
@@ -96,7 +96,7 @@ namespace corsika::sibyll {
*/
template <typename TSecondaries>
- void doInteraction(TSecondaries&, Code const projectile, Code const target,
+ void doInteraction(TSecondaries& view, Code const projectile, Code const target,
FourMomentum const& projectileP4, FourMomentum const& targetP4);
private:
diff --git a/corsika/modules/sibyll/NuclearInteractionModel.hpp b/corsika/modules/sibyll/NuclearInteractionModel.hpp
index 753cabf98f471832c8e2b315aecf62a7cf9bf0d3..2d7d81c2ce0717e33567715b554b93fb93ce15c9 100644
--- a/corsika/modules/sibyll/NuclearInteractionModel.hpp
+++ b/corsika/modules/sibyll/NuclearInteractionModel.hpp
@@ -29,7 +29,7 @@ namespace corsika::sibyll {
NuclearInteractionModel(TNucleonModel&, TEnvironment const&);
~NuclearInteractionModel();
- void constexpr isValid(Code const projectileId, Code const targetId,
+ bool constexpr isValid(Code const projectileId, Code const targetId,
HEPEnergyType const sqrtSnn) const;
void initializeNuclearCrossSections();
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
index 86d238e0fa736027803c8712d4a68df47b6da140..1ae0ff80af53f99ce31bbb383a459510e362a3e1 100644
--- a/corsika/modules/urqmd/UrQMD.hpp
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -35,7 +35,7 @@ namespace corsika::urqmd {
UrQMD(boost::filesystem::path const path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat"),
int const retryFlag = 0);
- void isValid(Code const projectileId, Code const targetId) const;
+ bool isValid(Code const projectileId, Code const targetId) const;
CrossSectionType getTabulatedCrossSection(Code const, Code const,
HEPEnergyType const) const;
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
index 06a3756f4b0adbdc8161ad122987b737476f956a..cf01b33ca499cf8d1df3d61bcc142f0aee921cdf 100644
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -146,11 +146,8 @@ int main() {
BetheBlochPDG eLoss{showerAxis};
// assemble all processes into an ordered process list
- auto sequence = make_sequence(
- stackInspect,
- make_select([](auto const& particle) { return is_nucleus(particle.getPID()); },
- sibyllNuc, sibyll),
- decay, eLoss, cut, trackWriter);
+ auto sequence = make_sequence(stackInspect, make_sequence(sibyllNuc, sibyll), decay,
+ eLoss, cut, trackWriter);
// define air shower object, run simulation
Cascade EAS(env, tracking, sequence, output, stack);
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index dbf99a5cba85272435cfd0827862981698a23e39..7df6dcbde54c0d40190da9a377c166d9472c588b 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -286,8 +286,7 @@ int main(int argc, char** argv) {
InteractionCounter sibyllCounted(sibyll);
corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
InteractionCounter sibyllNucCounted(sibyllNuc);
- auto heModelCounted = make_select([](auto const& p) { return is_nucleus(p.getPID()); },
- sibyllNucCounted, sibyllCounted);
+ auto heModelCounted = make_sequence(sibyllNucCounted, sibyllCounted);
corsika::pythia8::Decay decayPythia;
diff --git a/examples/mars.cpp b/examples/mars.cpp
index a284f1ecd85580ace4ec82867bf5ad8573ea5b05..afa8c71b70b8e2526785312acc38e1c7b684484c 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -313,8 +313,7 @@ int main(int argc, char** argv) {
InteractionCounter sibyllCounted(sibyll);
corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
InteractionCounter sibyllNucCounted(sibyllNuc);
- auto heModelCounted = make_select([](auto const& p) { return is_nucleus(p.getPID()); },
- sibyllNucCounted, sibyllCounted);
+ auto heModelCounted = make_sequence(sibyllNucCounted, sibyllCounted);
corsika::pythia8::Decay decayPythia;
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index d6c1af1dd4dada9e68f4adb4d2cec4f92facb42b..6bb97291aa2b2d7d3106fcf68a6b887aaf99c7d9 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -238,10 +238,8 @@ int main(int argc, char** argv) {
: cutE_(cutE) {}
bool operator()(const Particle& p) const { return (p.getEnergy() < cutE_); }
};
- auto hadronSequence =
- make_select(EnergySwitch(55_GeV), urqmdCounted,
- make_select([](auto const& p) { return is_nucleus(p.getPID()); },
- sibyllNucCounted, sibyllCounted));
+ auto hadronSequence = make_select(EnergySwitch(55_GeV), urqmdCounted,
+ make_sequence(sibyllNucCounted, sibyllCounted));
auto decaySequence = make_sequence(decayPythia, decaySibyll);
// directory for outputs
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index 06121eeec3e6c317871320e05a6cfb5b0dcf28c6..d54324be5dc00409defa1fc183c2f545f5263350 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -231,11 +231,13 @@ TEST_CASE("Pythia8Interface", "modules") {
corsika::pythia8::Interaction collision;
- CHECK_THROWS(collision.getCrossSection(
- Code::Proton, Code::Iron,
- {sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
- {rootCS, {0_eV, 0_eV, 100_GeV}}},
- {Iron::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
+ CHECK(collision.getCrossSection(
+ Code::Proton, Code::Iron,
+ {sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
+ {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {Iron::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}) /
+ 1_mb ==
+ Approx(0));
CHECK_THROWS(collision.doInteraction(
view, Code::Proton, Code::Iron,
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 2a8cd57a8d93f938f75fd0447eb6abae465663b8..40f544798157268b20e02bef6eff04b3f28ec323 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -96,8 +96,8 @@ TEST_CASE("QgsjetII", "[processes]") {
corsika::qgsjetII::InteractionModel model;
- CHECK_THROWS(model.isValid(Code::Electron, Code::Proton));
- CHECK_THROWS(model.isValid(Code::Proton, Code::Electron));
+ CHECK_FALSE(model.isValid(Code::Electron, Code::Proton));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Electron));
}
}
@@ -193,9 +193,12 @@ TEST_CASE("QgsjetIIInterface", "interaction,processes") {
FourMomentum const aP4(100_GeV, {cs, 99_GeV, 0_GeV, 0_GeV});
FourMomentum const bP4(1_TeV, {cs, 0.9_TeV, 0_GeV, 0_GeV});
- CHECK_THROWS(model.getCrossSection(get_nucleus_code(10, 5),
- get_nucleus_code(1000, 500), aP4, bP4));
- CHECK_THROWS(model.getCrossSection(Code::Nucleus, Code::Nucleus, aP4, bP4));
+ CHECK(model.getCrossSection(get_nucleus_code(10, 5), get_nucleus_code(1000, 500), aP4,
+ bP4) /
+ 1_mb ==
+ Approx(0));
+ CHECK(model.getCrossSection(Code::Nucleus, Code::Nucleus, aP4, bP4) / 1_mb ==
+ Approx(0));
CHECK_THROWS(
model.doInteraction(view, get_nucleus_code(1000, 500), Code::Oxygen, aP4, bP4));
}
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 957ddac0fe3cde92357934ac0d9f60961b6f0b87..2a19a030b7272d88c2912a7ce163bb0a81e64645 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -137,22 +137,22 @@ TEST_CASE("SibyllInterface", "modules") {
corsika::sibyll::InteractionModel model;
// sibyll only accepts protons or nuclei with 4<=A<=18 as targets
- CHECK_THROWS(model.isValid(Code::Proton, Code::Electron, 100_GeV));
- CHECK_NOTHROW(model.isValid(Code::Proton, Code::Hydrogen, 100_GeV));
- CHECK_THROWS(model.isValid(Code::Proton, Code::Deuterium, 100_GeV));
- CHECK_NOTHROW(model.isValid(Code::Proton, Code::Helium, 100_GeV));
- CHECK_THROWS(model.isValid(Code::Proton, Code::Helium3, 100_GeV));
- CHECK_THROWS(model.isValid(Code::Proton, Code::Iron, 100_GeV));
- CHECK_NOTHROW(model.isValid(Code::Proton, Code::Oxygen, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 100_GeV));
+ CHECK(model.isValid(Code::Proton, Code::Hydrogen, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Deuterium, 100_GeV));
+ CHECK(model.isValid(Code::Proton, Code::Helium, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Helium3, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Iron, 100_GeV));
+ CHECK(model.isValid(Code::Proton, Code::Oxygen, 100_GeV));
// beam particles
- CHECK_THROWS(model.isValid(Code::Electron, Code::Oxygen, 100_GeV));
- CHECK_THROWS(model.isValid(Code::Iron, Code::Oxygen, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Electron, Code::Oxygen, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Iron, Code::Oxygen, 100_GeV));
// energy too low
- CHECK_THROWS(model.isValid(Code::Proton, Code::Proton, 9_GeV));
- CHECK_NOTHROW(model.isValid(Code::Proton, Code::Proton, 11_GeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Proton, 9_GeV));
+ CHECK(model.isValid(Code::Proton, Code::Proton, 11_GeV));
// energy too high
- CHECK_THROWS(model.isValid(Code::Proton, Code::Proton, 1000001_GeV));
- CHECK_NOTHROW(model.isValid(Code::Proton, Code::Proton, 999999_GeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Proton, 1000001_GeV));
+ CHECK(model.isValid(Code::Proton, Code::Proton, 999999_GeV));
// hydrogen target == proton target == neutron target
FourMomentum const aP4(100_GeV, {cs, 99_GeV, 0_GeV, 0_GeV});
@@ -262,6 +262,11 @@ TEST_CASE("SibyllInterface", "modules") {
MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
corsika::sibyll::InteractionModel hmodel;
NuclearInteractionModel model(hmodel, *env);
+
+ CHECK(model.isValid(Code::Helium, Code::Oxygen, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::PiPlus, Code::Oxygen, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Electron, Code::Oxygen, 100_GeV));
+
Code const pid = Code::Oxygen;
HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid)));
FourMomentum const P4(Elab, plab);
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 0cabc8b726a16a133f0c40e8a12f4ac0a689c23b..3b3da197bf501a9fdf73b11fcf851c39fe953182 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -75,13 +75,13 @@ TEST_CASE("UrQMD") {
SECTION("valid") {
// this is how it is currently done
- CHECK_THROWS(urqmd.isValid(Code::K0, Code::Proton));
- CHECK_THROWS(urqmd.isValid(Code::DPlus, Code::Proton));
- CHECK_THROWS(urqmd.isValid(Code::Electron, Code::Proton));
- CHECK_THROWS(urqmd.isValid(Code::Proton, Code::Electron));
- CHECK_THROWS(urqmd.isValid(Code::Oxygen, Code::Oxygen));
- CHECK_THROWS(urqmd.isValid(Code::PiPlus, Code::Omega));
- CHECK_THROWS(
+ CHECK_FALSE(urqmd.isValid(Code::K0, Code::Proton));
+ CHECK_FALSE(urqmd.isValid(Code::DPlus, Code::Proton));
+ CHECK_FALSE(urqmd.isValid(Code::Electron, Code::Proton));
+ CHECK_FALSE(urqmd.isValid(Code::Proton, Code::Electron));
+ CHECK_FALSE(urqmd.isValid(Code::Oxygen, Code::Oxygen));
+ CHECK_FALSE(urqmd.isValid(Code::PiPlus, Code::Omega));
+ CHECK_FALSE(
urqmd.isValid(Code::PiPlus, Code::Proton)); // Proton is not a valid target....
CHECK_NOTHROW(urqmd.isValid(Code::Proton, Code::Oxygen));