From 8051172de392723296032997779efcfd83f73b3b Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Thu, 28 Oct 2021 06:33:00 +0200
Subject: [PATCH] fixed P8 test
---
corsika/detail/framework/core/Cascade.inl | 6 +-
.../framework/process/ProcessSequence.inl | 1 +
.../process/SwitchProcessSequence.inl | 228 +++++++++------
corsika/detail/media/NuclearComposition.inl | 8 +-
.../detail/modules/pythia8/Interaction.inl | 275 ++++++++----------
corsika/framework/process/ProcessSequence.hpp | 52 ++--
corsika/media/IMediumPropertyModel.hpp | 3 +-
corsika/media/MediumPropertyModel.hpp | 1 -
corsika/modules/pythia8/Interaction.hpp | 14 +-
.../sibyll/NuclearInteractionModel.hpp | 1 -
corsika/stack/VectorStack.hpp | 8 +-
tests/modules/testPythia8.cpp | 14 +-
12 files changed, 314 insertions(+), 297 deletions(-)
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index b003e83d4..f3e5938fe 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -78,9 +78,10 @@ namespace corsika {
NuclearComposition const& composition =
currentLogicalNode->getModelProperties().getNuclearComposition();
- // determine sqrtS per nucleon pair, sqrtS_NN
+ // determine projectile
HEPEnergyType const Elab = particle.getEnergy();
FourMomentum const projectileP4{Elab, particle.getMomentum()};
+ // determine cross section in material
CrossSectionType const sigma =
composition.getWeightedSum([=](Code const targetId) -> CrossSectionType {
FourMomentum const targetP4(
@@ -110,7 +111,7 @@ namespace corsika {
NuclearComposition const& composition =
currentLogicalNode->getModelProperties().getNuclearComposition();
- // determine sqrtS per nucleon pair, sqrtS_NN
+ // determine projectile
HEPEnergyType const Elab = particle.getEnergy();
FourMomentum const projectileP4{Elab, particle.getMomentum()};
@@ -291,7 +292,6 @@ namespace corsika {
and Decay!
*/
if (distance_interact < distance_decay) {
- // define boost of NUCLEON-NUCLEON frame
interaction(secondaries, projectileP4, composition, total_cx);
} else {
[[maybe_unused]] auto projectile = secondaries.getProjectile();
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index e27dbc9fd..25a592d52 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -482,6 +482,7 @@ namespace corsika {
decltype(projectile) const&, // template argument
decltype(projectile) const&, // parameters
Code, FourMomentum const&>;
+
static_assert((has_signature_cx1 || has_signature_cx2),
"TProcess1 has no method with correct signature \"CrossSectionType "
"getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index fdf1e685e..0d26014fd 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -217,18 +217,39 @@ namespace corsika {
if (select_(projectile)) {
if constexpr (is_interaction_process_v<process1_type>) {
- return A_.getCrossSection(projectile.getPID(), targetId,
- {projectile.getEnergy(), projectile.getMomentum()},
- targetP4);
+ bool constexpr has_signature_cx1 =
+ has_method_getCrossSection_v<TSequence, // process object
+ CrossSectionType, // return type
+ Code, Code, // parameters
+ FourMomentum const&, FourMomentum const&>;
+ if constexpr (has_signature_cx1) {
+
+ return A_.getCrossSection(projectile.getPID(), targetId,
+ {projectile.getEnergy(), projectile.getMomentum()},
+ targetP4);
+ } else {
+ return A_.getCrossSection(projectile, projectile.getPID(),
+ {projectile.getEnergy(), projectile.getMomentum()});
+ }
} else if (process1_type::is_process_sequence) {
return A_.getCrossSection(projectile, targetId, targetP4);
}
} else {
if constexpr (is_interaction_process_v<process2_type>) {
- return B_.getCrossSection(projectile.getPID(), targetId,
- {projectile.getEnergy(), projectile.getMomentum()},
- targetP4);
+ bool constexpr has_signature_cx1 =
+ has_method_getCrossSection_v<USequence, // process object
+ CrossSectionType, // return type
+ Code, Code, // parameters
+ FourMomentum const&, FourMomentum const&>;
+ if constexpr (has_signature_cx1) {
+
+ return B_.getCrossSection(projectile.getPID(), targetId,
+ {projectile.getEnergy(), projectile.getMomentum()},
+ targetP4);
+ } else {
+ return B_.getCrossSection(projectile, targetId, targetP4);
+ }
} else if (process2_type::is_process_sequence) {
return B_.getCrossSection(projectile, targetId, targetP4);
}
@@ -259,50 +280,72 @@ namespace corsika {
// get cross section vector for all material components
// for selected process A
- static_assert(
+ bool constexpr has_signature_cx1 =
has_method_getCrossSection_v<TSequence, // process object
CrossSectionType, // return type
Code, Code, // parameters
- FourMomentum const&, FourMomentum const&>,
- "TSequence has no method with correct signature \"CrossSectionType "
- "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
- "const&)\" required by "
- "InteractionProcess<TSequence>. ");
-
- std::vector<CrossSectionType> const weightedCrossSections =
- composition.getWeighted([=](Code const targetId) -> CrossSectionType {
- FourMomentum const targetP4(
- get_mass(targetId),
- MomentumVector(projectile.getMomentum().getCoordinateSystem(),
- {0_GeV, 0_GeV, 0_GeV}));
- return A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
- });
-
- cx_sum += std::accumulate(weightedCrossSections.cbegin(),
- weightedCrossSections.cend(), CrossSectionType::zero());
- if (cx_select < cx_sum) {
+ FourMomentum const&, FourMomentum const&>;
+ bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+ has_method_getCrossSectionTemplate_v<
+ TSequence, // process object
+ CrossSectionType, // return type
+ decltype(projectile) const&, // template argument
+ decltype(projectile) const&, // parameters
+ Code, FourMomentum const&>;
+
+ static_assert((has_signature_cx1 || has_signature_cx2),
+ "TSequence has no method with correct signature \"CrossSectionType "
+ "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+ "const&)\" required by "
+ "InteractionProcess<TSequence>. ");
+
+ std::vector<CrossSectionType> weightedCrossSections;
+ if constexpr (has_signature_cx1) {
+ /*std::vector<CrossSectionType> const*/ weightedCrossSections =
+ composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+ return A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+ });
+
+ cx_sum +=
+ std::accumulate(weightedCrossSections.cbegin(),
+ weightedCrossSections.cend(), CrossSectionType::zero());
+ } else { // this is for PROPOSAL
+ cx_sum += A_.template getCrossSection(projectile, projectileId, projectileP4);
+ }
- // now also sample targetId from weighted cross sections
- Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
- FourMomentum const targetP4(
- get_mass(targetId),
- MomentumVector(projectile.getMomentum().getCoordinateSystem(),
- {0_GeV, 0_GeV, 0_GeV}));
+ if (cx_select < cx_sum) {
- // interface checking on TProcess1
- static_assert(
- has_method_doInteract_v<TSequence, // process object
- void, // return type
- TSecondaryView, // template argument
- TSecondaryView&, // method parameters
- Code, Code, FourMomentum const&,
- FourMomentum const&>,
- "TSequence has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, "
- "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
- "InteractionProcess<TSequence>. ");
-
- A_.template doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+ if constexpr (has_signature_cx1) {
+
+ // now also sample targetId from weighted cross sections
+ Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+
+ // interface checking on TProcess1
+ static_assert(
+ has_method_doInteract_v<TSequence, // process object
+ void, // return type
+ TSecondaryView, // template argument
+ TSecondaryView&, // method parameters
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>,
+ "TSequence has no method with correct signature \"void "
+ "doInteraction<TSecondaryView>(TSecondaryView&, "
+ "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+ "InteractionProcess<TSequence>. ");
+
+ A_.template doInteraction(view, projectileId, targetId, projectileP4,
+ targetP4);
+ } else { // this is for PROPOSAL
+ A_.template doInteraction(view, projectileId, projectileP4);
+ }
return ProcessReturn::Interacted;
} // end collision branch A
@@ -320,51 +363,72 @@ namespace corsika {
Code const projectileId = projectile.getPID();
// get cross section vector for all material components, for selected process B
- static_assert(has_method_getCrossSection_v<USequence, // process object
- CrossSectionType, // return type
- Code, Code, FourMomentum const&,
- FourMomentum const&>, // parameters
+ bool constexpr has_signature_cx1 =
+ has_method_getCrossSection_v<USequence, // process object
+ CrossSectionType, // return type
+ Code, Code, // parameters
+ FourMomentum const&, FourMomentum const&>;
+ bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+ has_method_getCrossSectionTemplate_v<
+ USequence, // process object
+ CrossSectionType, // return type
+ decltype(projectile) const&, // template argument
+ decltype(projectile) const&, // parameters
+ Code, FourMomentum const&>;
+
+ static_assert((has_signature_cx1 || has_signature_cx2),
"USequence has no method with correct signature \"CrossSectionType "
"getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
"const&)\" required by "
"InteractionProcess<USequence>. ");
- std::vector<CrossSectionType> const weightedCrossSections =
- composition.getWeighted([=](Code const targetId) -> CrossSectionType {
- FourMomentum const targetP4(
- get_mass(targetId),
- MomentumVector(projectile.getMomentum().getCoordinateSystem(),
- {0_GeV, 0_GeV, 0_GeV}));
- return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
- });
-
- cx_sum += std::accumulate(weightedCrossSections.begin(),
- weightedCrossSections.end(), CrossSectionType::zero());
+ std::vector<CrossSectionType> weightedCrossSections;
+ if constexpr (has_signature_cx1) {
+ /* std::vector<CrossSectionType> const*/ weightedCrossSections =
+ composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+ return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+ });
+
+ cx_sum +=
+ std::accumulate(weightedCrossSections.begin(), weightedCrossSections.end(),
+ CrossSectionType::zero());
+ } else { // this is for PROPOSAL
+ cx_sum += B_.template getCrossSection(projectile, projectileId, projectileP4);
+ }
// check if we should execute THIS process and then EXIT
if (cx_select <= cx_sum) {
- // now also sample targetId from weighted cross sections
- Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
- FourMomentum const targetP4(
- get_mass(targetId),
- MomentumVector(projectile.getMomentum().getCoordinateSystem(),
- {0_GeV, 0_GeV, 0_GeV}));
-
- // interface checking on TProcess2
- static_assert(
- has_method_doInteract_v<USequence, // process object
- void, // return type
- TSecondaryView, // template argument
- TSecondaryView&, // method parameters
- Code, Code, FourMomentum const&,
- FourMomentum const&>,
- "USequence has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, "
- "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
- "InteractionProcess<USequence>. ");
-
- B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+ if constexpr (has_signature_cx1) {
+
+ // now also sample targetId from weighted cross sections
+ Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+
+ // interface checking on TProcess2
+ static_assert(
+ has_method_doInteract_v<USequence, // process object
+ void, // return type
+ TSecondaryView, // template argument
+ TSecondaryView&, // method parameters
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>,
+ "USequence has no method with correct signature \"void "
+ "doInteraction<TSecondaryView>(TSecondaryView&, "
+ "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+ "InteractionProcess<USequence>. ");
+
+ B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+ } else { // this is for PROPOSAL
+ B_.doInteraction(view, projectileId, projectileP4);
+ }
return ProcessReturn::Interacted;
} // end collision in branch B
@@ -372,7 +436,7 @@ namespace corsika {
} // end branch B_
return ProcessReturn::Ok;
- }
+ } // namespace corsika
template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
diff --git a/corsika/detail/media/NuclearComposition.inl b/corsika/detail/media/NuclearComposition.inl
index 0237d73e9..c499a593d 100644
--- a/corsika/detail/media/NuclearComposition.inl
+++ b/corsika/detail/media/NuclearComposition.inl
@@ -109,12 +109,8 @@ namespace corsika {
}
std::discrete_distribution channelDist(
- WeightProviderIterator //<decltype(numberFractions_.begin()),
- // decltype(sigma.begin())>
- (numberFractions_.begin(), sigma.begin()),
- WeightProviderIterator //<decltype(numberFractions_.begin()),
- // decltype(sigma.end())>
- (numberFractions_.end(), sigma.end()));
+ WeightProviderIterator(numberFractions_.begin(), sigma.begin()),
+ WeightProviderIterator(numberFractions_.end(), sigma.end()));
auto const iChannel = channelDist(randomStream);
return components_[iChannel];
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 5c495e693..a716b156d 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -80,6 +80,28 @@ namespace corsika::pythia8 {
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
}
+ inline void 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 (targetId != Code::Hydrogen && targetId != Code::Neutron &&
+ targetId != Code::Proton) {
+ throw std::runtime_error("wrong target for PYTHIA");
+ }
+
+ 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 (!canInteract(projectileId)) {
+ throw std::runtime_error("Projectile not supported by PYTHIA!");
+ }
+ }
+
inline void Interaction::configureLabFrameCollision(Code const projectileId,
Code const targetId,
HEPEnergyType const BeamEnergy) {
@@ -128,180 +150,117 @@ namespace corsika::pythia8 {
HEPEnergyType const CoMenergy = (projectileP4 + targetP4).getNorm();
- // interaction possible in pythia?
- if (targetId == Code::Proton || targetId == Code::Hydrogen) {
- if (canInteract(projectileId) && isValidCoMEnergy(CoMenergy)) {
- // input particle PDG
- auto const pdgCodeBeam = static_cast<int>(get_PDG(projectileId));
- auto const pdgCodeTarget = static_cast<int>(get_PDG(targetId));
- double const ecm = CoMenergy / 1_GeV;
-
- //! @todo: remove this const_cast, when Pythia8 becomes const-correct! CHECK!
- Pythia8::SigmaTotal& sigma = *const_cast<Pythia8::SigmaTotal*>(&sigma_);
-
- // calculate cross section
- sigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
- if (sigma.hasSigmaTot()) {
- double const sigEla = sigma.sigmaEl();
- double const sigProd = sigma.sigmaTot() - sigEla;
-
- return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
-
- } else {
- // we can't test pythia8 internals, LCOV_EXCL_START
- throw std::runtime_error("pythia cross section init failed");
- // we can't test pythia8 internals, LCOV_EXCL_STOP
- }
- } else {
- return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
- std::numeric_limits<double>::infinity() * 1_mb);
- }
+ isValid(projectileId, targetId, CoMenergy); // throws
+
+ // input particle PDG
+ auto const pdgCodeBeam = static_cast<int>(get_PDG(projectileId));
+ auto const pdgCodeTarget = static_cast<int>(get_PDG(targetId));
+ double const ecm = CoMenergy / 1_GeV;
+
+ //! @todo: remove this const_cast, when Pythia8 becomes const-correct! CHECK!
+ Pythia8::SigmaTotal& sigma = *const_cast<Pythia8::SigmaTotal*>(&sigma_);
+
+ // calculate cross section
+ sigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
+ if (sigma.hasSigmaTot()) {
+ double const sigEla = sigma.sigmaEl();
+ double const sigProd = sigma.sigmaTot() - sigEla;
+
+ return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
+
} else {
- throw std::runtime_error("invalid target for pythia");
+ // we can't test pythia8 internals, LCOV_EXCL_START
+ throw std::runtime_error("pythia cross section init failed");
+ // we can't test pythia8 internals, LCOV_EXCL_STOP
}
}
template <class TView>
inline void Interaction::doInteraction(TView& view, Code const projectileId,
- Code const targetId, FourMomentum const&,
- FourMomentum const&) {
+ Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) {
auto projectile = view.getProjectile();
CORSIKA_LOG_DEBUG(
"Pythia::Interaction: "
- "DoInteraction: {} interaction? ",
+ "doInteraction: {} interaction? ",
projectileId, corsika::pythia8::Interaction::canInteract(projectileId));
- if (is_nucleus(projectileId)) {
- // nuclei handled by different process, this should not happen
- throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
+ // define system
+ auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+ HEPEnergyType const sqrtS = sqrt(sqrtS2);
+ HEPEnergyType const eProjectileLab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+ static_pow<2>(get_mass(targetId))) /
+ (2 * get_mass(targetId));
+
+ isValid(projectileId, targetId, sqrtS); // throws
+
+ // position and time of interaction
+ Point pOrig = projectile.getPosition();
+ TimeType tOrig = projectile.getTime();
+
+ CORSIKA_LOG_DEBUG("Interaction: ebeam lab: {} GeV", eProjectileLab / 1_GeV);
+
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Pythia projectile: +z
+ COMBoost const boost(projectileP4, constants::nucleonMass);
+ auto const& labCS = boost.getOriginalCS();
+
+ CORSIKA_LOG_DEBUG("Interaction: position of interaction: ", pOrig.getCoordinates());
+ CORSIKA_LOG_DEBUG("Interaction: time: {}", tOrig);
+
+ CORSIKA_LOG_DEBUG(
+ "Interaction: "
+ " doInteraction: E(GeV): {}"
+ " Ecm(GeV): {}",
+ eProjectileLab / 1_GeV, sqrtS / 1_GeV);
+
+ count_++;
+
+ configureLabFrameCollision(projectileId, targetId, eProjectileLab);
+
+ // create event in pytia. LCOV_EXCL_START: we don't validate pythia8 internals
+ if (!Pythia8::Pythia::next())
+ throw std::runtime_error("Pythia::DoInteraction: failed!");
+ // LCOV_EXCL_STOP
+
+ // link to pythia stack
+ Pythia8::Event& event = Pythia8::Pythia::event;
+
+ // LCOV_EXCL_START, we don't validate pythia8 internals
+ if (print_listing_) {
+ // print final state
+ event.list();
}
+ // LCOV_EXCL_STOP
+
+ MomentumVector Plab_final(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV;
+ for (int i = 0; i < event.size(); ++i) {
+ Pythia8::Particle const& p8p = event[i];
+ // skip particles that have decayed in pythia
+ if (!p8p.isFinal()) continue;
- if (corsika::pythia8::Interaction::canInteract(projectileId)) {
-
- // define projectile
- HEPEnergyType const eProjectileLab = projectile.getEnergy();
- auto const pProjectileLab = projectile.getMomentum();
- CoordinateSystemPtr const& labCS = pProjectileLab.getCoordinateSystem();
-
- // position and time of interaction, not used in Sibyll
- Point pOrig = projectile.getPosition();
- TimeType tOrig = projectile.getTime();
-
- // define target
- // FOR NOW: target is always at rest
- auto const eTargetLab = 0_GeV + constants::nucleonMass;
- auto const pTargetLab = MomentumVector(labCS, 0_GeV, 0_GeV, 0_GeV);
- FourVector const PtargLab(eTargetLab, pTargetLab);
-
- CORSIKA_LOG_DEBUG(
- "Interaction: ebeam lab: {} GeV"
- "Interaction: pbeam lab: {} GeV",
- eProjectileLab / 1_GeV, pProjectileLab.getComponents() / 1_GeV);
-
- CORSIKA_LOG_DEBUG(
- "Interaction: etarget lab: {} GeV"
- "Interaction: ptarget lab: {} GeV ",
- eTargetLab / 1_GeV, pTargetLab.getComponents() / 1_GeV);
-
- FourVector const PprojLab(eProjectileLab, pProjectileLab);
-
- // define target kinematics in lab frame
- // define boost to and from CoM frame
- // CoM frame definition in Pythia projectile: +z
- COMBoost const boost(PprojLab, constants::nucleonMass);
-
- // just for show:
- // boost projecticle
- auto const PprojCoM = boost.toCoM(PprojLab);
-
- // boost target
- auto const PtargCoM = boost.toCoM(PtargLab);
-
- CORSIKA_LOG_DEBUG(
- "Interaction: ebeam CoM: {} GeV"
- "Interaction: pbeam CoM: {} GeV",
- PprojCoM.getTimeLikeComponent() / 1_GeV,
- PprojCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
-
- CORSIKA_LOG_DEBUG(
- "Interaction: etarget CoM: {} GeV"
- "Interaction: ptarget CoM: {} GeV",
- PtargCoM.getTimeLikeComponent() / 1_GeV,
- PtargCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
-
- CORSIKA_LOG_DEBUG("Interaction: position of interaction: ", pOrig.getCoordinates());
- CORSIKA_LOG_DEBUG("Interaction: time: {}", tOrig);
-
- HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = projectile.getMomentum();
- // invariant mass, i.e. cm. energy
- HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.getSquaredNorm());
-
- if (targetId != Code::Hydrogen && targetId != Code::Neutron &&
- targetId != Code::Proton)
- throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
-
- CORSIKA_LOG_DEBUG(
- "Interaction: "
- " DoInteraction: E(GeV): {}"
- " Ecm(GeV): {}",
- eProjectileLab / 1_GeV, Ecm / 1_GeV);
-
- if (eProjectileLab < 8.5_GeV || !isValidCoMEnergy(Ecm)) {
- CORSIKA_LOG_DEBUG(
- "Interaction: "
- " DoInteraction: should have dropped particle.. "
- "THIS IS AN ERROR");
- throw std::runtime_error("energy too low for PYTHIA");
-
- } else {
- count_++;
-
- configureLabFrameCollision(projectileId, targetId, eProjectileLab);
-
- // create event in pytia. LCOV_EXCL_START: we don't validate pythia8 internals
- if (!Pythia8::Pythia::next())
- throw std::runtime_error("Pythia::DoInteraction: failed!");
- // LCOV_EXCL_STOP
-
- // link to pythia stack
- Pythia8::Event& event = Pythia8::Pythia::event;
-
- // LCOV_EXCL_START, we don't validate pythia8 internals
- if (print_listing_) {
- // print final state
- event.list();
- }
- // LCOV_EXCL_STOP
-
- MomentumVector Plab_final(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV;
- for (int i = 0; i < event.size(); ++i) {
- Pythia8::Particle& p8p = event[i];
- // skip particles that have decayed in pythia
- if (!p8p.isFinal()) continue;
-
- auto const pyId = convert_from_PDG(static_cast<PDGCode>(p8p.id()));
-
- MomentumVector const pyPlab(
- labCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
-
- // add to corsika stack
- auto pnew =
- projectile.addSecondary(std::make_tuple(pyId, pyPlab, pOrig, tOrig));
-
- Plab_final += pnew.getMomentum();
- Elab_final += pnew.getEnergy();
- }
- CORSIKA_LOG_DEBUG(
- "conservation (all GeV): "
- "Elab_final= {}"
- ", Plab_final= {}",
- Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
- }
+ auto const pyId = convert_from_PDG(static_cast<PDGCode>(p8p.id()));
+
+ MomentumVector const pyPlab(labCS,
+ {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+
+ // add to corsika stack
+ auto pnew = projectile.addSecondary(std::make_tuple(pyId, pyPlab, pOrig, tOrig));
+
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
}
+
+ CORSIKA_LOG_DEBUG(
+ "conservation (all GeV): "
+ "Elab_final= {}"
+ ", Plab_final= {}",
+ Elab_final / 1_GeV, (Plab_final / 1_GeV).getComponents());
}
} // namespace corsika::pythia8
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 3b647ec82..74dad4f2f 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -141,7 +141,7 @@ namespace corsika {
* the static version of the
* ProcessSequence and all Process
* types are based on the CRTP C++
- * design pattern)
+ * design pattern).
*
* Template parameters:
* @tparam TProcess1 is of type BaseProcess, either a dedicatd process, or a
@@ -177,7 +177,7 @@ namespace corsika {
static bool const is_process_sequence = true;
/**
- * Only valid user constructor will create fully initialized object
+ * Only valid user constructor will create fully initialized object.
*
* ProcessSequence supports and encourages move semantics. You can
* use object, l-value references or r-value references to
@@ -339,10 +339,10 @@ namespace corsika {
* @fn make_sequence
* @ingroup Processes
*
- * Factory function to create a ProcessSequence
+ * Factory function to create a ProcessSequence.
*
* to construct ProcessSequences in a flexible and dynamic way the
- * `sequence` factory functions are provided
+ * `sequence` factory functions are provided.
*
* Any objects of type
* - BaseProcess
@@ -373,34 +373,34 @@ namespace corsika {
}
/**
- * @fn make_sequence
- * @ingroup Processes
- *
- * Factory function to create ProcessSequence
- *
- specialization for two input objects (no paramter pack in vB).
-
- @tparam TProcess1 another BaseProcess
- @tparam TProcess2 another BaseProcess
- @param vA needs to derive from BaseProcess
- @param vB needs to derive BaseProcess
- */
+ * @fn make_sequence
+ * @ingroup Processes
+ *
+ * Factory function to create ProcessSequence.
+ *
+ * specialization for two input objects (no paramter pack in vB).
+ *
+ * @tparam TProcess1 another BaseProcess
+ * @tparam TProcess2 another BaseProcess
+ * @param vA needs to derive from BaseProcess
+ * @param vB needs to derive BaseProcess
+ */
template <typename TProcess1, typename TProcess2>
ProcessSequence<TProcess1, TProcess2> make_sequence(TProcess1&& vA, TProcess2&& vB) {
return ProcessSequence<TProcess1, TProcess2>(vA, vB);
}
/**
- @fn make_sequence
- @ingroup Processes
-
- Factory function to create ProcessSequence from a single BaseProcess
-
- also allow a single Process in ProcessSequence, accompany by
- `NullModel`
-
- @tparam TProcess1 another BaseProcess
- @param vA needs to derive from BaseProcess
+ * @fn make_sequence
+ * @ingroup Processes
+ *
+ * Factory function to create ProcessSequence from a single BaseProcess.
+ *
+ * also allow a single Process in ProcessSequence, accompany by
+ * `NullModel`.
+ *
+ * @tparam TProcess1 another BaseProcess
+ * @param vA needs to derive from BaseProcess
*/
template <typename TProcess>
ProcessSequence<TProcess, NullModel> make_sequence(TProcess&& vA) {
diff --git a/corsika/media/IMediumPropertyModel.hpp b/corsika/media/IMediumPropertyModel.hpp
index 860fa22da..94e80fb51 100644
--- a/corsika/media/IMediumPropertyModel.hpp
+++ b/corsika/media/IMediumPropertyModel.hpp
@@ -16,10 +16,9 @@
namespace corsika {
/**
- * An interface for type of media, needed e.g. to determine energy losses
+ * An interface for type of media, needed e.g. to determine energy losses.
*
* This is the base interface for media types.
- *
*/
template <typename TModel>
class IMediumPropertyModel : public TModel {
diff --git a/corsika/media/MediumPropertyModel.hpp b/corsika/media/MediumPropertyModel.hpp
index ab89a8a54..97bce6ff8 100644
--- a/corsika/media/MediumPropertyModel.hpp
+++ b/corsika/media/MediumPropertyModel.hpp
@@ -14,7 +14,6 @@ namespace corsika {
/**
* A model for the energy loss property of a medium.
- *
*/
template <typename T>
class MediumPropertyModel : public T {
diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp
index 63376ab17..02a5313b7 100644
--- a/corsika/modules/pythia8/Interaction.hpp
+++ b/corsika/modules/pythia8/Interaction.hpp
@@ -21,20 +21,22 @@ namespace corsika::pythia8 {
class Interaction : public InteractionProcess<Interaction>, public Pythia8::Pythia {
public:
- Interaction(const bool print_listing = false);
+ Interaction(bool const print_listing = false);
~Interaction();
void setStable(std::vector<Code> const&);
- void setUnstable(const Code);
- void setStable(const Code);
+ void setUnstable(Code const);
+ void setStable(Code const);
bool isValidCoMEnergy(HEPEnergyType const ecm) const {
return (10_GeV < ecm) && (ecm < 1_PeV);
}
- bool canInteract(const Code) const;
- void configureLabFrameCollision(const Code, const Code, const HEPEnergyType);
+ bool canInteract(Code const) const;
+ void configureLabFrameCollision(Code const, Code const, HEPEnergyType const);
+ void isValid(Code const projectileId, Code const targetId,
+ HEPEnergyType const sqrtS) const;
/**
* Returns inelastic AND elastic cross sections.
*
@@ -88,7 +90,7 @@ namespace corsika::pythia8 {
private:
default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("pythia");
Pythia8::SigmaTotal sigma_;
- const bool internalDecays_ = true;
+ bool const internalDecays_ = true;
int count_ = 0;
bool print_listing_ = false;
};
diff --git a/corsika/modules/sibyll/NuclearInteractionModel.hpp b/corsika/modules/sibyll/NuclearInteractionModel.hpp
index 1c4d9c7e9..753cabf98 100644
--- a/corsika/modules/sibyll/NuclearInteractionModel.hpp
+++ b/corsika/modules/sibyll/NuclearInteractionModel.hpp
@@ -10,7 +10,6 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
namespace corsika::sibyll {
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index df1b020a2..22873790c 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -66,7 +66,6 @@ namespace corsika {
* Set data of new particle.
*
* @param v tuple containing: PID, kinetic Energy, Direction Vector, Position, Time
- *
*/
void setParticleData(particle_data_momentum_type const& v);
@@ -75,7 +74,6 @@ namespace corsika {
*
* @param parent parent particle
* @param v tuple containing: PID, kinetic Energy, Direction Vector, Position, Time
- *
*/
void setParticleData(ParticleInterface<TStackIterator> const& parent,
particle_data_momentum_type const& v);
@@ -97,9 +95,9 @@ namespace corsika {
}
/**
- The MomentumVector v is used to determine the DirectionVector, and to update the
- particle energy.
- */
+ * The MomentumVector v is used to determine the DirectionVector, and to update the
+ * particle energy.
+ */
void setMomentum(MomentumVector const& v) {
HEPMomentumType const P = v.getNorm();
if (P == 0_eV) {
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index 2ccc5407d..06121eeec 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -183,17 +183,17 @@ TEST_CASE("Pythia8Interface", "modules") {
std::tuple<CrossSectionType, CrossSectionType> xs_test =
collision.getCrossSectionInelEla(
Code::Proton, Code::Hydrogen,
- {sqrt(static_pow<2>(Iron::mass) + static_pow<2>(100_GeV)),
+ {sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
{rootCS, {0_eV, 0_eV, 100_GeV}}},
{Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
- CHECK(std::get<0>(xs_test) / 1_mb == Approx(353).margin(2));
- CHECK(std::get<1>(xs_test) / 1_mb == Approx(75).margin(2));
+ CHECK(std::get<0>(xs_test) / 1_mb == Approx(314).margin(2));
+ CHECK(std::get<1>(xs_test) / 1_mb == Approx(69).margin(2));
collision.doInteraction(view, Code::Proton, Code::Hydrogen,
- {sqrt(static_pow<2>(Iron::mass) + static_pow<2>(100_GeV)),
+ {sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
{rootCS, {0_eV, 0_eV, 100_GeV}}},
{Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
- CHECK(view.getSize() == 38);
+ CHECK(view.getSize() == 12);
}
SECTION("pythia too low energy") {
@@ -207,8 +207,8 @@ TEST_CASE("Pythia8Interface", "modules") {
CHECK_THROWS(collision.doInteraction(
view, Code::Neutron, Code::Hydrogen,
- {sqrt(static_pow<2>(Neutron::mass) + static_pow<2>(100_GeV)),
- {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {sqrt(static_pow<2>(Neutron::mass) + static_pow<2>(1_GeV)),
+ {rootCS, {0_eV, 0_eV, 1_GeV}}},
{Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
}
--
GitLab