From 2af74c694d1ea2603a13bed7cfe96615593a12eb Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 19 Oct 2021 00:01:33 +0200
Subject: [PATCH] sibyll works...
---
corsika/detail/framework/core/Cascade.inl | 39 ++++----
.../detail/framework/geometry/FourVector.inl | 22 ++---
.../framework/process/ProcessSequence.inl | 77 +++++++++------
.../process/SwitchProcessSequence.inl | 98 ++++++++++---------
corsika/detail/framework/utility/COMBoost.inl | 52 +++++++---
.../modules/sibyll/InteractionModel.inl | 31 +++---
.../sibyll/NuclearInteractionModel.inl | 80 ++++++++-------
corsika/framework/core/Cascade.hpp | 10 +-
corsika/framework/geometry/FourVector.hpp | 83 +++++++++-------
corsika/framework/process/ProcessSequence.hpp | 33 ++++++-
.../process/SwitchProcessSequence.hpp | 4 +-
corsika/framework/utility/COMBoost.hpp | 56 ++++++-----
corsika/modules/sibyll/InteractionModel.hpp | 13 ++-
.../sibyll/NuclearInteractionModel.hpp | 8 +-
tests/framework/testCascade.cpp | 10 +-
tests/framework/testProcessSequence.cpp | 44 +++++----
tests/modules/testSibyll.cpp | 37 ++++---
17 files changed, 404 insertions(+), 293 deletions(-)
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index ae11cf09a..0e00520bb 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -80,18 +80,18 @@ namespace corsika {
// determine sqrtS per nucleon pair, sqrtS_NN
Code const projectileId = particle.getPID();
- unsigned int const projectileA =
- (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1);
- HEPEnergyType const ElabNN = particle.getEnergy() / projectileA;
- HEPEnergyType const sqrtSnn = sqrt(2 * ElabNN * constants::nucleonMass);
-
- COMBoost const boost{{ElabNN, particle.getMomentum() / projectileA},
- constants::nucleonMass};
+ HEPEnergyType const Elab = particle.getEnergy();
+ FourMomentum const projectileP4{Elab, particle.getMomentum()};
CrossSectionType const sigma =
composition.getWeightedSum([=](Code const targetId) -> CrossSectionType {
- return sequence_.getCrossSection(projectileId, targetId, sqrtSnn);
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(particle.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+ return sequence_.getCrossSection(projectileId, targetId, projectileP4,
+ targetP4);
});
- interaction(secondaries, boost, sqrtSnn, composition, sigma);
+ interaction(secondaries, projectileP4, composition, sigma);
sequence_.doSecondaries(secondaries);
particle.erase(); // primary particle is done
}
@@ -114,16 +114,19 @@ namespace corsika {
// determine sqrtS per nucleon pair, sqrtS_NN
Code const projectileId = particle.getPID();
- unsigned int const projectileA =
- (is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1);
- HEPEnergyType const ElabNN = particle.getEnergy() / projectileA;
- HEPEnergyType const sqrtSnn = sqrt(2 * ElabNN * constants::nucleonMass);
+ HEPEnergyType const Elab = particle.getEnergy();
+ FourMomentum const projectileP4{Elab, particle.getMomentum()};
// determine combined full inelastic cross section of the particles in the material
CrossSectionType const total_cx =
composition.getWeightedSum([=](Code const targetId) -> CrossSectionType {
- return sequence_.getCrossSection(projectileId, targetId, sqrtSnn);
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(particle.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+ return sequence_.getCrossSection(projectileId, targetId, projectileP4,
+ targetP4);
});
// calculate interaction length in medium
@@ -293,9 +296,7 @@ namespace corsika {
*/
if (distance_interact < distance_decay) {
// define boost of NUCLEON-NUCLEON frame
- COMBoost const boost({ElabNN, particle.getMomentum() / projectileA},
- constants::nucleonMass);
- interaction(secondaries, boost, sqrtSnn, composition, total_cx);
+ interaction(secondaries, projectileP4, composition, total_cx);
} else {
[[maybe_unused]] auto projectile = secondaries.getProjectile();
@@ -334,7 +335,7 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
inline ProcessReturn Cascade<TTracking, TProcessList, TOutput, TStack>::interaction(
- stack_view_type& view, COMBoost const& boost, HEPEnergyType const sqrtSnn,
+ stack_view_type& view, FourMomentum const& projectileP4,
NuclearComposition const& composition,
CrossSectionType const initial_cross_section) {
@@ -346,7 +347,7 @@ namespace corsika {
UniformRealDistribution<CrossSectionType> uniDist(initial_cross_section);
CrossSectionType const sample_process_by_cx = uniDist(rng_);
- auto const returnCode = sequence_.selectInteraction(view, boost, sqrtSnn, composition,
+ auto const returnCode = sequence_.selectInteraction(view, projectileP4, composition,
rng_, sample_process_by_cx);
if (returnCode != ProcessReturn::Interacted) {
CORSIKA_LOG_DEBUG("Particle did not interact!");
diff --git a/corsika/detail/framework/geometry/FourVector.inl b/corsika/detail/framework/geometry/FourVector.inl
index 5d871430d..62542aae5 100644
--- a/corsika/detail/framework/geometry/FourVector.inl
+++ b/corsika/detail/framework/geometry/FourVector.inl
@@ -54,8 +54,8 @@ namespace corsika {
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
- operator+=(FourVector const& b) {
+ inline FourVector<TTimeType, TSpaceVecType>&
+ FourVector<TTimeType, TSpaceVecType>::operator+=(FourVector const& b) {
timeLike_ += b.timeLike_;
spaceLike_ += b.spaceLike_;
@@ -63,39 +63,39 @@ namespace corsika {
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
- operator-=(FourVector const& b) {
+ inline FourVector<TTimeType, TSpaceVecType>&
+ FourVector<TTimeType, TSpaceVecType>::operator-=(FourVector const& b) {
timeLike_ -= b.timeLike_;
spaceLike_ -= b.spaceLike_;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
- operator*=(double const b) {
+ inline FourVector<TTimeType, TSpaceVecType>&
+ FourVector<TTimeType, TSpaceVecType>::operator*=(double const b) {
timeLike_ *= b;
spaceLike_ *= b;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
- operator/=(double const b) {
+ inline FourVector<TTimeType, TSpaceVecType>&
+ FourVector<TTimeType, TSpaceVecType>::operator/=(double const b) {
timeLike_ /= b;
spaceLike_.getComponents() /= b;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
- inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
- operator/(double const b) {
+ inline FourVector<TTimeType, TSpaceVecType>&
+ FourVector<TTimeType, TSpaceVecType>::operator/(double const b) {
*this /= b;
return *this;
}
template <typename TTimeType, typename TSpaceVecType>
inline typename FourVector<TTimeType, TSpaceVecType>::norm_type
- FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
+ FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
second)>::value)
return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index e7718f885..fb977f252 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -314,7 +314,8 @@ namespace corsika {
ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
IndexProcess2>::getCrossSection(Code const projectileId,
Code const targetId,
- HEPEnergyType const sqrtSnn) const {
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
CrossSectionType tot = CrossSectionType::zero();
@@ -322,14 +323,14 @@ namespace corsika {
// errors if process1_type is invalid
if constexpr (is_interaction_process_v<process1_type> ||
process1_type::is_process_sequence) {
- tot += A_.getCrossSection(projectileId, targetId, sqrtSnn);
+ tot += A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
}
}
if constexpr (is_process_v<process2_type>) { // to protect from further compiler
// errors if process2_type is invalid
if constexpr (is_interaction_process_v<process2_type> ||
process2_type::is_process_sequence) {
- tot += B_.getCrossSection(projectileId, targetId, sqrtSnn);
+ tot += B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
}
}
return tot;
@@ -413,8 +414,7 @@ namespace corsika {
template <typename TSecondaryView>
inline ProcessReturn
ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
- selectInteraction(TSecondaryView&& view, COMBoost const& boost,
- [[maybe_unused]] HEPEnergyType const sqrtSnn,
+ selectInteraction(TSecondaryView&& view, FourMomentum const& projectileP4,
[[maybe_unused]] NuclearComposition const& composition,
[[maybe_unused]] TRNG&& rng,
[[maybe_unused]] CrossSectionType const cx_select,
@@ -426,8 +426,8 @@ namespace corsika {
// errors if process1_type is invalid
if constexpr (process1_type::is_process_sequence) {
// if A is a process sequence --> check inside
- ProcessReturn const ret = A_.selectInteraction(view, boost, sqrtSnn, composition,
- rng, cx_select, cx_sum);
+ ProcessReturn const ret =
+ A_.selectInteraction(view, projectileP4, composition, rng, cx_select, cx_sum);
// if A_ did succeed, stop routine. Not checking other static branch B_.
if (ret != ProcessReturn::Ok) { return ret; }
} else if constexpr (is_interaction_process_v<process1_type>) {
@@ -436,17 +436,23 @@ namespace corsika {
Code const projectileId = projectile.getPID();
// get cross section vector for all material components
- static_assert(has_method_getCrossSection_v<TProcess1, // process object
- CrossSectionType, // return type
- Code, // parameters
- Code, HEPEnergyType>,
- "TProcess1 has no method with correct signature \"CrossSectionType "
- "getCrossSection(Code, Code, HEPEnergyType)\" required by "
- "InteractionProcess<TProcess1>. ");
+ static_assert(
+ has_method_getCrossSection_v<TProcess1, // process object
+ CrossSectionType, // return type
+ Code, Code, // parameters
+ FourMomentum const&, FourMomentum const&>,
+ "TProcess1 has no method with correct signature \"CrossSectionType "
+ "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+ "const&)\" required by "
+ "InteractionProcess<TProcess1>. ");
std::vector<CrossSectionType> const weightedCrossSections =
composition.getWeighted([=](Code const targetId) -> CrossSectionType {
- return A_.getCrossSection(projectileId, targetId, sqrtSnn);
+ 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(),
@@ -457,6 +463,10 @@ namespace corsika {
// 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(
@@ -464,13 +474,14 @@ namespace corsika {
void, // return type
TSecondaryView, // template argument
TSecondaryView&, // method parameters
- COMBoost const&, Code, Code, HEPEnergyType>,
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>,
"TProcess1 has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, COMBoost&, Code, "
- "Code, HEPEnergyType)\" required for "
+ "doInteraction<TSecondaryView>(TSecondaryView&, "
+ "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
"InteractionProcess<TProcess1>. ");
- A_.template doInteraction(view, boost, projectileId, targetId, sqrtSnn);
+ A_.template doInteraction(view, projectileId, targetId, projectileP4, targetP4);
return ProcessReturn::Interacted;
}
@@ -482,7 +493,7 @@ namespace corsika {
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
- return B_.selectInteraction(view, boost, sqrtSnn, composition, rng, cx_select,
+ return B_.selectInteraction(view, projectileP4, composition, rng, cx_select,
cx_sum);
} else if constexpr (is_interaction_process_v<process2_type>) {
@@ -492,15 +503,20 @@ namespace corsika {
// get cross section vector for all material components
static_assert(has_method_getCrossSection_v<TProcess2, // process object
CrossSectionType, // return type
- Code, // parameters
- Code, HEPEnergyType>,
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>, // parameters
"TProcess2 has no method with correct signature \"CrossSectionType "
- "getCrossSection(Code, Code, HEPEnergyType)\" required by "
+ "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+ "const&)\" required by "
"InteractionProcess<TProcess1>. ");
std::vector<CrossSectionType> const weightedCrossSections =
composition.getWeighted([=](Code const targetId) -> CrossSectionType {
- return B_.getCrossSection(projectileId, targetId, sqrtSnn);
+ 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(),
@@ -511,6 +527,10 @@ namespace corsika {
// 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(
@@ -518,13 +538,14 @@ namespace corsika {
void, // return type
TSecondaryView, // template argument
TSecondaryView&, // method parameters
- COMBoost const&, Code, Code, HEPEnergyType>,
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>,
"TProcess1 has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, COMBoost&, Code, "
- "Code, HEPEnergyType)\" required for "
+ "doInteraction<TSecondaryView>(TSecondaryView&, "
+ "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
"InteractionProcess<TProcess2>. ");
- B_.doInteraction(view, boost, projectileId, targetId, sqrtSnn);
+ B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
return ProcessReturn::Interacted;
}
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index 46e9bfbd2..b7cf27db0 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -243,8 +243,8 @@ namespace corsika {
template <typename TSecondaryView, typename TRNG>
inline ProcessReturn SwitchProcessSequence<
TCondition, TSequence, USequence, IndexStart, IndexProcess1,
- IndexProcess2>::selectInteraction(TSecondaryView& view, COMBoost const& boost,
- HEPEnergyType const sqrtSnn,
+ IndexProcess2>::selectInteraction(TSecondaryView& view,
+ FourMomentum const& projectileP4,
NuclearComposition const& composition, TRNG& rng,
[[maybe_unused]] CrossSectionType const cx_select,
[[maybe_unused]] CrossSectionType cx_sum) {
@@ -252,7 +252,7 @@ namespace corsika {
if (select_(view.parent())) {
if constexpr (process1_type::is_process_sequence) {
// if A_ is a process sequence --> check inside
- return A_.selectInteraction(view, boost, sqrtSnn, composition, rng, cx_select,
+ return A_.selectInteraction(view, projectileP4, composition, rng, cx_select,
cx_sum);
} else if constexpr (is_interaction_process_v<process1_type>) {
@@ -265,16 +265,18 @@ namespace corsika {
has_method_getCrossSection_v<TSequence, // process object
CrossSectionType, // return type
Code, // parameters
- Code, HEPEnergyType, unsigned int, unsigned int>,
+ Code, FourMomentum const&, FourMomentum const&>,
"TSequence has no method with correct signature \"CrossSectionType "
- "getCrossSection(Code, Code, HEPEnergyType, unsigned int, unsigned int"
- ")\" required by InteractionProcess<TSequence>. ");
+ "getCrossSection(Code, Code, FourMomentum const&, FourMomntum const&)\" "
+ "required by InteractionProcess<TSequence>. ");
std::vector<CrossSectionType> const weightedCrossSections =
composition.getWeighted([=](Code const targetId) -> CrossSectionType {
- return A_.getCrossSection(
- projectileId, targetId, sqrtSnn, projectileA,
- is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ 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(),
@@ -283,23 +285,24 @@ namespace corsika {
// 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
- COMBoost const&, Code, Code, HEPEnergyType,
- unsigned int, unsigned int>,
- "USequence has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, COMBoost&, Code, "
- "Code, HEPEnergyType, unsigned int, unsigned int)\" required for "
- "InteractionProcess<USequence>. ");
-
- A_.template doInteraction(view, boost, projectileId, targetId, sqrtSnn,
- projectileA,
- is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ static_assert(has_method_doInteract_v<TSequence, // 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&, FourMomntum const&)\" required for "
+ "InteractionProcess<USequence>. ");
+
+ A_.template doInteraction(view, projectileId, targetId, projectileP4, targetP4);
return ProcessReturn::Interacted;
@@ -310,7 +313,7 @@ namespace corsika {
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
- return B_.selectInteraction(view, boost, sqrtSnn, composition, rng, cx_select,
+ return B_.selectInteraction(view, projectileP4, composition, rng, cx_select,
cx_sum);
} else if constexpr (is_interaction_process_v<process2_type>) {
@@ -323,16 +326,18 @@ namespace corsika {
has_method_getCrossSection_v<USequence, // process object
CrossSectionType, // return type
Code, // parameters
- Code, HEPEnergyType, unsigned int, unsigned int>,
+ Code, FourMomentum const&, FourMomentum const&>,
"USequence has no method with correct signature \"CrossSectionType "
- "getCrossSection(Code, Code, HEPEnergyType, unsigned int, unsigned int"
- ")\" required by InteractionProcess<USequence>. ");
+ "getCrossSection(Code, Code, FourMomentum const&, FourMomentum const&)\" "
+ "required by InteractionProcess<USequence>. ");
std::vector<CrossSectionType> const weightedCrossSections =
composition.getWeighted([=](Code const targetId) -> CrossSectionType {
- return B_.getCrossSection(
- projectileId, targetId, sqrtSnn, projectileA,
- is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ 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.cbegin(),
@@ -342,23 +347,24 @@ namespace corsika {
// 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<USequence, // process object
- void, // return type
- TSecondaryView, // template argument
- TSecondaryView&, // method parameters
- COMBoost const&, Code, Code, HEPEnergyType,
- unsigned int, unsigned int>,
- "USequence has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, COMBoost&, Code, "
- "Code, HEPEnergyType, unsigned int, unsigned int)\" required for "
- "InteractionProcess<USequence>. ");
-
- B_.template doInteraction(view, boost, projectileId, targetId, sqrtSnn,
- projectileA,
- is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ 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_.template doInteraction(view, projectileId, targetId, projectileP4, targetP4);
return ProcessReturn::Interacted;
} // end collision in branch B
diff --git a/corsika/detail/framework/utility/COMBoost.inl b/corsika/detail/framework/utility/COMBoost.inl
index 842a6f99e..7551b6e0a 100644
--- a/corsika/detail/framework/utility/COMBoost.inl
+++ b/corsika/detail/framework/utility/COMBoost.inl
@@ -19,17 +19,17 @@
namespace corsika {
- inline COMBoost::COMBoost(FourVector<HEPEnergyType, MomentumVector> const& Pprojectile,
+ inline COMBoost::COMBoost(FourMomentum const& P4projectile,
HEPMassType const massTarget)
- : originalCS_{Pprojectile.getSpaceLikeComponents().getCoordinateSystem()}
+ : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()}
, rotatedCS_{
- make_rotationToZ(Pprojectile.getSpaceLikeComponents().getCoordinateSystem(),
- Pprojectile.getSpaceLikeComponents())} {
- auto const pProjectile = Pprojectile.getSpaceLikeComponents();
+ make_rotationToZ(P4projectile.getSpaceLikeComponents().getCoordinateSystem(),
+ P4projectile.getSpaceLikeComponents())} {
+ auto const pProjectile = P4projectile.getSpaceLikeComponents();
auto const pProjNormSquared = pProjectile.getSquaredNorm();
auto const pProjNorm = sqrt(pProjNormSquared);
- auto const eProjectile = Pprojectile.getTimeLikeComponent();
+ auto const eProjectile = P4projectile.getTimeLikeComponent();
auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
auto const s =
massTarget * massTarget + massProjectileSquared + 2 * eProjectile * massTarget;
@@ -43,7 +43,29 @@ namespace corsika {
coshEta, boost_.determinant() - 1);
}
- inline COMBoost::COMBoost(MomentumVector const& momentum, HEPEnergyType mass)
+ inline COMBoost::COMBoost(FourMomentum const& P4projectile,
+ FourMomentum const& P4target)
+ : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()}
+ , rotatedCS_{make_rotationToZ(
+ P4projectile.getSpaceLikeComponents().getCoordinateSystem(),
+ P4projectile.getSpaceLikeComponents() + P4target.getSpaceLikeComponents())} {
+ // this is the center-of-momentum CM frame
+ auto const pCM =
+ P4projectile.getSpaceLikeComponents() + P4target.getSpaceLikeComponents();
+ auto const pCM2 = pCM.getSquaredNorm();
+ auto const pCMnorm = sqrt(pCM2);
+
+ auto const s = (P4projectile + P4target).getNormSqr();
+ auto const sqrtS = sqrt(s);
+ auto const sinhEta = -pCMnorm / sqrtS;
+ auto const coshEta = sqrt(1 + pCM2 / s);
+
+ setBoost(coshEta, sinhEta);
+ CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+ coshEta, boost_.determinant() - 1);
+ }
+
+ inline COMBoost::COMBoost(MomentumVector const& momentum, HEPEnergyType const mass)
: originalCS_{momentum.getCoordinateSystem()}
, rotatedCS_{make_rotationToZ(momentum.getCoordinateSystem(), momentum)} {
auto const squaredNorm = momentum.getSquaredNorm();
@@ -56,12 +78,12 @@ namespace corsika {
}
template <typename FourVector>
- inline FourVector COMBoost::toCoM(FourVector const& p) const {
- auto pComponents = p.getSpaceLikeComponents().getComponents(rotatedCS_);
+ inline FourVector COMBoost::toCoM(FourVector const& p4) const {
+ auto pComponents = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
Eigen::Vector3d eVecRotated = pComponents.getEigenVector();
Eigen::Vector2d lab;
- lab << (p.getTimeLikeComponent() * (1 / 1_GeV)),
+ lab << (p4.getTimeLikeComponent() * (1 / 1_GeV)),
(eVecRotated(2) * (1 / 1_GeV).magnitude());
auto const boostedZ = boost_ * lab;
@@ -73,9 +95,9 @@ namespace corsika {
}
template <typename FourVector>
- inline FourVector COMBoost::fromCoM(FourVector const& p) const {
- auto pCM = p.getSpaceLikeComponents().getComponents(rotatedCS_);
- auto const Ecm = p.getTimeLikeComponent();
+ inline FourVector COMBoost::fromCoM(FourVector const& p4) const {
+ auto pCM = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
+ auto const Ecm = p4.getTimeLikeComponent();
Eigen::Vector2d com;
com << (Ecm * (1 / 1_GeV)), (pCM.getEigenVector()(2) * (1 / 1_GeV).magnitude());
@@ -83,7 +105,7 @@ namespace corsika {
CORSIKA_LOG_TRACE(
"COMBoost::fromCoM Ecm={} GeV"
" pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
- Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV, p.getNorm() / 1_GeV);
+ Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV, p4.getNorm() / 1_GeV);
auto const boostedZ = inverseBoost_ * com;
auto const E_lab = boostedZ(0) * 1_GeV;
@@ -102,7 +124,7 @@ namespace corsika {
return FourVector{E_lab, pLab};
}
- inline void COMBoost::setBoost(double coshEta, double sinhEta) {
+ inline void COMBoost::setBoost(double const coshEta, double const sinhEta) {
boost_ << coshEta, sinhEta, sinhEta, coshEta;
inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
}
diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
index 604ea7b9e..c0716be3a 100644
--- a/corsika/detail/modules/sibyll/InteractionModel.inl
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -60,9 +60,13 @@ namespace corsika::sibyll {
inline std::tuple<CrossSectionType, CrossSectionType>
InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
- HEPEnergyType const sqrtSnn) const {
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
- isValid(projectileId, targetId, sqrtSnn); // throws
+ int targetSibCode = 1; // nucleon or particle count
+ if (is_nucleus(targetId)) { targetSibCode = get_nucleus_A(targetId); }
+ // sqrtS per target nucleon
+ HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
double dummy, dum1, dum3, dum4, dumdif[3]; // dummies needed for fortran call
int const iBeam = corsika::sibyll::getSibyllXSCode(
@@ -92,24 +96,23 @@ namespace corsika::sibyll {
template <typename TSecondaryView>
inline void InteractionModel::doInteraction(TSecondaryView& secondaries,
- COMBoost const& boost,
Code const projectileId,
Code const targetId,
- HEPEnergyType const sqrtSnn) {
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) {
+
+ int targetSibCode = 1; // nucleon or particle count
+ if (is_nucleus(targetId)) { targetSibCode = get_nucleus_A(targetId); }
+ CORSIKA_LOG_DEBUG("sibyll code: {} (nucleon/particle count)", targetSibCode);
+
+ // sqrtS per target nucleon
+ HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
+ COMBoost const boost(projectileP4, targetP4 / targetSibCode);
isValid(projectileId, targetId, sqrtSnn); // throws
CORSIKA_LOG_DEBUG("pId={} tId={} sqrtSnn={}GeV", projectileId, targetId, sqrtSnn);
- int targetSibCode = -1;
- if (is_nucleus(targetId)) {
- targetSibCode = get_nucleus_A(targetId);
- } else {
- // nucleon target: p or n
- targetSibCode = 1;
- }
- CORSIKA_LOG_DEBUG("sibyll code: {}", targetSibCode);
-
// beam id for sibyll
int const projectileSibyllCode = corsika::sibyll::convertToSibyllRaw(projectileId);
@@ -177,6 +180,6 @@ namespace corsika::sibyll {
Elab_final / 1_GeV, (Elab_final - Elab_initial) / Elab_initial * 100,
constants::nucleonMass * get_nwounded() / 1_GeV,
(Plab_final / 1_GeV).getComponents());
- }
+ } // namespace corsika::sibyll
} // namespace corsika::sibyll
diff --git a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
index 7c84c79e5..906cfa8b6 100644
--- a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
@@ -116,8 +116,7 @@ namespace corsika::sibyll {
if (ptarg == Code::Argon) continue; // NEED TO IGNORE Argon ....
++k;
CORSIKA_LOG_DEBUG("init target component: {}", ptarg);
- int const ib =
- get_nucleus_A(ptarg); // this assumes the universe is only made out of "Nuclei"
+ int const ib = get_nucleus_A(ptarg);
hadronicInteraction_.isValid(Code::Proton, ptarg, 100_GeV); // throws
targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
// loop over energies, fNEnBins log. energy bins
@@ -125,10 +124,15 @@ namespace corsika::sibyll {
// hard coded energy grid, has to be aligned to definition in signuc2!!, no
// comment..
HEPEnergyType const Ecm = pow(10., 1. + 1. * i) * 1_GeV;
+ 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});
// get p-p cross sections
auto const protonId = Code::Proton;
- auto const [siginel, sigela] =
- hadronicInteraction_.getCrossSectionInelEla(protonId, protonId, Ecm);
+ auto const [siginel, sigela] = hadronicInteraction_.getCrossSectionInelEla(
+ protonId, protonId, projectileP4, targetP4);
const double dsig = siginel / 1_mb;
const double dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
@@ -170,8 +174,10 @@ namespace corsika::sibyll {
CrossSectionType inline NuclearInteractionModel<
TEnvironment, TNucleonModel>::getCrossSection(Code const projectileId,
Code const targetId,
- HEPEnergyType const sqrtSnn) const {
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+ HEPEnergyType const sqrtSnn = (projectileP4 + targetP4).getNorm();
isValid(projectileId, targetId, sqrtSnn); // throws
HEPEnergyType const LabEnergyPerNuc =
static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
@@ -184,26 +190,29 @@ namespace corsika::sibyll {
template <typename TEnvironment, typename TNucleonModel>
template <typename TSecondaryView>
inline void NuclearInteractionModel<TEnvironment, TNucleonModel>::doInteraction(
- TSecondaryView& view, COMBoost const& boost, Code const projectileId,
- Code const targetId, HEPEnergyType const sqrtSnn) {
+ TSecondaryView& view, Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4) {
- isValid(projectileId, targetId, sqrtSnn); // throws
+ // model is only designed for projectile nuclei. Collisions are broken down into
+ // "nucleon-target" collisions.
+ size_t const projectileA = get_nucleus_A(projectileId);
+
+ // 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
// projectile is always nucleus!
- unsigned int const projectileA = get_nucleus_A(projectileId);
+ // Elab corresponding to sqrtSnucleon -> fixed target projectile
+ COMBoost const boost(nucleonP4, targetP4);
- CORSIKA_LOG_DEBUG("pId={} tId={} sqrtSnn={}GeV Aproj={}", projectileId, targetId,
- sqrtSnn / 1_GeV, projectileA);
+ CORSIKA_LOG_DEBUG("pId={} tId={} sqrtSnucleon={}GeV Aproj={}", projectileId, targetId,
+ sqrtSnucleon / 1_GeV, projectileA);
count_++;
- HEPEnergyType const ProjMass = get_mass(projectileId);
-
- // lab. Energy per projectile nucleon
- HEPEnergyType const eProjectileLab =
- static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
- HEPMomentumType const pProjectileLab =
- sqrt(static_pow<2>(eProjectileLab) - static_pow<2>(constants::nucleonMass));
- MomentumVector const p3ProjectileLab(boost.getRotatedCS(),
- {0_GeV, 0_GeV, pProjectileLab});
+ // lab. momentum per projectile nucleon
+ HEPMomentumType const pNucleonLab = nucleonP4.getSpaceLikeComponents().getNorm();
+ // nucleon momentum in direction of CM motion (lab system)
+ MomentumVector const p3NucleonLab(boost.getRotatedCS(), {0_GeV, 0_GeV, pNucleonLab});
/*
FOR NOW: allow nuclei with A<18 or protons/nucleon only.
@@ -211,8 +220,10 @@ namespace corsika::sibyll {
allowed air in atmosphere also contains some Argon.
*/
int kATarget = -1;
+ size_t targetA = 1;
if (is_nucleus(targetId)) {
kATarget = get_nucleus_A(targetId);
+ targetA = kATarget;
} else if (targetId == Code::Proton || targetId == Code::Neutron ||
targetId == Code::Hydrogen) {
kATarget = 1;
@@ -227,7 +238,9 @@ namespace corsika::sibyll {
// (needed to determine number of nucleon-nucleon scatterings)
auto const protonId = Code::Proton;
auto const [prodCrossSection, elaCrossSection] =
- hadronicInteraction_.getCrossSectionInelEla(protonId, protonId, sqrtSnn);
+ hadronicInteraction_.getCrossSectionInelEla(
+ protonId, protonId, nucleonP4,
+ targetP4 / targetA); // todo check, wrong RU
double const sigProd = prodCrossSection / 1_mb;
double const sigEla = elaCrossSection / 1_mb;
// sample number of interactions (only input variables, output in common cnucms)
@@ -270,8 +283,10 @@ namespace corsika::sibyll {
// (LCOV_EXCL_STOP)
// position and time of interaction, not used in NUCLIB
- Point pOrig{boost.getOriginalCS(), {0_m, 0_m, 0_m}};
- TimeType delay = 0_s; // there is no time in sibyll
+ auto const& projectile = view.parent();
+ // position and time of interaction, not used in NUCLI
+ Point const& pOrig = projectile.getPosition();
+ TimeType const delay = projectile.getTime();
CORSIKA_LOG_DEBUG("Interaction: position of interaction: {} {}",
pOrig.getCoordinates(), delay / 1_s);
@@ -298,12 +313,8 @@ namespace corsika::sibyll {
// CORSIKA 7 way
// spectators inherit momentum from original projectile
- double const mass_ratio = mass / ProjMass;
- auto const p3lab = p3ProjectileLab * mass_ratio;
-
- CORSIKA_LOG_DEBUG("mass ratio {}, fragment momentum {}", mass_ratio,
- p3lab.getComponents() / 1_GeV);
-
+ auto const p3lab = p3NucleonLab * nuclA;
+ CORSIKA_LOG_DEBUG("fragment momentum {}", p3lab.getComponents() / 1_GeV);
view.addSecondary(std::make_tuple(specCode, p3lab, pOrig, delay));
}
@@ -318,8 +329,7 @@ namespace corsika::sibyll {
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
- double const mass_ratio = get_mass(elaNucCode) / ProjMass;
- auto const p3lab = p3ProjectileLab * mass_ratio;
+ auto const p3lab = p3NucleonLab;
view.addSecondary(std::make_tuple(elaNucCode, p3lab, pOrig, delay));
}
@@ -327,20 +337,20 @@ namespace corsika::sibyll {
CORSIKA_LOG_DEBUG("calculate inelastic nucleon-nucleon interactions..");
for (int j = 0; j < nInelNucleons; ++j) {
// TODO: sample neutron or proton
- auto pCode = Code::Proton;
+ auto const pCode = Code::Proton;
// temporarily add to stack, will be removed after interaction in DoInteraction
CORSIKA_LOG_DEBUG("inelastic interaction no. {}", j);
typename TSecondaryView::inner_stack_value_type nucleonStack;
auto inelasticNucleon =
- nucleonStack.addParticle(std::make_tuple(pCode, p3ProjectileLab, pOrig, delay));
+ nucleonStack.addParticle(std::make_tuple(pCode, p3NucleonLab, pOrig, delay));
inelasticNucleon.setNode(view.getProjectile().getNode());
// create inelastic interaction for each nucleon
CORSIKA_LOG_TRACE("calling HadronicInteraction...");
// create new StackView for each of the nucleons
TSecondaryView nucleon_secondaries(inelasticNucleon);
// all inner hadronic event generator
- hadronicInteraction_.doInteraction(nucleon_secondaries, boost, pCode, targetId,
- sqrtSnn);
+ hadronicInteraction_.doInteraction(nucleon_secondaries, pCode, targetId, nucleonP4,
+ targetP4);
for (const auto& pSec : nucleon_secondaries) {
view.addSecondary(std::make_tuple(pSec.getPID(), pSec.getMomentum(),
pSec.getPosition(), pSec.getTime()));
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index 796655a03..ac7a899eb 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -12,11 +12,12 @@
#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/random/ExponentialDistribution.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/random/UniformRealDistribution.hpp>
#include <corsika/framework/stack/SecondaryView.hpp>
-#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/media/Environment.hpp>
@@ -33,8 +34,6 @@
namespace corsika {
- class COMBoost; // fwd-decl
-
/**
*
* The Cascade class is constructed from template arguments making
@@ -98,7 +97,7 @@ namespace corsika {
/**
* set the nodes for all particles on the stack according to their numerical
- * position
+ * position.
*/
void setNodes();
@@ -129,8 +128,7 @@ namespace corsika {
void step(particle_type& vParticle);
ProcessReturn decay(stack_view_type& view, InverseTimeType initial_inv_decay_time);
- ProcessReturn interaction(stack_view_type& view, COMBoost const& boost,
- HEPEnergyType const sqrtSnn,
+ ProcessReturn interaction(stack_view_type& view, FourMomentum const& projectileP4,
NuclearComposition const& composition,
CrossSectionType const initial_cross_section);
void setEventType(stack_view_type& view, history::EventType);
diff --git a/corsika/framework/geometry/FourVector.hpp b/corsika/framework/geometry/FourVector.hpp
index 0b87accec..0a23f4e94 100644
--- a/corsika/framework/geometry/FourVector.hpp
+++ b/corsika/framework/geometry/FourVector.hpp
@@ -9,36 +9,42 @@
#pragma once
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
#include <type_traits>
+/**
+ * @file FourVector.hpp
+ * @author Ralf Ulrich
+ * @brief General FourVector object.
+ * @date 2021-10-16
+ */
+
namespace corsika {
/**
- Description of physical four-vectors
-
- FourVector fully supports units, e.g. E in [GeV/c] and p in [GeV],
- or also t in [s] and r in [m], etc.
-
- However, for HEP applications it is also possible to use E and p
- both in [GeV].
-
- Thus, the input units of time-like and space-like coordinates
- must either be idential (e.g. GeV) or scaled by "c" as in
- [E/c]=[p].
-
-
- The FourVector can return its squared-norm \ref getNormSqr and its
- norm \ref getNorm, whereas norm is sqrt(abs(norm-squared)). The
- physical units are always calculated and returned properly.
-
- FourVector can also return if it is TimeLike, SpaceLike or PhotonLike.
-
- When a FourVector is initialized with a lvalue references,
- e.g. as `FourVector<TimeType&, Vector<length_d>&>`, references
- are also used as internal data types, which should lead to
- complete disappearance of the FourVector class during
- optimization.
+ * Description of physical four-vectors
+ *
+ * FourVector fully supports units, e.g. E in [GeV/c] and p in [GeV],
+ * or also t in [s] and r in [m], etc.
+ *
+ * However, for HEP applications it is also possible to use E and p
+ * both in [GeV].
+ *
+ * Thus, the input units of time-like and space-like coordinates
+ * must either be idential (e.g. GeV) or scaled by "c" as in
+ * [E/c]=[p].
+ *
+ * The FourVector can return its squared-norm \ref getNormSqr and its
+ * norm \ref getNorm, whereas norm is sqrt(abs(norm-squared)). The
+ * physical units are always calculated and returned properly.
+ *
+ * FourVector can also return if it is TimeLike, SpaceLike or PhotonLike.
+ *
+ * When a FourVector is initialized with a lvalue references,
+ * e.g. as `FourVector<TimeType&, Vector<length_d>&>`, references
+ * are also used as internal data types, which should lead to
+ * complete disappearance of the FourVector class during
+ * optimization.
*/
template <typename TTimeType, typename TSpaceVecType>
@@ -73,13 +79,11 @@ namespace corsika {
, spaceLike_(eS) {}
/**
- *
* @return timeLike_
*/
TTimeType getTimeLikeComponent() const;
/**
- *
* @return spaceLike_
*/
TSpaceVecType& getSpaceLikeComponents();
@@ -124,12 +128,12 @@ namespace corsika {
FourVector& operator/(double const);
/**
- Scalar product of two FourVectors
-
- Note that the product between two 4-vectors assumes that you use
- the same "c" convention for both. Only the LHS vector is checked
- for this. You cannot mix different conventions due to
- unit-checking.
+ * Scalar product of two FourVectors.
+ *
+ * Note that the product between two 4-vectors assumes that you use
+ * the same "c" convention for both. Only the LHS vector is checked
+ * for this. You cannot mix different conventions due to
+ * unit-checking.
*/
norm_type operator*(FourVector const& b);
@@ -152,7 +156,7 @@ namespace corsika {
* value-copies.
* @{
*
- **/
+ */
friend FourVector<time_type, space_vec_type> operator+(FourVector const& a,
FourVector const& b) {
return FourVector<time_type, space_vec_type>(a.timeLike_ + b.timeLike_,
@@ -179,20 +183,25 @@ namespace corsika {
private:
/**
- This function is there to automatically remove the eventual
- extra factor of "c" for the time-like quantity.
+ * This function is there to automatically remove the eventual
+ * extra factor of "c" for the time-like quantity.
*/
norm_square_type getTimeSquared() const;
};
/**
* streaming operator
- **/
+ */
template <typename TTimeType, typename TSpaceVecType>
std::ostream& operator<<(std::ostream& os,
corsika::FourVector<TTimeType, TSpaceVecType> const& qv);
+ /**
+ * @typedef FourMomentum A FourVector with HEPEnergyType and MomentumVector.
+ */
+ typedef FourVector<HEPEnergyType, MomentumVector> FourMomentum;
+
} // namespace corsika
#include <corsika/detail/framework/geometry/FourVector.inl>
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index b14c9d8d3..4ecf3d0cd 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -29,6 +29,8 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+
namespace corsika {
class COMBoost; // fwd-decl
@@ -242,16 +244,43 @@ namespace corsika {
ContinuousProcessStepLength getMaxStepLength(TParticle&& particle, TTrack&& vTrack);
CrossSectionType getCrossSection(Code const projectileId, Code const targetId,
- HEPEnergyType const sqrtSnn) const;
+ FourMomentum const& projectilP4,
+ FourMomentum const& targetP4) const;
template <typename TParticle>
TimeType getLifetime(TParticle&& particle) {
return 1. / getInverseLifetime(particle);
}
+ /**
+ * @brief Selects one concrete InteractionProcess and samples a target nucleus from
+ * the material.
+ *
+ * The selectInteraction method statically loops over all active InteractionProcess
+ * and calculates the material-weighted cross section for all of them. In an iterative
+ * way those cross sections are summed up. The random number cx_select, uniformely
+ * drawn from the cross section before energy losses, is used to discriminate the
+ * selected sub-process here. If the cross section after the step smaller than it was
+ * before, there is a non-zero probability that the particle survives and no
+ * interaction takes place. This method becomes imprecise when cross section rise with
+ * falling energies.
+ *
+ * If a sub-process was selected, the target nucleus is selected from the material
+ * (weighted with cross section). The interaction is then executed.
+ *
+ * @tparam TSecondaryView Object type as storage for new secondary particles.
+ * @tparam TRNG Object type to produce random numbers.
+ * @param view Object to store new secondary particles.
+ * @param projectileP4 The four momentum of the projectile.
+ * @param composition The environment/material composition.
+ * @param rng Random number object.
+ * @param cx_select Drawn random numer, uniform between [0, cx_initial]
+ * @param cx_sum For interal use, to sum up cross section contributions.
+ * @return ProcessReturn
+ */
template <typename TSecondaryView, typename TRNG>
inline ProcessReturn selectInteraction(
- TSecondaryView&& view, COMBoost const& boost, HEPEnergyType const sqrtSnn,
+ TSecondaryView&& view, FourMomentum const& projectileP4,
NuclearComposition const& composition, TRNG&& rng,
CrossSectionType const cx_select,
CrossSectionType cx_sum = CrossSectionType::zero());
diff --git a/corsika/framework/process/SwitchProcessSequence.hpp b/corsika/framework/process/SwitchProcessSequence.hpp
index d3f01120c..74f07fff1 100644
--- a/corsika/framework/process/SwitchProcessSequence.hpp
+++ b/corsika/framework/process/SwitchProcessSequence.hpp
@@ -150,8 +150,8 @@ namespace corsika {
HEPEnergyType const sqrtSnn) const;
template <typename TSecondaryView, typename TRNG>
- ProcessReturn selectInteraction(TSecondaryView& view, COMBoost const& boost,
- HEPEnergyType const sqrtSnn,
+ ProcessReturn selectInteraction(TSecondaryView& view,
+ FourMomentum const& projectileP4,
NuclearComposition const& composition, TRNG& rng,
CrossSectionType const cx_select,
CrossSectionType cx_sum = CrossSectionType::zero());
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index d1f9ee785..8bba6b43d 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -19,29 +19,29 @@
namespace corsika {
/**
- @defgroup Utilities
-
- Collection of classes and methods to perform recurring tasks.
- **/
+ * @defgroup Utilities
+ *
+ * Collection of classes and methods to perform recurring tasks.
+ */
/**
- @class COMBoost
- @ingroup Utilities
-
- This utility class handles Lorentz boost between different
- referenence frames, using FourVector.
-
- The class is initialized with projectile and optionally target
- energy/momentum data. During initialization, a rotation matrix is
- calculated to represent the projectile movement (and thus the
- boost) along the z-axis. Also the inverse of this rotation is
- calculated. The Lorentz boost matrix and its inverse are
- determined as 2x2 matrices considering the energy and
- pz-momentum.
-
- Different constructors are offered with different specialization
- for the cases of collisions (projectile-target) or just decays
- (projectile only).
+ * @class COMBoost
+ * @ingroup Utilities
+ *
+ * This utility class handles Lorentz boost between different
+ * referenence frames, using FourVector.
+ *
+ * The class is initialized with projectile and optionally target
+ * energy/momentum data. During initialization, a rotation matrix is
+ * calculated to represent the projectile movement (and thus the
+ * boost) along the z-axis. Also the inverse of this rotation is
+ * calculated. The Lorentz boost matrix and its inverse are
+ * determined as 2x2 matrices considering the energy and
+ * pz-momentum.
+ *
+ * Different constructors are offered with different specialization
+ * for the cases of collisions (projectile-target) or just decays
+ * (projectile only).
*/
class COMBoost {
@@ -49,19 +49,21 @@ namespace corsika {
public:
//! construct a COMBoost given four-vector of projectile and mass of target (target at
//! rest)
- COMBoost(FourVector<HEPEnergyType, MomentumVector> const& Pprojectile,
- HEPEnergyType const massTarget);
+ COMBoost(FourMomentum const& P4projectile, HEPEnergyType const massTarget);
+
+ //! construct a COMBoost given two four-vectors of projectile target
+ COMBoost(FourMomentum const& P4projectile, FourMomentum const& P4target);
//! construct a COMBoost to boost into the rest frame given a 3-momentum and mass
- COMBoost(MomentumVector const& momentum, HEPEnergyType mass);
+ COMBoost(MomentumVector const& momentum, HEPEnergyType const mass);
//! transforms a 4-momentum from lab frame to the center-of-mass frame
template <typename FourVector>
- FourVector toCoM(FourVector const& p) const;
+ FourVector toCoM(FourVector const& p4) const;
//! transforms a 4-momentum from the center-of-mass frame back to lab frame
template <typename FourVector>
- FourVector fromCoM(FourVector const& p) const;
+ FourVector fromCoM(FourVector const& p4) const;
//! returns the rotated coordinate system
CoordinateSystemPtr getRotatedCS() const;
@@ -71,7 +73,7 @@ namespace corsika {
protected:
//! internal method
- void setBoost(double coshEta, double sinhEta);
+ void setBoost(double const coshEta, double const sinhEta);
private:
Eigen::Matrix2d boost_;
diff --git a/corsika/modules/sibyll/InteractionModel.hpp b/corsika/modules/sibyll/InteractionModel.hpp
index bf4b13081..2b4191368 100644
--- a/corsika/modules/sibyll/InteractionModel.hpp
+++ b/corsika/modules/sibyll/InteractionModel.hpp
@@ -64,7 +64,8 @@ namespace corsika::sibyll {
* @return a tuple of: inelastic cross section, elastic cross section
*/
std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
- Code const projectile, Code const target, HEPEnergyType const sqrtSnn) const;
+ Code const projectile, Code const target, FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const;
/**
* @brief returns inelastic (production) cross section.
@@ -82,8 +83,10 @@ namespace corsika::sibyll {
* elastic cross section
*/
CrossSectionType getCrossSection(Code const projectile, Code const target,
- HEPEnergyType const sqrtSnn) const {
- return std::get<0>(getCrossSectionInelEla(projectile, target, sqrtSnn));
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+ return std::get<0>(
+ getCrossSectionInelEla(projectile, target, projectileP4, targetP4));
}
/**
@@ -92,8 +95,8 @@ namespace corsika::sibyll {
*/
template <typename TSecondaries>
- void doInteraction(TSecondaries&, COMBoost const& boost, Code const projectile,
- Code const target, HEPEnergyType const sqrtSnn);
+ void doInteraction(TSecondaries&, Code const projectile, Code const target,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
private:
HEPEnergyType constexpr getMinEnergyCoM() const { return minEnergyCoM_; }
diff --git a/corsika/modules/sibyll/NuclearInteractionModel.hpp b/corsika/modules/sibyll/NuclearInteractionModel.hpp
index 4606f80e2..807214574 100644
--- a/corsika/modules/sibyll/NuclearInteractionModel.hpp
+++ b/corsika/modules/sibyll/NuclearInteractionModel.hpp
@@ -44,11 +44,13 @@ namespace corsika::sibyll {
unsigned int constexpr getMaxNFragments() const { return gMaxNFragments_; }
unsigned int constexpr getNEnergyBins() const { return gNEnBins_; }
- CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const) const;
+ CrossSectionType getCrossSection(Code const, Code const,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const;
template <typename TSecondaryView>
- void doInteraction(TSecondaryView&, COMBoost const& boost, Code const, Code const,
- HEPEnergyType const);
+ void doInteraction(TSecondaryView&, Code const, Code const,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
private:
int count_ = 0;
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index d2d9ef068..f2ac504e3 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -17,11 +17,10 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-
#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
@@ -89,12 +88,13 @@ public:
class ProcessSplit : public InteractionProcess<ProcessSplit> {
public:
- CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const) const {
+ CrossSectionType getCrossSection(Code const, Code const, FourMomentum const&,
+ FourMomentum const&) const {
return 1_mb;
}
template <typename TView>
- void doInteraction(TView& view, COMBoost const&, Code, Code, HEPEnergyType) {
+ void doInteraction(TView& view, Code, Code, FourMomentum const&, FourMomentum const&) {
++calls_;
auto vP = view.getProjectile();
const HEPEnergyType Ekin = vP.getKineticEnergy();
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index 312e8495a..ac4c6392d 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -58,7 +58,10 @@ struct DummyData {
double data_[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
typedef DummyNode node_type; // for BoundaryCrossingProcess
Code getPID() const { return Code::Proton; }
- // MomentumVector getMomentum() const {}
+ MomentumVector getMomentum() const {
+ // only need the coordinate system
+ return MomentumVector{get_root_CoordinateSystem(), 0_eV, 0_eV, 0_eV};
+ }
HEPEnergyType getEnergy() const { return 10_GeV; }
unsigned int getNuclearA() const { return 1; }
};
@@ -208,13 +211,14 @@ public:
}
template <typename TView>
- void doInteraction(TView& v, COMBoost const&, Code const, Code const,
- HEPEnergyType const) const {
+ void doInteraction(TView& v, Code const, Code const, FourMomentum const&,
+ FourMomentum const&) const {
checkInteract |= 1;
for (int i = 0; i < nData; ++i) v.parent().data_[i] += 1 + i;
}
- CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const) const {
+ CrossSectionType getCrossSection(Code const, Code const, FourMomentum const&,
+ FourMomentum const&) const {
return 10_mb;
}
@@ -234,14 +238,15 @@ public:
}
template <typename TView>
- void doInteraction(TView& v, COMBoost const&, Code const, Code const,
- HEPEnergyType const) const {
+ void doInteraction(TView& v, Code const, Code const, FourMomentum const&,
+ FourMomentum const&) const {
checkInteract |= 2;
for (int i = 0; i < nData; ++i) v.parent().data_[i] /= 1.1;
CORSIKA_LOG_DEBUG("Process2::doInteraction");
}
- CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const) const {
+ CrossSectionType getCrossSection(Code const, Code const, FourMomentum const&,
+ FourMomentum const&) const {
CORSIKA_LOG_DEBUG("Process2::getCrossSection");
return 20_mb;
}
@@ -262,14 +267,15 @@ public:
}
template <typename TView>
- void doInteraction(TView& v, COMBoost const&, Code const, Code const,
- HEPEnergyType const) const {
+ void doInteraction(TView& v, Code const, Code const, FourMomentum const&,
+ FourMomentum const&) const {
checkInteract |= 4;
for (int i = 0; i < nData; ++i) v.parent().data_[i] *= 1.01;
CORSIKA_LOG_DEBUG("Process3::doInteraction");
}
- CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const) const {
+ CrossSectionType getCrossSection(Code const, Code const, FourMomentum const&,
+ FourMomentum const&) const {
CORSIKA_LOG_DEBUG("Process3::getCrossSection");
return 30_mb;
}
@@ -297,8 +303,8 @@ public:
return ProcessReturn::Ok;
}
template <typename TView>
- void doInteraction(TView&, COMBoost const&, Code const, Code const,
- HEPEnergyType const) const {
+ void doInteraction(TView&, Code const, Code const, FourMomentum const&,
+ FourMomentum const&) const {
checkInteract |= 8;
}
@@ -613,7 +619,7 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
- logging::set_level(logging::level::trace);
+ logging::set_level(logging::level::info);
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
@@ -717,9 +723,9 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
particle.data_[0] = 100; // data positive --> sequence1
DummyRNG rng;
- COMBoost const noBoost({10_GeV, {rootCS, {0_eV, 0_eV, 0_eV}}}, 0_GeV);
+ FourMomentum const projectileP4{10_GeV, {rootCS, {0_eV, 0_eV, 0_eV}}};
NuclearComposition const noComposition({Code::Nitrogen}, {1});
- sequence3.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
sequence3.selectDecay(view, time_select);
CHECK(checkInteract == 0b100); // this is Process3
CHECK(checkDecay == 0b001); // this is Decay1
@@ -727,7 +733,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
CHECK(checkSec == 0);
cx_select = 1.01 * 30_mb;
checkInteract = 0;
- sequence3.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
CHECK(checkInteract == 0b001); // this is Process1
checkDecay = 0;
@@ -735,7 +741,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkSec = 0;
checkCont = 0;
particle.data_[0] = -100; // data negative --> sequence2
- sequence3.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
sequence3.selectDecay(view, time_select);
CHECK(checkInteract == 0b010); // this is Process2
CHECK(checkDecay == 0b010); // this is Decay2
@@ -763,7 +769,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkSec = 0;
checkCont = 0;
particle.data_[0] = -100; // data negative --> sequence1
- sequence4.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
+ sequence4.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
sequence4.doSecondaries(view);
sequence4.selectDecay(view, time_select);
sequence4.doSecondaries(view);
@@ -777,7 +783,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
time_select = 1e5 / second;
checkDecay = 0;
checkInteract = 0;
- sequence3.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
+ sequence3.selectInteraction(view, projectileP4, noComposition, rng, cx_select);
sequence3.selectDecay(view, time_select);
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 73e812b5c..a08e33d8c 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -161,12 +161,14 @@ TEST_CASE("SibyllInterface", "modules") {
CHECK_NOTHROW(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});
+ FourMomentum const bP4(1_GeV, {cs, 0_GeV, 0_GeV, 0_GeV});
auto const [xs_prod_pp, xs_ela_pp] =
- model.getCrossSectionInelEla(Code::Proton, Code::Proton, 100_GeV);
+ model.getCrossSectionInelEla(Code::Proton, Code::Proton, aP4, bP4);
auto const [xs_prod_pn, xs_ela_pn] =
- model.getCrossSectionInelEla(Code::Proton, Code::Neutron, 100_GeV);
+ model.getCrossSectionInelEla(Code::Proton, Code::Neutron, aP4, bP4);
auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] =
- model.getCrossSectionInelEla(Code::Proton, Code::Hydrogen, 100_GeV);
+ model.getCrossSectionInelEla(Code::Proton, Code::Hydrogen, aP4, bP4);
CHECK(xs_prod_pp == xs_prod_pHydrogen);
CHECK(xs_prod_pp == xs_prod_pn);
CHECK(xs_ela_pp == xs_ela_pHydrogen);
@@ -178,16 +180,15 @@ TEST_CASE("SibyllInterface", "modules") {
SECTION("InteractionInterface - low energy") {
const HEPEnergyType P0 = 60_GeV;
- MomentumVector plab =
- MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
+ MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
// also print particles after sibyll was called
corsika::sibyll::InteractionModel model;
model.setVerbose(true);
HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(Proton::mass));
- HEPEnergyType const sqrtSnn = sqrt(2 * Elab * constants::nucleonMass);
+ FourMomentum const projectileP4(Elab, plab);
+ FourMomentum const nucleonP4(Oxygen::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
view.clear();
- COMBoost boost = getCOMboost(Elab, plab, cs);
- model.doInteraction(view, boost, Code::Proton, Code::Oxygen, sqrtSnn);
+ model.doInteraction(view, Code::Proton, Code::Oxygen, projectileP4, nucleonP4);
auto const pSum = sumMomentum(view, cs);
/*
@@ -254,7 +255,7 @@ TEST_CASE("SibyllInterface", "modules") {
Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] CrossSectionType const cx =
- model.getCrossSection(Code::Proton, Code::Oxygen, sqrtSnn);
+ model.getCrossSection(Code::Proton, Code::Oxygen, projectileP4, nucleonP4);
CHECK(cx / 1_mb == Approx(300).margin(1));
// CHECK(view.getEntries() == 9); //! \todo: this was 20 before refactory-2020: check
// "also sibyll not stable wrt. to compiler
@@ -264,24 +265,22 @@ TEST_CASE("SibyllInterface", "modules") {
SECTION("NuclearInteractionInterface") {
HEPMomentumType const P0 = 50_TeV;
- MomentumVector plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
+ MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
corsika::sibyll::InteractionModel hmodel;
NuclearInteractionModel model(hmodel, *env);
Code const pid = Code::Oxygen;
- size_t const A = get_nucleus_A(pid);
HEPEnergyType const Elab = sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid)));
- HEPEnergyType const sqrtSnn = sqrt((Elab / A + constants::nucleonMass + P0 / A) *
- (Elab / A + constants::nucleonMass - P0 / A));
- model.doInteraction(view, getCOMboost(Elab, plab, cs), Code::Oxygen, Code::Oxygen,
- sqrtSnn);
- CrossSectionType const cx =
- model.getCrossSection(Code::Oxygen, Code::Oxygen, sqrtSnn);
+ FourMomentum const P4(Elab, plab);
+ FourMomentum const targetP4(get_mass(Code::Oxygen),
+ MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
+ model.doInteraction(view, pid, Code::Oxygen, P4, targetP4);
+ CrossSectionType const cx = model.getCrossSection(pid, Code::Oxygen, P4, targetP4);
// Felix, are those changes OK? Below are the checks before refactory-2020
// CHECK(length / 1_g * 1_cm * 1_cm == Approx(44.2).margin(.1));
// CHECK(view.getSize() == 11);
- CHECK(cx / 1_mb == Approx(1300).margin(300)); // this is not physics validation
+ CHECK(cx / 1_mb == Approx(1100).margin(100)); // this is not physics validation
// CHECK(view.getSize() == 20); // also sibyll not stable wrt. to compiler changes
- CHECK(view.getSize() == Approx(100).margin(90)); // this is not physics validation
+ CHECK(view.getSize() == Approx(200).margin(90)); // this is not physics validation
}
}
--
GitLab