From c9c1540f1abf818e77f50a34ef734f38e11c79f9 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sat, 9 Oct 2021 09:55:15 +0200
Subject: [PATCH] large overhaul, hadint interface
---
corsika/detail/framework/core/Cascade.inl | 176 +++---
.../framework/process/InteractionProcess.hpp | 58 +-
.../framework/process/ProcessSequence.inl | 170 +++--
.../process/SwitchProcessSequence.inl | 179 ++++--
corsika/detail/framework/utility/COMBoost.inl | 13 +-
corsika/detail/media/NuclearComposition.inl | 79 ++-
corsika/detail/modules/sibyll/Interaction.inl | 351 -----------
.../modules/sibyll/InteractionModel.inl | 183 ++++++
.../modules/sibyll/NuclearInteraction.inl | 586 ------------------
.../sibyll/NuclearInteractionModel.inl | 363 +++++++++++
corsika/detail/stack/VectorStack.inl | 9 +-
corsika/framework/core/Cascade.hpp | 8 +-
corsika/framework/process/BaseProcess.hpp | 4 +-
corsika/framework/process/DecayProcess.hpp | 4 +-
.../framework/process/InteractionProcess.hpp | 83 ++-
corsika/framework/process/ProcessSequence.hpp | 356 +++++------
.../process/SwitchProcessSequence.hpp | 21 +-
corsika/framework/utility/COMBoost.hpp | 3 +
corsika/media/NuclearComposition.hpp | 69 ++-
corsika/modules/Sibyll.hpp | 39 +-
corsika/modules/sibyll/Interaction.hpp | 73 ---
corsika/modules/sibyll/NuclearInteraction.hpp | 70 ---
corsika/stack/VectorStack.hpp | 8 +-
tests/common/SetupTestEnvironment.hpp | 2 +-
tests/framework/testCascade.cpp | 32 +-
tests/framework/testCascade.hpp | 5 +-
tests/framework/testProcessSequence.cpp | 85 ++-
tests/media/CMakeLists.txt | 1 +
tests/media/testEnvironment.cpp | 7 +-
tests/media/testMagneticField.cpp | 3 +-
tests/media/testMedium.cpp | 3 +-
tests/media/testNuclearComposition.cpp | 68 ++
tests/media/testRefractiveIndex.cpp | 6 +-
tests/media/testShowerAxis.cpp | 3 +-
tests/modules/testCONEX.cpp | 2 +-
tests/modules/testSibyll.cpp | 170 +++--
tests/modules/testTracking.cpp | 12 +-
37 files changed, 1575 insertions(+), 1729 deletions(-)
delete mode 100644 corsika/detail/modules/sibyll/Interaction.inl
create mode 100644 corsika/detail/modules/sibyll/InteractionModel.inl
delete mode 100644 corsika/detail/modules/sibyll/NuclearInteraction.inl
create mode 100644 corsika/detail/modules/sibyll/NuclearInteractionModel.inl
delete mode 100644 corsika/modules/sibyll/Interaction.hpp
delete mode 100644 corsika/modules/sibyll/NuclearInteraction.hpp
create mode 100644 tests/media/testNuclearComposition.cpp
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 64ca0a3aa..383d22265 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -9,14 +9,21 @@
#pragma once
#include <corsika/framework/core/PhysicalUnits.hpp>
+
#include <corsika/framework/process/ProcessReturn.hpp>
#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
#include <corsika/framework/process/ContinuousProcessIndex.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/utility/COMBoost.hpp>
+
#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
#include <cassert>
#include <cmath>
@@ -59,41 +66,82 @@ namespace corsika {
inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceInteraction() {
CORSIKA_LOG_TRACE("forced interaction!");
setNodes();
- auto vParticle = stack_.getNextParticle();
- stack_view_type secondaries(vParticle);
- interaction(secondaries, sequence_.getInverseInteractionLength(vParticle));
+ auto particle = stack_.getNextParticle();
+ stack_view_type secondaries(particle);
+
+ auto const* currentLogicalNode = particle.getNode();
+ // assert that particle stays outside void Universe if it has no
+ // model properties set
+ assert((currentLogicalNode != &*environment_.getUniverse() ||
+ environment_.getUniverse()->hasModelProperties()) &&
+ "FATAL: The environment model has no valid properties set!");
+ NuclearComposition const& composition =
+ currentLogicalNode->getModelProperties().getNuclearComposition();
+
+ // determine sqrtS per nucleon pair, sqrtS_NN
+ Code const projectileId = particle.getPID();
+ unsigned int const projectileA =
+ (is_nucleus(projectileId) ? particle.getNuclearA() : 1);
+ HEPEnergyType const ElabNN = particle.getEnergy() / projectileA;
+ HEPEnergyType const sqrtSnn = sqrt(2 * ElabNN * constants::nucleonMass);
+
+ COMBoost const boost{{ElabNN, particle.getMomentum() / projectileA},
+ constants::nucleonMass};
+ CrossSectionType const sigma =
+ composition.getWeightedSum([=](Code const targetId) -> CrossSectionType {
+ return sequence_.getCrossSection(particle, targetId, sqrtSnn);
+ });
+ interaction(secondaries, boost, sqrtSnn, composition, sigma);
sequence_.doSecondaries(secondaries);
- vParticle.erase(); // primary particle is done
+ particle.erase(); // primary particle is done
}
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
inline void Cascade<TTracking, TProcessList, TOutput, TStack>::step(
- particle_type& vParticle) {
+ particle_type& particle) {
+
+ // determine the volume where the particle is (last) known to be
+ auto const* currentLogicalNode = particle.getNode();
+
+ // assert that particle stays outside void Universe if it has no
+ // model properties set
+ assert((currentLogicalNode != &*environment_.getUniverse() ||
+ environment_.getUniverse()->hasModelProperties()) &&
+ "FATAL: The environment model has no valid properties set!");
+
+ NuclearComposition const& composition =
+ currentLogicalNode->getModelProperties().getNuclearComposition();
+
+ // determine sqrtS per nucleon pair, sqrtS_NN
+ Code const projectileId = particle.getPID();
+ unsigned int const projectileA =
+ (is_nucleus(projectileId) ? particle.getNuclearA() : 1);
+ HEPEnergyType const ElabNN = particle.getEnergy() / projectileA;
+ HEPEnergyType const sqrtSnn = sqrt(2 * ElabNN * constants::nucleonMass);
- // determine combined total interaction length (inverse)
- InverseGrammageType const total_inv_lambda =
- sequence_.getInverseInteractionLength(vParticle);
+ // 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(particle, targetId, sqrtSnn);
+ });
+
+ // calculate interaction length in medium
+ GrammageType const total_lambda =
+ (composition.getAverageMassNumber() * constants::u) / total_cx;
// sample random exponential step length in grammage
- ExponentialDistribution expDist(1 / total_inv_lambda);
+ ExponentialDistribution expDist(total_lambda);
GrammageType const next_interact = expDist(rng_);
CORSIKA_LOG_DEBUG(
"total_lambda={} g/cm2, "
", next_interact={} g/cm2",
- double((1. / total_inv_lambda) / 1_g * 1_cm * 1_cm),
+ double(total_lambda / 1_g * 1_cm * 1_cm),
double(next_interact / 1_g * 1_cm * 1_cm));
- auto const* currentLogicalNode = vParticle.getNode();
-
- // assert that particle stays outside void Universe if it has no
- // model properties set
- assert((currentLogicalNode != &*environment_.getUniverse() ||
- environment_.getUniverse()->hasModelProperties()) &&
- "FATAL: The environment model has no valid properties set!");
-
// determine combined total inverse decay time
- InverseTimeType const total_inv_lifetime = sequence_.getInverseLifetime(vParticle);
+ InverseTimeType const total_inv_lifetime = sequence_.getInverseLifetime(particle);
// sample random exponential decay time
ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
@@ -105,11 +153,11 @@ namespace corsika {
(1 / total_inv_lifetime) / 1_s, next_decay / 1_s);
// convert next_decay from time to length [m]
- LengthType const distance_decay = next_decay * vParticle.getMomentum().getNorm() /
- vParticle.getEnergy() * constants::c;
+ LengthType const distance_decay = next_decay * particle.getMomentum().getNorm() /
+ particle.getEnergy() * constants::c;
// determine geometric tracking
- auto [step, nextVol] = tracking_.getTrack(vParticle);
+ auto [step, nextVol] = tracking_.getTrack(particle);
auto geomMaxLength = step.getLength(1);
// convert next_step from grammage to length
@@ -119,7 +167,7 @@ namespace corsika {
// determine the maximum geometric step length
ContinuousProcessStepLength const continuousMaxStep =
- sequence_.getMaxStepLength(vParticle, step);
+ sequence_.getMaxStepLength(particle, step);
LengthType const continuous_max_dist = continuousMaxStep;
// take minimum of geometry, interaction, decay for next step
@@ -150,22 +198,22 @@ namespace corsika {
// move particle along the trajectory to new position
// also update momentum/direction/time
step.setLength(min_distance);
- vParticle.setPosition(step.getPosition(1));
+ particle.setPosition(step.getPosition(1));
// assumption: tracking does not change absolute momentum (continuous physics can and
// will):
- vParticle.setMomentum(step.getDirection(1) * vParticle.getMomentum().getNorm());
+ particle.setMomentum(step.getDirection(1) * particle.getMomentum().getNorm());
// apply all continuous processes on particle + track
- if (sequence_.doContinuous(vParticle, step, limitingId) ==
+ if (sequence_.doContinuous(particle, step, limitingId) ==
ProcessReturn::ParticleAbsorbed) {
CORSIKA_LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
- vParticle.getPID(), vParticle.getEnergy() / 1_GeV);
- if (vParticle.isErased()) {
+ particle.getPID(), particle.getEnergy() / 1_GeV);
+ if (particle.isErased()) {
CORSIKA_LOG_WARN(
"Particle marked as Absorbed in doContinuous, but prematurely erased. This "
"may be bug. Check.");
} else {
- vParticle.erase();
+ particle.erase();
}
return;
}
@@ -188,28 +236,29 @@ namespace corsika {
if (nextVol == environment_.getUniverse().get()) {
CORSIKA_LOG_DEBUG(
"particle left physics world, is now in unknown space -> delete");
- vParticle.erase();
+ particle.erase();
}
- vParticle.setNode(nextVol);
+ particle.setNode(nextVol);
/*
doBoundary may delete the particle (or not)
- caveat: any changes to vParticle, or even the production
+ caveat: any changes to particle, or even the production
of new secondaries is currently not passed to ParticleCut,
thus, particles outside the desired phase space may be produced.
\todo: this must be fixed.
*/
- sequence_.doBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+ sequence_.doBoundaryCrossing(particle, *currentLogicalNode, *nextVol);
return; // step finished
}
CORSIKA_LOG_DEBUG("step limit reached (e.g. deflection). nothing further happens.");
+ // final sanity check, no actions
{
auto const* numericalNodeAfterStep =
- environment_.getUniverse()->getContainingNode(vParticle.getPosition());
+ environment_.getUniverse()->getContainingNode(particle.getPosition());
CORSIKA_LOG_TRACE(
"Geometry check: numericalNodeAfterStep={} currentLogicalNode={}",
fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
@@ -231,23 +280,25 @@ namespace corsika {
// secondaries, b) the projectile particle deleted (or
// changed)
- stack_view_type secondaries(vParticle);
+ stack_view_type secondaries(particle);
/*
Create SecondaryView object on Stack. The data container
remains untouched and identical, and 'projectile' is identical
- to 'vParticle' above this line. However,
- projectile.AddSecondaries populate the SecondaryView, which can
+ to 'particle' above this line. However,
+ projectile.addSecondaries populate the SecondaryView, which can
then be used afterwards for further processing. Thus: it is
- important to use projectile/view (and not vParticle) for Interaction,
+ important to use projectile/view (and not particle) for Interaction,
and Decay!
*/
-
- [[maybe_unused]] auto projectile = secondaries.getProjectile();
-
if (distance_interact < distance_decay) {
- interaction(secondaries, total_inv_lambda);
+ // define boost of NUCLEON-NUCLEON frame
+ COMBoost const boost({ElabNN, particle.getMomentum() / projectileA},
+ constants::nucleonMass);
+ interaction(secondaries, boost, sqrtSnn, composition, total_cx);
} else {
+ [[maybe_unused]] auto projectile = secondaries.getProjectile();
+
if (decay(secondaries, total_inv_lifetime) == ProcessReturn::Decayed) {
if (secondaries.getSize() == 1 &&
projectile.getPID() == secondaries.getNextParticle().getPID()) {
@@ -258,7 +309,7 @@ namespace corsika {
}
sequence_.doSecondaries(secondaries);
- vParticle.erase();
+ particle.erase();
} // namespace corsika
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
@@ -266,19 +317,6 @@ namespace corsika {
stack_view_type& view, InverseTimeType initial_inv_decay_time) {
CORSIKA_LOG_DEBUG("decay");
-#ifdef DEBUG
- InverseTimeType const actual_decay_time = sequence_.getInverseLifetime(view.parent());
- if (actual_decay_time * 0.99 > initial_inv_decay_time) {
- CORSIKA_LOG_WARN(
- "Decay time decreased during step! This leads to un-physical step length. "
- "delta_inverse_decay_time={}",
- (actual_decay_time != InverseTimeType::zero() &&
- initial_inv_decay_time != InverseTimeType::zero()
- ? 1 / initial_inv_decay_time - 1 / actual_decay_time
- : TimeType::zero()));
- }
-#endif
-
// one option is that decay_time is now larger (less
// probability for decay) than it was before the step, thus,
// no decay might actually occur and is allowed
@@ -296,28 +334,20 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TOutput, typename TStack>
inline ProcessReturn Cascade<TTracking, TProcessList, TOutput, TStack>::interaction(
- stack_view_type& view, InverseGrammageType initial_inv_int_length) {
- CORSIKA_LOG_DEBUG("collide");
+ stack_view_type& view, COMBoost const& boost, HEPEnergyType const sqrtSnn,
+ NuclearComposition const& composition,
+ CrossSectionType const initial_cross_section) {
-#ifdef DEBUG
- InverseGrammageType const actual_inv_length = sequence_.getInverseInteractionLength(
- view.parent()); // 1/lambda_int after step, -dE/dX etc.
-
- if (actual_inv_length * 0.99 > initial_inv_int_length) {
- CORSIKA_LOG_WARN(
- "Interaction length decreased during step! This leads to un-physical step "
- "length. delta_inverse_interaction_length={}",
- 1 / initial_inv_int_length - 1 / actual_inv_length);
- }
-#endif
+ CORSIKA_LOG_DEBUG("collide");
- // one option is that interaction_length is now larger (less
+ // one option is that cross section is now smaller (less
// probability for collision) than it was before the step, thus,
// no interaction might actually occur and is allowed
- UniformRealDistribution<InverseGrammageType> uniDist(initial_inv_int_length);
- const auto sample_process = uniDist(rng_);
- auto const returnCode = sequence_.selectInteraction(view, sample_process);
+ UniformRealDistribution<CrossSectionType> uniDist(initial_cross_section);
+ CrossSectionType const sample_process_by_cx = uniDist(rng_);
+ auto const returnCode = sequence_.selectInteraction(view, boost, sqrtSnn, composition,
+ rng_, sample_process_by_cx);
if (returnCode != ProcessReturn::Interacted) {
CORSIKA_LOG_DEBUG("Particle did not interact!");
}
diff --git a/corsika/detail/framework/process/InteractionProcess.hpp b/corsika/detail/framework/process/InteractionProcess.hpp
index 2446a385d..baa15a8bc 100644
--- a/corsika/detail/framework/process/InteractionProcess.hpp
+++ b/corsika/detail/framework/process/InteractionProcess.hpp
@@ -14,10 +14,12 @@
namespace corsika {
/**
- traits test for InteractionProcess::doInteraction method
- */
+ * @file InteractionProcess.hpp
+ *
+ * traits test for InteractionProcess::doInteraction methods etc.
+ */
- template <class TProcess, typename TReturn, typename... TArgs>
+ template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
struct has_method_doInteract : public detail::has_method_signature<TReturn, TArgs...> {
///! method signature
@@ -29,7 +31,7 @@ namespace corsika {
//! signature of templated method
template <class T>
- static decltype(testSignature(&T::template doInteraction<TArgs...>)) test(
+ static decltype(testSignature(&T::template doInteraction<TTemplate>)) test(
std::nullptr_t);
//! signature of non-templated method
@@ -46,11 +48,10 @@ namespace corsika {
//! @}
};
- //! @file InteractionProcess.hpp
//! value traits type
- template <class TProcess, typename TReturn, typename... TArgs>
+ template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
bool constexpr has_method_doInteract_v =
- has_method_doInteract<TProcess, TReturn, TArgs...>::value;
+ has_method_doInteract<TProcess, TReturn, TTemplate, TArgs...>::value;
/**
traits test for InteractionProcess::getInteractionLength method
@@ -86,11 +87,48 @@ namespace corsika {
//! @}
};
- //! @file InteractionProcess.hpp
- //! value traits type
-
+ //! value traits type shortcut
template <class TProcess, typename TReturn, typename... TArgs>
bool constexpr has_method_getInteractionLength_v =
has_method_getInteractionLength<TProcess, TReturn, TArgs...>::value;
+ /**
+ traits test for InteractionProcess::getCrossSection method
+ */
+
+ template <class TProcess, typename TReturn, typename... TArgs>
+ struct has_method_getCrossSection
+ : public detail::has_method_signature<TReturn, TArgs...> {
+
+ ///! method signature
+ using detail::has_method_signature<TReturn, TArgs...>::testSignature;
+
+ //! the default value
+ template <class T>
+ static std::false_type test(...);
+
+ //! templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::template getCrossSection<TArgs...>)) test(
+ std::nullptr_t);
+
+ //! non templated parameter option
+ template <class T>
+ static decltype(testSignature(&T::getCrossSection)) test(std::nullptr_t);
+
+ public:
+ /**
+ @name traits results
+ @{
+ */
+ using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+ static const bool value = type::value;
+ //! @}
+ };
+
+ //! value traits type shortcut
+ template <class TProcess, typename TReturn, typename... TArgs>
+ bool constexpr has_method_getCrossSection_v =
+ has_method_getCrossSection<TProcess, TReturn, TArgs...>::value;
+
} // namespace corsika
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 4d381abb1..70e0a6f7b 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -9,6 +9,7 @@
#pragma once
#include <corsika/framework/core/PhysicalUnits.hpp>
+
#include <corsika/framework/process/BaseProcess.hpp>
#include <corsika/framework/process/BoundaryCrossingProcess.hpp>
#include <corsika/framework/process/ContinuousProcess.hpp>
@@ -259,7 +260,8 @@ namespace corsika {
template <typename TParticle, typename TTrack>
inline ContinuousProcessStepLength
ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
- IndexProcess2>::getMaxStepLength(TParticle& particle, TTrack& vTrack) {
+ IndexProcess2>::getMaxStepLength(TParticle&& particle,
+ TTrack&& vTrack) {
// if no other process in the sequence implements it
ContinuousProcessStepLength max_length(std::numeric_limits<double>::infinity() *
meter);
@@ -309,24 +311,34 @@ namespace corsika {
template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
int IndexProcess2>
template <typename TParticle>
- inline InverseGrammageType
+ inline CrossSectionType
ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
- IndexProcess2>::getInverseInteractionLength(TParticle&& particle) {
+ IndexProcess2>::getCrossSection(TParticle&& projectile,
+ Code const targetId,
+ HEPEnergyType const sqrtSnn) const {
- InverseGrammageType tot = 0 * meter * meter / gram; // default value
+ CrossSectionType tot = CrossSectionType::zero();
if constexpr (is_process_v<process1_type>) { // to protect from further compiler
// errors if process1_type is invalid
- if constexpr (is_interaction_process_v<process1_type> ||
- process1_type::is_process_sequence) {
- tot += A_.getInverseInteractionLength(particle);
+ if constexpr (is_interaction_process_v<process1_type>) {
+ Code const projectileId = projectile.getPID();
+ tot += A_.getCrossSection(projectileId, targetId, sqrtSnn,
+ is_nucleus(projectileId) ? projectile.getNuclearA() : 1,
+ is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+ } else if constexpr (process1_type::is_process_sequence) {
+ tot += A_.getCrossSection(projectile, targetId, sqrtSnn);
}
}
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_.getInverseInteractionLength(particle);
+ if constexpr (is_interaction_process_v<process2_type>) {
+ Code const projectileId = projectile.getPID();
+ tot += B_.getCrossSection(projectileId, targetId, sqrtSnn,
+ is_nucleus(projectileId) ? projectile.getNuclearA() : 1,
+ is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
+ } else if constexpr (process2_type::is_process_sequence) {
+ tot += B_.getCrossSection(projectile, targetId, sqrtSnn);
}
}
return tot;
@@ -410,62 +422,138 @@ namespace corsika {
template <typename TSecondaryView>
inline ProcessReturn
ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1, IndexProcess2>::
- selectInteraction(TSecondaryView& view,
- [[maybe_unused]] InverseGrammageType lambda_inv_select,
- [[maybe_unused]] InverseGrammageType lambda_inv_sum) {
+ selectInteraction(TSecondaryView&& view, COMBoost const& boost,
+ [[maybe_unused]] HEPEnergyType const sqrtSnn,
+ [[maybe_unused]] NuclearComposition const& composition,
+ [[maybe_unused]] TRNG&& rng,
+ [[maybe_unused]] CrossSectionType const cx_select,
+ [[maybe_unused]] CrossSectionType cx_sum) {
- // TODO: add check for lambda_inv_select > lambda_inv_tot
+ // TODO: add check for cx_select > cx_tot
if constexpr (is_process_v<process1_type>) { // to protect from further compiler
// 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, lambda_inv_select, lambda_inv_sum);
+ ProcessReturn const ret = A_.selectInteraction(view, boost, sqrtSnn, 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>) {
- // if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
+
+ auto const& projectile = view.parent();
+ Code const projectileId = projectile.getPID();
+ unsigned int const projectileA =
+ is_nucleus(projectileId) ? projectile.getNuclearA() : 1;
+
+ // get cross section vector for all material components
+ static_assert(
+ has_method_getCrossSection_v<TProcess1, // process object
+ CrossSectionType, // return type
+ Code, // parameters
+ Code, HEPEnergyType, unsigned int, unsigned int>,
+ "TProcess1 has no method with correct signature \"CrossSectionType "
+ "getCrossSection(Code, Code, HEPEnergyType, unsigned int, unsigned int"
+ ")\" required by InteractionProcess<TProcess1>. ");
+
+ 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) : 1);
+ });
+
+ cx_sum += std::accumulate(weightedCrossSections.cbegin(),
+ weightedCrossSections.cend(), CrossSectionType::zero());
+
// check if we should execute THIS process and then EXIT
- if (lambda_inv_select <= lambda_inv_sum) {
+ if (cx_select <= cx_sum) {
+
+ // now also sample targetId from weighted cross sections
+ Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
// interface checking on TProcess1
- static_assert(has_method_doInteract_v<TProcess1, void, TSecondaryView&>,
- "TDerived has no method with correct signature \"void "
- "doInteraction(TSecondaryView&)\" required for "
- "InteractionProcess<TDerived>. ");
+ static_assert(
+ has_method_doInteract_v<TProcess1, // process object
+ void, // return type
+ TSecondaryView, // template argument
+ TSecondaryView&, // method parameters
+ COMBoost const&, Code, Code, HEPEnergyType,
+ unsigned int, unsigned int>,
+ "TProcess1 has no method with correct signature \"void "
+ "doInteraction<TSecondaryView>(TSecondaryView&, COMBoost&, Code, "
+ "Code, HEPEnergyType, unsigned int, unsigned int)\" required for "
+ "InteractionProcess<TProcess1>. ");
+
+ A_.template doInteraction(view, boost, projectileId, targetId, sqrtSnn,
+ projectileA,
+ is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
- A_.template doInteraction(view);
return ProcessReturn::Interacted;
}
- } // end branch A
- }
+ }
+ } // end branch A
if constexpr (is_process_v<process2_type>) { // to protect from further compiler
// errors if process2_type is invalid
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
- return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
+ return B_.selectInteraction(view, boost, sqrtSnn, composition, rng, cx_select,
+ cx_sum);
} else if constexpr (is_interaction_process_v<process2_type>) {
- // if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
- // soon as SecondaryView::parent() is migrated!
- // check if we should execute THIS process and then EXIT
- if (lambda_inv_select <= lambda_inv_sum) {
- // interface checking on TProcess1
- static_assert(has_method_doInteract_v<TProcess2, void, TSecondaryView&>,
- "TDerived has no method with correct signature \"void "
- "doInteraction(TSecondaryView&)\" required for "
- "InteractionProcess<TDerived>. ");
+ auto const& projectile = view.parent();
+ Code const projectileId = projectile.getPID();
+ unsigned int const projectileA =
+ is_nucleus(projectileId) ? projectile.getNuclearA() : 1;
+
+ // get cross section vector for all material components
+ static_assert(
+ has_method_getCrossSection_v<TProcess2, // process object
+ CrossSectionType, // return type
+ Code, // parameters
+ Code, HEPEnergyType, unsigned int, unsigned int>,
+ "TProcess2 has no method with correct signature \"CrossSectionType "
+ "getCrossSection(Code, Code, HEPEnergyType, unsigned int, unsigned int"
+ ")\" required by InteractionProcess<TProcess1>. ");
+
+ 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) : 1);
+ });
+
+ cx_sum += std::accumulate(weightedCrossSections.begin(),
+ weightedCrossSections.end(), CrossSectionType::zero());
+
+ // 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);
+
+ // interface checking on TProcess2
+ static_assert(
+ has_method_doInteract_v<TProcess2, // process object
+ void, // return type
+ TSecondaryView, // template argument
+ TSecondaryView&, // method parameters
+ COMBoost const&, Code, Code, HEPEnergyType,
+ unsigned int, unsigned int>,
+ "TProcess1 has no method with correct signature \"void "
+ "doInteraction<TSecondaryView>(TSecondaryView&, COMBoost&, Code, "
+ "Code, HEPEnergyType, unsigned int, unsigned int)\" required for "
+ "InteractionProcess<TProcess2>. ");
+
+ B_.doInteraction(view, boost, projectileId, targetId, sqrtSnn, projectileA,
+ is_nucleus(targetId) ? get_nucleus_A(targetId) : 1);
- B_.doInteraction(view);
return ProcessReturn::Interacted;
}
- } // end branch B_
- }
+ }
+ } // end branch B_
return ProcessReturn::Ok;
}
@@ -501,7 +589,7 @@ namespace corsika {
template <typename TSecondaryView>
inline ProcessReturn ProcessSequence<
TProcess1, TProcess2, IndexStart, IndexProcess1,
- IndexProcess2>::selectDecay(TSecondaryView& view,
+ IndexProcess2>::selectDecay(TSecondaryView&& view,
[[maybe_unused]] InverseTimeType decay_inv_select,
[[maybe_unused]] InverseTimeType decay_inv_sum) {
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index 634395f55..46e9bfbd2 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -210,80 +210,161 @@ namespace corsika {
template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TParticle>
- inline InverseGrammageType SwitchProcessSequence<
+ CrossSectionType SwitchProcessSequence<
TCondition, TSequence, USequence, IndexStart, IndexProcess1,
- IndexProcess2>::getInverseInteractionLength(TParticle&& particle) {
-
- if (select_(particle)) {
- if constexpr (is_interaction_process_v<process1_type> ||
- process1_type::is_process_sequence) {
- return A_.getInverseInteractionLength(particle);
+ IndexProcess2>::getCrossSection(TParticle const& projectile, Code const targetId,
+ HEPEnergyType const sqrtSnn) const {
+
+ if (select_(projectile.parent())) {
+ if constexpr (is_interaction_process_v<process1_type>) {
+ return A_.getCrossSection(projectile.getPID(), targetId, sqrtSnn,
+ projectile.getNuclearA(),
+ is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ } else if (process1_type::is_process_sequence) {
+ return A_.getCrossSection(projectile, targetId, sqrtSnn,
+ is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
}
} else {
-
- if constexpr (is_interaction_process_v<process2_type> ||
- process2_type::is_process_sequence) {
- return B_.getInverseInteractionLength(particle);
+ if constexpr (is_interaction_process_v<process2_type>) {
+ return B_.getCrossSection(projectile.getPID(), targetId, sqrtSnn,
+ projectile.getNuclearA(),
+ is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ } else if (process2_type::is_process_sequence) {
+ return B_.getCrossSection(projectile, targetId, sqrtSnn,
+ is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
}
}
- return 0 * meter * meter / gram; // default value
+ return CrossSectionType::zero(); // default value
}
template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
- template <typename TSecondaryView>
- inline ProcessReturn SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
- IndexProcess1, IndexProcess2>::
- selectInteraction(TSecondaryView& view,
- [[maybe_unused]] InverseGrammageType lambda_inv_select,
- [[maybe_unused]] InverseGrammageType lambda_inv_sum) {
+ template <typename TSecondaryView, typename TRNG>
+ inline ProcessReturn SwitchProcessSequence<
+ TCondition, TSequence, USequence, IndexStart, IndexProcess1,
+ IndexProcess2>::selectInteraction(TSecondaryView& view, COMBoost const& boost,
+ HEPEnergyType const sqrtSnn,
+ NuclearComposition const& composition, TRNG& rng,
+ [[maybe_unused]] CrossSectionType const cx_select,
+ [[maybe_unused]] CrossSectionType cx_sum) {
+
if (select_(view.parent())) {
if constexpr (process1_type::is_process_sequence) {
// if A_ is a process sequence --> check inside
- ProcessReturn const ret =
- A_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
- // if A_ did succeed, stop routine. Not checking other static branch B_.
- if (ret != ProcessReturn::Ok) { return ret; }
+ return A_.selectInteraction(view, boost, sqrtSnn, composition, rng, cx_select,
+ cx_sum);
} else if constexpr (is_interaction_process_v<process1_type>) {
- // if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
- // check if we should execute THIS process and then EXIT
- if (lambda_inv_select < lambda_inv_sum) {
- // interface checking on TSequence
- static_assert(has_method_doInteract_v<TSequence, void, TSecondaryView&>,
- "TDerived has no method with correct signature \"void "
- "doInteraction(TSecondaryView&)\" required for "
- "InteractionProcess<TDerived>. ");
+ auto const& projectile = view.parent();
+ Code const projectileId = projectile.getPID();
+ unsigned int const projectileA = projectile.getNuclearA();
+
+ // get cross section vector for all material components
+ static_assert(
+ has_method_getCrossSection_v<TSequence, // process object
+ CrossSectionType, // return type
+ Code, // parameters
+ Code, HEPEnergyType, unsigned int, unsigned int>,
+ "TSequence has no method with correct signature \"CrossSectionType "
+ "getCrossSection(Code, Code, HEPEnergyType, unsigned int, unsigned int"
+ ")\" 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);
+ });
+
+ cx_sum += std::accumulate(weightedCrossSections.cbegin(),
+ weightedCrossSections.cend(), CrossSectionType::zero());
+ if (cx_select < cx_sum) {
+
+ // now also sample targetId from weighted cross sections
+ Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+
+ // 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);
- A_.doInteraction(view);
return ProcessReturn::Interacted;
- }
- } // end branch A_
- } else {
+ } // end collision branch A
+ }
+
+ } else { // selection: end branch A, start branch B
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
- return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
+ return B_.selectInteraction(view, boost, sqrtSnn, composition, rng, cx_select,
+ cx_sum);
} else if constexpr (is_interaction_process_v<process2_type>) {
- // if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
- // check if we should execute THIS process and then EXIT
- if (lambda_inv_select < lambda_inv_sum) {
- // interface checking on TSequence
- static_assert(has_method_doInteract_v<USequence, void, TSecondaryView&>,
- "TDerived has no method with correct signature \"void "
- "doInteraction(TSecondaryView&)\" required for "
- "InteractionProcess<TDerived>. ");
+ auto const& projectile = view.parent();
+ Code const projectileId = projectile.getPID();
+ unsigned int const projectileA = projectile.getNuclearA();
+
+ // get cross section vector for all material components
+ static_assert(
+ has_method_getCrossSection_v<USequence, // process object
+ CrossSectionType, // return type
+ Code, // parameters
+ Code, HEPEnergyType, unsigned int, unsigned int>,
+ "USequence has no method with correct signature \"CrossSectionType "
+ "getCrossSection(Code, Code, HEPEnergyType, unsigned int, unsigned int"
+ ")\" 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);
+ });
+
+ cx_sum += std::accumulate(weightedCrossSections.cbegin(),
+ weightedCrossSections.cend(), CrossSectionType::zero());
+
+ if (cx_select < cx_sum) {
+
+ // now also sample targetId from weighted cross sections
+ Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+
+ // 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);
- B_.doInteraction(view);
return ProcessReturn::Interacted;
- }
- } // end branch B_
- }
+ } // end collision in branch B
+ }
+ } // end branch B_
+
return ProcessReturn::Ok;
}
diff --git a/corsika/detail/framework/utility/COMBoost.inl b/corsika/detail/framework/utility/COMBoost.inl
index 9b05bef0e..842a6f99e 100644
--- a/corsika/detail/framework/utility/COMBoost.inl
+++ b/corsika/detail/framework/utility/COMBoost.inl
@@ -39,7 +39,6 @@ namespace corsika {
auto const coshEta = sqrt(1 + pProjNormSquared / s);
setBoost(coshEta, sinhEta);
-
CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
coshEta, boost_.determinant() - 1);
}
@@ -52,6 +51,8 @@ namespace corsika {
auto const sinhEta = -norm / mass;
auto const coshEta = sqrt(1 + squaredNorm / (mass * mass));
setBoost(coshEta, sinhEta);
+ CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+ coshEta, boost_.determinant() - 1);
}
template <typename FourVector>
@@ -92,15 +93,13 @@ namespace corsika {
Vector<typename decltype(pCM)::dimension_type> pLab{rotatedCS_, pCM};
pLab.rebase(originalCS_);
- FourVector f(E_lab, pLab);
-
CORSIKA_LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
" plab={} GeV (norm={} GeV) "
" GeV), invariant mass = {}",
- E_lab / 1_GeV, f.getNorm() / 1_GeV, pLab.getComponents(),
- pLab.getNorm() / 1_GeV);
+ E_lab / 1_GeV, FourVector{E_lab, pLab}.getNorm() / 1_GeV,
+ pLab.getComponents(), pLab.getNorm() / 1_GeV);
- return f;
+ return FourVector{E_lab, pLab};
}
inline void COMBoost::setBoost(double coshEta, double sinhEta) {
@@ -110,4 +109,6 @@ namespace corsika {
inline CoordinateSystemPtr COMBoost::getRotatedCS() const { return rotatedCS_; }
+ inline CoordinateSystemPtr COMBoost::getOriginalCS() const { return originalCS_; }
+
} // namespace corsika
diff --git a/corsika/detail/media/NuclearComposition.inl b/corsika/detail/media/NuclearComposition.inl
index 4995799d5..0237d73e9 100644
--- a/corsika/detail/media/NuclearComposition.inl
+++ b/corsika/detail/media/NuclearComposition.inl
@@ -23,31 +23,52 @@
namespace corsika {
inline NuclearComposition::NuclearComposition(std::vector<Code> const& pComponents,
- std::vector<float> const& pFractions)
+ std::vector<double> const& pFractions)
: numberFractions_(pFractions)
, components_(pComponents)
- , avgMassNumber_(std::inner_product(
- pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
- std::plus<double>(), [](auto const compID, auto const fraction) -> double {
- if (is_nucleus(compID)) {
- return get_nucleus_A(compID) * fraction;
- } else {
- return get_mass(compID) / convert_SI_to_HEP(constants::u) * fraction;
- }
- })) {
- assert(pComponents.size() == pFractions.size());
- auto const sumFractions =
- std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
-
- if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
+ , avgMassNumber_(getWeightedSum([](Code const compID) -> double {
+ if (is_nucleus(compID)) {
+ return get_nucleus_A(compID);
+ } else {
+ return get_mass(compID) / convert_SI_to_HEP(constants::u);
+ }
+ })) {
+ if (pComponents.size() != pFractions.size()) {
+ throw std::runtime_error(
+ "Cannot construct NuclearComposition from vectors of different sizes.");
+ }
+ auto const sumFractions = std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.);
+
+ if (!(0.999 < sumFractions && sumFractions < 1.001)) {
throw std::runtime_error("element fractions do not add up to 1");
}
this->updateHash();
}
template <typename TFunction>
- inline auto NuclearComposition::getWeightedSum(TFunction const& func) const {
- using ResultQuantity = decltype(func(*components_.cbegin()));
+ inline auto NuclearComposition::getWeighted(TFunction const& func) const {
+ using ResultQuantity = decltype(func(std::declval<Code>()));
+ auto const product = [&](auto const compID, auto const fraction) {
+ return func(compID) * fraction;
+ };
+
+ if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
+ std::vector<ResultQuantity> result(components_.size(), ResultQuantity::zero());
+ std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+ result.begin(), product);
+ return result;
+ } else {
+ std::vector<ResultQuantity> result(components_.size(), ResultQuantity(0));
+ std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+ result.begin(), product);
+ return result;
+ }
+ } // namespace corsika
+
+ template <typename TFunction>
+ inline auto NuclearComposition::getWeightedSum(TFunction const& func) const
+ -> decltype(func(std::declval<Code>())) {
+ using ResultQuantity = decltype(func(std::declval<Code>()));
auto const prod = [&](auto const compID, auto const fraction) {
return func(compID) * fraction;
@@ -68,7 +89,7 @@ namespace corsika {
inline size_t NuclearComposition::getSize() const { return numberFractions_.size(); }
- inline std::vector<float> const& NuclearComposition::getFractions() const {
+ inline std::vector<double> const& NuclearComposition::getFractions() const {
return numberFractions_;
}
@@ -82,16 +103,18 @@ namespace corsika {
template <class TRNG>
inline Code NuclearComposition::sampleTarget(std::vector<CrossSectionType> const& sigma,
- TRNG& randomStream) const {
-
- assert(sigma.size() == numberFractions_.size());
+ TRNG&& randomStream) const {
+ if (sigma.size() != numberFractions_.size()) {
+ throw std::runtime_error("incompatible vector sigma as input");
+ }
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 //<decltype(numberFractions_.begin()),
+ // decltype(sigma.begin())>
+ (numberFractions_.begin(), sigma.begin()),
+ WeightProviderIterator //<decltype(numberFractions_.begin()),
+ // decltype(sigma.end())>
+ (numberFractions_.end(), sigma.end()));
auto const iChannel = channelDist(randomStream);
return components_[iChannel];
@@ -99,6 +122,7 @@ namespace corsika {
// Note: when this class ever modifies its internal data, the hash
// must be updated, too!
+ // the hash value is important to find tables, etc.
inline size_t NuclearComposition::getHash() const { return hash_; }
inline bool NuclearComposition::operator==(NuclearComposition const& v) const {
@@ -107,7 +131,8 @@ namespace corsika {
inline void NuclearComposition::updateHash() {
std::vector<std::size_t> hashes;
- for (float ifrac : this->getFractions()) hashes.push_back(std::hash<float>{}(ifrac));
+ for (double ifrac : this->getFractions())
+ hashes.push_back(std::hash<double>{}(ifrac));
for (Code icode : this->getComponents())
hashes.push_back(std::hash<int>{}(static_cast<int>(icode)));
std::size_t h = std::hash<double>{}(this->getAverageMassNumber());
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
deleted file mode 100644
index 8c424a7db..000000000
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ /dev/null
@@ -1,351 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/modules/sibyll/Interaction.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/modules/sibyll/ParticleConversion.hpp>
-#include <corsika/modules/sibyll/SibStack.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-
-#include <sibyll2.3d.hpp>
-
-#include <tuple>
-
-namespace corsika::sibyll {
-
- inline Interaction::Interaction(const bool sibyll_printout_on)
- : sibyll_listing_(sibyll_printout_on) {
- // initialize Sibyll
- static bool initialized = false;
- if (!initialized) {
- sibyll_ini_();
- initialized = true;
- }
- }
-
- inline Interaction::~Interaction() {
- CORSIKA_LOG_DEBUG("Sibyll::Interaction n={}, Nnuc={}", count_, nucCount_);
- }
-
- inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
- Interaction::getCrossSection(corsika::Code const BeamId, corsika::Code const TargetId,
- corsika::HEPEnergyType const CoMenergy) const {
- double sigProd, sigEla, dummy, dum1, dum3, dum4;
- double dumdif[3];
- int const iBeam = corsika::sibyll::getSibyllXSCode(
- BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
- if (!iBeam) {
- throw std::runtime_error(
- fmt::format("Interaction of beam {} not defined in "
- "Sibyll!",
- BeamId));
- }
- if (!isValidCoMEnergy(CoMenergy)) {
- throw std::runtime_error(
- "Interaction: getCrossSection: CoM energy outside range for Sibyll!");
- }
- double const dEcm = CoMenergy / 1_GeV;
- // single nucleon target (p,n, hydrogen) or 4<=A<=18
- if (isValidTarget(TargetId)) {
- // single nucleon target
- if (TargetId == corsika::Code::Proton || TargetId == Code::Hydrogen ||
- TargetId == Code::Neutron) {
- sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
- } else {
- // nuclear target
- int const iTarget = corsika::get_nucleus_A(TargetId);
- sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
- }
- } else {
- // throw std::runtime_error(
- // "Sibyll nuclear target outside range. Only nuclei with 4<=A<18 are
- // allowed.");
-
- // no interaction in sibyll possible, return infinite cross section? or throw?
- sigProd = std::numeric_limits<double>::infinity();
- sigEla = std::numeric_limits<double>::infinity();
- }
- return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
- }
-
- template <typename TParticle>
- inline corsika::GrammageType Interaction::getInteractionLength(
- TParticle const& projectile) const {
-
- corsika::Code const corsikaBeamId = projectile.getPID();
-
- // beam corsika for sibyll : 1, 2, 3 for p, pi, k
- // read from cross section code table
- bool const kInteraction = corsika::sibyll::canInteract(corsikaBeamId);
-
- MomentumVector const& pLab = projectile.getMomentum();
- CoordinateSystemPtr const& labCS = pLab.getCoordinateSystem();
-
- // FOR NOW: assume target is at rest
- MomentumVector pTarget(labCS, {0_GeV, 0_GeV, 0_GeV});
-
- // total momentum and energy
- HEPEnergyType Elab = projectile.getEnergy() + constants::nucleonMass;
- MomentumVector pTotLab(labCS, {0_GeV, 0_GeV, 0_GeV});
- pTotLab += pLab;
- pTotLab += pTarget;
- auto const pTotLabNorm = pTotLab.getNorm();
- // calculate cm. energy
- HEPEnergyType const ECoM = sqrt(
- (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
-
- CORSIKA_LOG_DEBUG(
- "Interaction: LambdaInt: \n"
- " input energy: {} GeV "
- " beam can interact: {} "
- " beam pid: {}",
- projectile.getEnergy() / 1_GeV, kInteraction, projectile.getPID());
-
- // TODO: move limits into variables
- // FR: removed && Elab >= 8.5_GeV
- if (kInteraction && isValidCoMEnergy(ECoM)) {
-
- // get target from environment
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
-
- auto const* currentNode = projectile.getNode();
- auto const& mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
-
- si::CrossSectionType weightedProdCrossSection = mediumComposition.getWeightedSum(
- [=](corsika::Code targetID) -> si::CrossSectionType {
- // Argon needs special handling ....
- return targetID == Code::Argon ? CrossSectionType::zero()
- : std::get<0>(this->getCrossSection(
- corsikaBeamId, targetID, ECoM));
- });
-
- CORSIKA_LOG_DEBUG(
- "Interaction: "
- "IntLength: weighted CrossSection (mb): {} ",
- weightedProdCrossSection / 1_mb);
-
- // calculate interaction length in medium
- GrammageType const int_length = mediumComposition.getAverageMassNumber() *
- constants::u / weightedProdCrossSection;
- CORSIKA_LOG_DEBUG(
- "Interaction: "
- "interaction length (g/cm2): {} ",
- int_length / (0.001_kg) * 1_cm * 1_cm);
-
- return int_length;
- }
-
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
-
- /*
- * In this function SIBYLL is called to produce one event. The
- * event is copied (and boosted) into the shower lab frame.
- */
-
- template <typename TSecondaryView>
- inline void Interaction::doInteraction(TSecondaryView& view) {
-
- auto const projectile = view.getProjectile();
- auto const corsikaBeamId = projectile.getPID();
-
- if (corsika::is_nucleus(corsikaBeamId)) {
- // nuclei handled by different process, this should not happen
- throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
- }
-
- // position and time of interaction, not used in Sibyll
- Point const pOrig = projectile.getPosition();
- TimeType const tOrig = projectile.getTime();
-
- // define projectile
- HEPEnergyType const eProjectileLab = projectile.getEnergy();
- auto const pProjectileLab = projectile.getMomentum();
- CoordinateSystemPtr const& originalCS = pProjectileLab.getCoordinateSystem();
-
- CORSIKA_LOG_DEBUG(
- "ProcessSibyll: "
- "DoInteraction: pid {} interaction ",
- corsikaBeamId);
-
- // define target
- // for Sibyll is always a single nucleon
- // FOR NOW: target is always at rest
- auto const eTargetLab = 0_GeV + constants::nucleonMass;
- auto const pTargetLab = MomentumVector(originalCS, 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());
- 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 Sibyll projectile: +z
- COMBoost const boost(PprojLab, constants::nucleonMass);
- auto const& csPrime = boost.getRotatedCS();
-
- // just for show:
- // boost projecticle
- [[maybe_unused]] auto const PprojCoM = boost.toCoM(PprojLab);
- // boost target
- [[maybe_unused]] auto const PtargCoM = boost.toCoM(PtargLab);
- CORSIKA_LOG_DEBUG(
- "Interaction: ebeam CoM: {} GeV "
- "Interaction: pbeam CoM: {} GeV ",
- PprojCoM.getTimeLikeComponent() / 1_GeV,
- PprojCoM.getSpaceLikeComponents().getComponents(csPrime) / 1_GeV);
- CORSIKA_LOG_DEBUG(
- "Interaction: etarget CoM: {} GeV "
- "Interaction: ptarget CoM: {} GeV ",
- PtargCoM.getTimeLikeComponent() / 1_GeV,
- PtargCoM.getSpaceLikeComponents().getComponents(csPrime) / 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());
-
- // sample target mass number
- auto const* currentNode = projectile.getNode();
- auto const& mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
- // get cross sections for target materials
- /*
- Here we read the cross section from the interaction model again,
- should be passed from getInteractionLength if possible
- */
- //#warning reading interaction cross section again, should not be necessary
- auto const& compVec = mediumComposition.getComponents();
- std::vector<CrossSectionType> cross_section_of_components(compVec.size());
-
- for (size_t i = 0; i < compVec.size(); ++i) {
- auto const targetId = compVec[i];
- if (targetId == Code::Argon) continue; // skip Argon ....
- const auto [sigProd, sigEla] = getCrossSection(corsikaBeamId, targetId, Ecm);
- [[maybe_unused]] const auto& dummy_sigEla = sigEla;
- cross_section_of_components[i] = sigProd;
- }
-
- auto const targetCode =
- mediumComposition.sampleTarget(cross_section_of_components, RNG_);
- CORSIKA_LOG_DEBUG("Interaction: target selected: {} ", targetCode);
- /*
- FOR NOW: allow nuclei with A<18 or protons only.
- when medium composition becomes more complex, approximations will have to be
- allowed air in atmosphere also contains some Argon.
- */
- int targetSibCode = -1;
- if (is_nucleus(targetCode)) targetSibCode = get_nucleus_A(targetCode);
- if (targetCode == Proton::code) targetSibCode = 1;
- CORSIKA_LOG_DEBUG("Interaction: sibyll code: {}", targetSibCode);
- if (targetSibCode > int(maxTargetMassNumber_) || targetSibCode < 1)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 or protons are "
- "allowed.");
-
- // beam id for sibyll
- int const kBeam = corsika::sibyll::convertToSibyllRaw(corsikaBeamId);
-
- CORSIKA_LOG_DEBUG(
- "Interaction: "
- " DoInteraction: E(GeV): {} "
- " Ecm(GeV): {} ",
- eProjectileLab / 1_GeV, Ecm / 1_GeV);
- if (Ecm > getMaxEnergyCoM())
- throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
- // FR: removed eProjectileLab < 8.5_GeV ||
- if (Ecm < getMinEnergyCoM()) {
- CORSIKA_LOG_DEBUG(
- "Interaction: "
- " DoInteraction: should have dropped particle.. "
- "THIS IS AN ERROR");
- throw std::runtime_error("energy too low for SIBYLL");
- } else {
- count_++;
- // Sibyll does not know about units..
- const double sqs = Ecm / 1_GeV;
- // running sibyll, filling stack
- sibyll_(kBeam, targetSibCode, sqs);
-
- if (sibyll_listing_) {
- // print final state
- int print_unit = 6;
- sib_list_(print_unit);
- nucCount_ += get_nwounded() - 1;
- }
-
- // add particles from sibyll to stack
- // link to sibyll stack
- SibStack ss;
-
- MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
- for (auto& psib : ss) {
-
- // abort on particles that have decayed in Sibyll. Should not happen!
- if (psib.hasDecayed())
- throw std::runtime_error("found particle that decayed in SIBYLL!");
-
- // transform 4-momentum to lab. frame
- // note that the momentum needs to be rotated back
- auto const tmp = psib.getMomentum().getComponents();
- auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp);
- HEPEnergyType const eCoM = psib.getEnergy();
- auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
- auto const p3lab = Plab.getSpaceLikeComponents();
- assert(p3lab.getCoordinateSystem() == originalCS); // just to be sure!
-
- class KinematicParticle {};
-
- // add to corsika stack
- auto pnew = view.addSecondary(std::make_tuple(
- corsika::sibyll::convertFromSibyll(psib.getPID()), p3lab, pOrig, tOrig));
-
- Plab_final += pnew.getMomentum();
- Elab_final += pnew.getEnergy();
- Ecm_final += psib.getEnergy();
- }
- CORSIKA_LOG_DEBUG(
- "conservation (all GeV): "
- "Ecm_initial(per nucleon)={:.2f}, Ecm_final(per nucleon)={:.2f}, "
- "Elab_initial={:.2f}, Elab_final={:.2f}, "
- "Elab-diff (%)={:.2f}, "
- "m in target nucleons={:.2f}, "
- "Plab_initial={:.2f}, "
- "Plab_final={:.2f} ",
- Ecm / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Etot / 1_GeV,
- Elab_final / 1_GeV,
- (Elab_final / (Etot + get_nwounded() * constants::nucleonMass) - 1) * 100,
- constants::nucleonMass * get_nwounded() / 1_GeV,
- (pProjectileLab / 1_GeV).getComponents(), (Plab_final / 1_GeV).getComponents());
- }
- }
-
-} // namespace corsika::sibyll
diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
new file mode 100644
index 000000000..1af79a9d2
--- /dev/null
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -0,0 +1,183 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+
+#include <corsika/modules/sibyll/ParticleConversion.hpp>
+#include <corsika/modules/sibyll/SibStack.hpp>
+
+#include <sibyll2.3d.hpp>
+
+#include <tuple>
+
+namespace corsika::sibyll {
+
+ inline void InteractionModel::setVerbose(bool const flag) { sibyll_listing_ = flag; }
+
+ inline InteractionModel::InteractionModel()
+ : sibyll_listing_(false) {
+ // initialize Sibyll
+ static bool initialized = false;
+ if (!initialized) {
+ sibyll_ini_();
+ initialized = true;
+ }
+ }
+
+ inline InteractionModel::~InteractionModel() {
+ CORSIKA_LOG_DEBUG("Sibyll::Model n={}, Nnuc={}", count_, nucCount_);
+ }
+
+ inline void constexpr InteractionModel::isValid(Code const projectileId,
+ Code const targetId,
+ HEPEnergyType const sqrtSnn,
+ unsigned int const,
+ unsigned int const targetA) const {
+ if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) {
+ // i.e. nuclei handled by different process, this should not happen
+ throw std::runtime_error("CoM energy out of bounds for SIBYLL");
+ }
+
+ unsigned int targA = targetA;
+ if (is_nucleus(targetId) && targetId != Code::Nucleus) {
+ targA = get_nucleus_A(targetId);
+ }
+
+ if ((is_nucleus(targetId) &&
+ (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) &&
+ targetId != Code::Proton && targetId != Code::Hydrogen &&
+ targetId != Code::Neutron) {
+ throw std::runtime_error("Target cannot be handled by SIBYLL");
+ }
+
+ if (is_nucleus(projectileId) && !corsika::sibyll::canInteract(projectileId)) {
+ throw std::runtime_error("Projectile cannot be handled by SIBYLL");
+ }
+ }
+
+ inline std::tuple<CrossSectionType, CrossSectionType>
+ InteractionModel::getCrossSectionInelEla(Code const projectileId, Code const targetId,
+ HEPEnergyType const sqrtSnn,
+ unsigned int const projectileA,
+ unsigned int const targetA) const {
+
+ isValid(projectileId, targetId, sqrtSnn, projectileA, targetA); // throws
+
+ double dummy, dum1, dum3, dum4, dumdif[3]; // dummies needed for fortran call
+ int const iBeam = corsika::sibyll::getSibyllXSCode(
+ projectileId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
+
+ double const dEcm = sqrtSnn / 1_GeV;
+ // single nucleon target (p,n, hydrogen) or 4<=A<=18
+ double sigProd = 0;
+ double sigEla = 0;
+ // single nucleon target
+ if (targetId == Code::Proton || targetId == Code::Hydrogen ||
+ targetId == Code::Neutron) {
+ sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
+ } else {
+ // nuclear target
+ int const iTarget = targetA;
+ sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
+ }
+ return {sigProd * 1_mb, sigEla * 1_mb};
+ } // namespace corsika::sibyll
+
+ /**
+ * In this function SIBYLL is called to produce one event. The
+ * event is copied (and boosted) into the shower lab frame.
+ */
+
+ template <typename TSecondaryView>
+ inline void InteractionModel::doInteraction(
+ TSecondaryView& secondaries, COMBoost const& boost, Code const projectileId,
+ Code const targetId, HEPEnergyType const sqrtSnn, unsigned int const projectileA,
+ unsigned int const targetA) {
+
+ isValid(projectileId, targetId, sqrtSnn, projectileA, targetA); // throws
+
+ CORSIKA_LOG_DEBUG("pId={} tId={} sqrtSnn={}GeV", projectileId, targetId, sqrtSnn);
+
+ int targetSibCode = -1;
+ if (is_nucleus(targetId)) { targetSibCode = targetA; }
+ if (targetId == Proton::code) targetSibCode = 1;
+ CORSIKA_LOG_DEBUG("sibyll code: {}", targetSibCode);
+
+ // beam id for sibyll
+ int const projectileSibyllCode = corsika::sibyll::convertToSibyllRaw(projectileId);
+
+ count_++;
+ // Sibyll does not know about units..
+ double const sqs = sqrtSnn / 1_GeV;
+ // running sibyll, filling stack
+ sibyll_(projectileSibyllCode, targetSibCode, sqs);
+
+ if (sibyll_listing_) {
+ // print final state
+ int print_unit = 6;
+ sib_list_(print_unit);
+ nucCount_ += get_nwounded() - 1;
+ }
+
+ // ------ output and particle readout -----
+ auto const& csPrime = boost.getRotatedCS();
+
+ // add particles from sibyll to stack
+
+ // position and time of interaction, not used in Sibyll
+ Point const pOrig = Point(csPrime, {0_m, 0_m, 0_m});
+ TimeType const tOrig = 0_s; // no time in sibyll
+ CORSIKA_LOG_DEBUG("position of interaction: {}, time {} ", pOrig.getCoordinates(),
+ tOrig);
+
+ // link to sibyll stack
+ SibStack ss;
+
+ auto const& originalCS = boost.getOriginalCS();
+ MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
+ for (auto& psib : ss) {
+ // abort on particles that have decayed in Sibyll. Should not happen!
+ if (psib.hasDecayed())
+ throw std::runtime_error("found particle that decayed in SIBYLL!");
+
+ // transform 4-momentum to lab. frame
+ // note that the momentum needs to be rotated back
+ auto const tmp = psib.getMomentum().getComponents();
+ auto const pCoM = MomentumVector(csPrime, tmp);
+ HEPEnergyType const eCoM = psib.getEnergy();
+ auto const Plab = boost.fromCoM(FourVector{eCoM, pCoM});
+ auto const p3lab = Plab.getSpaceLikeComponents();
+
+ // add to corsika stack
+ auto pnew = secondaries.addSecondary(std::make_tuple(
+ corsika::sibyll::convertFromSibyll(psib.getPID()), p3lab, pOrig, tOrig));
+
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
+ Ecm_final += psib.getEnergy();
+ }
+ HEPEnergyType const Elab_initial =
+ static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
+ CORSIKA_LOG_DEBUG(
+ "conservation (all GeV): "
+ "sqrtSnn={}, sqrtSnn_final={}, "
+ "Elab_initial={}, Elab_final={}, "
+ "diff(%)={}, "
+ "E in nucleons={}, "
+ "Plab_final={} ",
+ sqrtSnn / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Elab_initial,
+ 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
diff --git a/corsika/detail/modules/sibyll/NuclearInteraction.inl b/corsika/detail/modules/sibyll/NuclearInteraction.inl
deleted file mode 100644
index 34d6ce22b..000000000
--- a/corsika/detail/modules/sibyll/NuclearInteraction.inl
+++ /dev/null
@@ -1,586 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/modules/sibyll/Interaction.hpp>
-#include <corsika/modules/sibyll/NuclearInteraction.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-#include <corsika/framework/core/Logging.hpp>
-
-#include <nuclib.hpp>
-
-namespace corsika::sibyll {
-
- template <typename TEnvironment>
- inline NuclearInteraction<TEnvironment>::NuclearInteraction(sibyll::Interaction& hadint,
- TEnvironment const& env)
- : environment_(env)
- , hadronicInteraction_(hadint) {
-
- // initialize hadronic interaction module
-
- // check compatibility of energy ranges, someone could try to use low-energy model..
- if (!hadronicInteraction_.isValidCoMEnergy(getMinEnergyPerNucleonCoM()) ||
- !hadronicInteraction_.isValidCoMEnergy(getMaxEnergyPerNucleonCoM())) {
- throw std::runtime_error(
- "NuclearInteraction: hadronic interaction model incompatible!");
- }
-
- // initialize nuclib
- // TODO: make sure this does not overlap with sibyll
- nuc_nuc_ini_();
-
- // initialize cross sections
- initializeNuclearCrossSections();
- }
-
- template <typename TEnvironment>
- inline NuclearInteraction<TEnvironment>::~NuclearInteraction() {
- CORSIKA_LOG_DEBUG("Nuclib::NuclearInteraction n={} Nnuc={}", count_, nucCount_);
- }
-
- template <typename TEnvironment>
- inline void NuclearInteraction<TEnvironment>::printCrossSectionTable(Code pCode) {
- if (pCode == Code::Argon) {
- CORSIKA_LOG_WARN("SIBYLL cannot handle Argon as target!");
- return;
- }
- const int k = targetComponentsIndex_.at(pCode);
- Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
- Code::Neon, Code::Argon, Code::Iron};
-
- std::ostringstream table;
- table << "Nuclear CrossSectionTable pCode=" << pCode << " :\n en/A ";
- for (auto& j : pNuclei) table << std::setw(9) << j;
- table << "\n";
-
- // loop over energy bins
- for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
- table << " " << i << " ";
-
- for (auto& n : pNuclei) {
- auto const j = get_nucleus_A(n);
- table << " " << std::setprecision(5) << std::setw(8)
- << cnucsignuc_.sigma[j - 1][k][i];
- }
- table << "\n";
- }
- CORSIKA_LOG_DEBUG(table.str());
- }
-
- template <typename TEnvironment>
- inline void NuclearInteraction<TEnvironment>::initializeNuclearCrossSections() {
-
- auto& universe = *(environment_.getUniverse());
-
- auto const allElementsInUniverse = std::invoke([&]() {
- std::set<Code> allElementsInUniverse;
- auto collectElements = [&](auto& vtn) {
- if (vtn.hasModelProperties()) {
- auto const& comp =
- vtn.getModelProperties().getNuclearComposition().getComponents();
- for (auto const c : comp) allElementsInUniverse.insert(c);
- }
- };
- universe.walk(collectElements);
- return allElementsInUniverse;
- });
-
- CORSIKA_LOG_DEBUG("NuclearInteraction: initializing nuclear cross sections...");
-
- // loop over target components, at most 4!!
- int k = -1;
- for (auto& ptarg : allElementsInUniverse) {
- if (ptarg == Code::Argon) continue; // NEED TO IGNORE Argon ....
- ++k;
- CORSIKA_LOG_DEBUG("NuclearInteraction: init target component: {}", ptarg);
- const int ib = get_nucleus_A(ptarg);
- if (!hadronicInteraction_.isValidTarget(ptarg)) {
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id={}",
- ptarg);
- throw std::runtime_error(
- " target can not be handled by hadronic interaction model! ");
- }
- targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
- // loop over energies, fNEnBins log. energy bins
- for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
- // hard coded energy grid, has to be aligned to definition in signuc2!!, no
- // comment..
- const HEPEnergyType Ecm = pow(10., 1. + 1. * i) * 1_GeV;
- // get p-p cross sections
- auto const protonId = Code::Proton;
- auto const [siginel, sigela] =
- hadronicInteraction_.getCrossSection(protonId, protonId, Ecm);
- const double dsig = siginel / 1_mb;
- const double dsigela = sigela / 1_mb;
- // loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
- for (unsigned int j = 1; j < gMaxNucleusAProjectile_; ++j) {
- const int jj = j + 1;
- double sig_out, dsig_out, sigqe_out, dsigqe_out;
- sigma_mc_(jj, ib, dsig, dsigela, gNSample_, sig_out, dsig_out, sigqe_out,
- dsigqe_out);
- // write to table
- cnucsignuc_.sigma[j][k][i] = sig_out;
- cnucsignuc_.sigqe[j][k][i] = sigqe_out;
- }
- }
- }
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: cross sections for {} "
- " components initialized!",
- targetComponentsIndex_.size());
- for (auto& ptarg : allElementsInUniverse) { printCrossSectionTable(ptarg); }
- }
-
- template <typename TEnvironment>
- inline CrossSectionType NuclearInteraction<TEnvironment>::readCrossSectionTable(
- const int ia, Code pTarget, HEPEnergyType elabnuc) {
-
- const int ib = targetComponentsIndex_.at(pTarget) + 1; // table index in fortran
- auto const ECoMNuc = sqrt(2. * constants::nucleonMass * elabnuc);
- if (ECoMNuc < getMinEnergyPerNucleonCoM() || ECoMNuc > getMaxEnergyPerNucleonCoM())
- throw std::runtime_error("NuclearInteraction: energy outside tabulated range!");
- const double e0 = elabnuc / 1_GeV;
- double sig;
- CORSIKA_LOG_DEBUG("ReadCrossSectionTable: {} {} {}", ia, ib, e0);
- signuc2_(ia, ib, e0, sig);
- CORSIKA_LOG_DEBUG("ReadCrossSectionTable: sig={}", sig);
- return sig * 1_mb;
- }
-
- // TODO: remove elastic cross section?
- template <typename TEnvironment>
- template <typename TParticle>
- std::tuple<CrossSectionType, CrossSectionType> inline NuclearInteraction<
- TEnvironment>::getCrossSection(TParticle const& projectile, Code const TargetId) {
-
- if (!is_nucleus(projectile.getPID())) {
- throw std::runtime_error(
- "NuclearInteraction: getCrossSection: particle not a nucleus!");
- }
-
- unsigned int const iBeamA = get_nucleus_A(projectile.getPID());
- HEPEnergyType LabEnergyPerNuc = projectile.getEnergy() / iBeamA;
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: getCrossSection: called with: beamNuclA={} "
- " TargetId={} LabEnergyPerNuc={}GeV ",
- iBeamA, TargetId, LabEnergyPerNuc / 1_GeV);
-
- // use nuclib to calc. nuclear cross sections
- // TODO: for now assumes air with hard coded composition
- // extend to arbitrary mixtures, requires smarter initialization
- // get nuclib projectile code: nucleon number
- if (iBeamA > getMaxNucleusAProjectile() || iBeamA < 2) {
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
- "A=" +
- std::to_string(iBeamA));
- throw std::runtime_error(
- "NuclearInteraction: getCrossSection: beam nucleus outside allowed range for "
- "NUCLIB!");
- }
-
- if (hadronicInteraction_.isValidTarget(TargetId)) {
- auto const sigProd = readCrossSectionTable(iBeamA, TargetId, LabEnergyPerNuc);
- CORSIKA_LOG_DEBUG("cross section (mb): " + std::to_string(sigProd / 1_mb));
- return std::make_tuple(sigProd, 0_mb);
- } else {
- throw std::runtime_error("target outside range.");
- }
- return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
- std::numeric_limits<double>::infinity() * 1_mb);
- }
-
- template <typename TEnvironment>
- template <typename TParticle>
- inline GrammageType NuclearInteraction<TEnvironment>::getInteractionLength(
- TParticle const& projectile) {
-
- // coordinate system, get global frame of reference
-
- const Code corsikaBeamId = projectile.getPID();
-
- if (!is_nucleus(corsikaBeamId)) {
- // no nuclear interaction
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
-
- // read from cross section code table
-
- MomentumVector pLab = projectile.getMomentum();
- CoordinateSystemPtr const& labCS = pLab.getCoordinateSystem();
-
- // FOR NOW: assume target is at rest
- MomentumVector pTarget(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
-
- // total momentum and energy
- HEPEnergyType Elab = projectile.getEnergy() + constants::nucleonMass;
- int const nuclA = get_nucleus_A(corsikaBeamId);
- auto const ElabNuc = projectile.getEnergy() / nuclA;
-
- MomentumVector pTotLab(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- pTotLab += pLab;
- pTotLab += pTarget;
- auto const pTotLabNorm = pTotLab.getNorm();
- // calculate cm. energy
- [[maybe_unused]] HEPEnergyType const ECoM = sqrt(
- (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
- auto const ECoMNN = sqrt(2. * ElabNuc * constants::nucleonMass);
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: LambdaInt: \n"
- " input energy: {}GeV\n"
- " input energy CoM: {}GeV\n"
- " beam pid: {}\n"
- " beam A: {}\n"
- " input energy per nucleon: {}GeV\n"
- " input energy CoM per nucleon: {}GeV ",
- Elab / 1_GeV, ECoM / 1_GeV, get_name(corsikaBeamId), nuclA, ElabNuc / 1_GeV,
- ECoMNN / 1_GeV);
-
- // energy limits
- // TODO: values depend on hadronic interaction model !! this is sibyll specific
- if (ElabNuc >= 8.5_GeV && ECoMNN >= gMinEnergyPerNucleonCoM_ &&
- ECoMNN < gMaxEnergyPerNucleonCoM_) {
-
- // get target from environment
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- auto const* const currentNode = projectile.getNode();
- auto const& mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
- // determine average interaction length
- // weighted sum
- int i = -1;
- CrossSectionType weightedProdCrossSection = 0_mb;
- // get weights of components from environment/medium
- const auto& w = mediumComposition.getFractions();
- // loop over components in medium
- for (auto const targetId : mediumComposition.getComponents()) {
- if (targetId == Code::Argon) continue; // NEED TO IGNORE Argon ....
- i++;
- CORSIKA_LOG_DEBUG("NuclearInteraction: get interaction length for target: {}",
- get_name(targetId));
- auto const [productionCrossSection, elaCrossSection] =
- getCrossSection(projectile, targetId);
- [[maybe_unused]] auto& dummy_elaCrossSection = elaCrossSection;
-
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: "
- "IntLength: nuclib return (mb): " +
- std::to_string(productionCrossSection / 1_mb));
- weightedProdCrossSection += w[i] * productionCrossSection;
- }
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: "
- "IntLength: weighted CrossSection (mb): {} ",
- weightedProdCrossSection / 1_mb);
-
- // calculate interaction length in medium
- GrammageType const int_length = mediumComposition.getAverageMassNumber() *
- constants::u / weightedProdCrossSection;
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: "
- "interaction length (g/cm2): {} ",
- int_length * (1_cm * 1_cm / (0.001_kg)));
-
- return int_length;
- } else {
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
- }
- }
-
- template <typename TEnvironment>
- template <typename TSecondaryView>
- inline void NuclearInteraction<TEnvironment>::doInteraction(TSecondaryView& view) {
-
- auto projectile = view.getProjectile();
-
- // this routine superimposes different nucleon-nucleon interactions
- // in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
-
- const auto ProjId = projectile.getPID();
- // TODO: calculate projectile mass in nuclearStackExtension
- // const auto ProjMass = projectile.getMass();
-
- CORSIKA_LOG_DEBUG("NuclearInteraction: DoInteraction: called with: {}",
- get_name(ProjId));
-
- // check if target-style nucleus (enum)
- if (!is_nucleus(ProjId)) {
- throw std::runtime_error(
- "NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
- "should use NuclearStackExtension!");
- }
-
- auto const ProjMass = get_mass(ProjId);
- CORSIKA_LOG_DEBUG("NuclearInteraction: projectile mass: {} ", ProjMass / 1_GeV);
-
- count_++;
-
- // position and time of interaction, not used in NUCLIB
- Point pOrig = projectile.getPosition();
- TimeType tOrig = projectile.getTime();
-
- CORSIKA_LOG_DEBUG("Interaction: position of interaction: {}", pOrig.getCoordinates());
- CORSIKA_LOG_DEBUG("Interaction: time: {} ", tOrig / 1_s);
-
- // projectile nucleon number
- const unsigned int kAProj = get_nucleus_A(ProjId);
- if (kAProj > getMaxNucleusAProjectile())
- throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
-
- // kinematics
- // define projectile nucleus
- HEPEnergyType const eProjectileLab = projectile.getEnergy();
- MomentumVector const pProjectileLab = projectile.getMomentum();
- FourVector const PprojLab(eProjectileLab, pProjectileLab);
- CoordinateSystemPtr const& labCS = pProjectileLab.getCoordinateSystem();
-
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: eProj lab: {} "
- "pProj lab: {} ",
- eProjectileLab / 1_GeV, pProjectileLab.getComponents() / 1_GeV);
-
- // define projectile nucleon
- HEPEnergyType const eProjectileNucLab = eProjectileLab / kAProj;
- MomentumVector const pProjectileNucLab = pProjectileLab / kAProj;
- FourVector const PprojNucLab(eProjectileNucLab, pProjectileNucLab);
-
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: eProjNucleon lab (GeV): {} "
- "pProjNucleon lab (GeV): {} ",
- eProjectileNucLab / 1_GeV, pProjectileNucLab.getComponents() / 1_GeV);
-
- // define target
- // always a nucleon
- // target is always at rest
- auto const eTargetNucLab = 0_GeV + constants::nucleonMass;
- auto const pTargetNucLab = MomentumVector(labCS, 0_GeV, 0_GeV, 0_GeV);
- FourVector const PtargNucLab(eTargetNucLab, pTargetNucLab);
-
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: etarget lab(GeV): {} "
- "NuclearInteraction: ptarget lab(GeV): {} ",
- eTargetNucLab / 1_GeV, pTargetNucLab.getComponents() / 1_GeV);
-
- // center-of-mass energy in nucleon-nucleon frame
- auto const PtotNN4 = PtargNucLab + PprojNucLab;
- HEPEnergyType EcmNN = PtotNN4.getNorm();
-
- CORSIKA_LOG_DEBUG("NuclearInteraction: nuc-nuc cm energy: {}", EcmNN / 1_GeV);
-
- if (!hadronicInteraction_.isValidCoMEnergy(EcmNN)) {
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
- "interaction model!");
- throw std::runtime_error("NuclearInteraction: DoInteraction: energy too low!");
- }
-
- // define boost to NUCLEON-NUCLEON frame
- COMBoost const boost(PprojNucLab, constants::nucleonMass);
- // boost projecticle
- auto const PprojNucCoM = boost.toCoM(PprojNucLab);
-
- // boost target
- auto const PtargNucCoM = boost.toCoM(PtargNucLab);
-
- CORSIKA_LOG_DEBUG(
- "Interaction: ebeam CoM: {} "
- ", pbeam CoM: {} ",
- PprojNucCoM.getTimeLikeComponent() / 1_GeV,
- PprojNucCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
- CORSIKA_LOG_DEBUG(
- "Interaction: etarget CoM: {}"
- ", ptarget CoM: {}",
- PtargNucCoM.getTimeLikeComponent() / 1_GeV,
- PtargNucCoM.getSpaceLikeComponents().getComponents() / 1_GeV);
-
- // sample target nucleon number
- //
- // proton stand-in for nucleon
- const auto beamId = Code::Proton;
- auto const* const currentNode = projectile.getNode();
- const auto& mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
- CORSIKA_LOG_DEBUG("get nucleon-nucleus cross sections for target materials..");
- // get cross sections for target materials
- // using nucleon-target-nucleus cross section!!!
- /*
- Here we read the cross section from the interaction model again,
- should be passed from getInteractionLength if possible
- */
- auto const& compVec = mediumComposition.getComponents();
- std::vector<CrossSectionType> cross_section_of_components(compVec.size());
-
- for (size_t i = 0; i < compVec.size(); ++i) {
- auto const targetId = compVec[i];
- if (targetId == Code::Argon) continue; // NEED TO IGNORE Argon ....
- CORSIKA_LOG_DEBUG("target component: {}", get_name(targetId));
- CORSIKA_LOG_DEBUG("beam id: {}", get_name(beamId));
- const auto [sigProd, sigEla] =
- hadronicInteraction_.getCrossSection(beamId, targetId, EcmNN);
- cross_section_of_components[i] = sigProd;
- [[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
- }
-
- const auto targetCode =
- mediumComposition.sampleTarget(cross_section_of_components, RNG_);
- CORSIKA_LOG_DEBUG("Interaction: target selected: {}", get_name(targetCode));
- /*
- FOR NOW: allow nuclei with A<18 or protons only.
- when medium composition becomes more complex, approximations will have to be
- allowed air in atmosphere also contains some Argon.
- */
- int kATarget = -1;
- if (is_nucleus(targetCode))
- kATarget = get_nucleus_A(targetCode);
- else if (targetCode == Code::Proton)
- kATarget = 1;
- CORSIKA_LOG_DEBUG("NuclearInteraction: nuclib target code: " +
- std::to_string(kATarget));
- if (!hadronicInteraction_.isValidTarget(targetCode))
- throw std::runtime_error("target outside range. ");
- // end of target sampling
-
- // superposition
- CORSIKA_LOG_DEBUG(
- "NuclearInteraction: sampling nuc. multiple interaction structure.. ");
- // get nucleon-nucleon cross section
- // (needed to determine number of nucleon-nucleon scatterings)
- const auto protonId = Code::Proton;
- const auto [prodCrossSection, elaCrossSection] =
- hadronicInteraction_.getCrossSection(protonId, protonId, EcmNN);
- const double sigProd = prodCrossSection / 1_mb;
- const double sigEla = elaCrossSection / 1_mb;
- // sample number of interactions (only input variables, output in common cnucms)
- // nuclear multiple scattering according to glauber (r.i.p.)
- int_nuc_(kATarget, kAProj, sigProd, sigEla);
-
- CORSIKA_LOG_DEBUG(
- "number of nucleons in target : {}\n"
- "number of wounded nucleons in target : {}\n"
- "number of nucleons in projectile : {}\n"
- "number of wounded nucleons in project. : {}\n"
- "number of inel. nuc.-nuc. interactions : {}\n"
- "number of elastic nucleons in target : {}\n"
- "number of elastic nucleons in project. : {}\n"
- "impact parameter: {}",
- kATarget, cnucms_.na, kAProj, cnucms_.nb, cnucms_.ni, cnucms_.nael, cnucms_.nbel,
- cnucms_.b);
-
- // calculate fragmentation
- CORSIKA_LOG_DEBUG("calculating nuclear fragments..");
- // number of interactions
- // include elastic
- const int nElasticNucleons = cnucms_.nbel;
- const int nInelNucleons = cnucms_.nb;
- const int nIntProj = nInelNucleons + nElasticNucleons;
- const double impactPar = cnucms_.b; // only needed to avoid passing common var.
- int nFragments = 0;
- // number of fragments is limited to 60
- int AFragments[60];
- // call fragmentation routine
- // input: target A, projectile A, number of int. nucleons in projectile, impact
- // parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
- // in pf in common fragments, neglected
- fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
-
- // this should not occur but well :) (LCOV_EXCL_START)
- if (nFragments > (int)getMaxNFragments())
- throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
- // (LCOV_EXCL_STOP)
-
- CORSIKA_LOG_DEBUG("number of fragments: " + std::to_string(nFragments));
- CORSIKA_LOG_DEBUG("adding nuclear fragments to particle stack..");
- // put nuclear fragments on corsika stack
- for (int j = 0; j < nFragments; ++j) {
- CORSIKA_LOG_DEBUG("fragment {}: A={} px={} py={} pz={}", j, AFragments[j],
- fragments_.ppp[j][0], fragments_.ppp[j][1], fragments_.ppp[j][2]);
- const auto nuclA = AFragments[j];
- // get Z from stability line
- const auto nuclZ = int(nuclA / 2.15 + 0.7);
-
- // TODO: do we need to catch single nucleons??
- Code specCode = Code::Neutron; // sample neutron or proton ?
- if (nuclA > 1) specCode = get_nucleus_code(nuclA, nuclZ);
- HEPMassType const mass = get_mass(specCode);
-
- CORSIKA_LOG_DEBUG("NuclearInteraction: adding fragment: {}", get_name(specCode));
- CORSIKA_LOG_DEBUG("NuclearInteraction: A,Z: {}, {}", nuclA, nuclZ);
- CORSIKA_LOG_DEBUG("NuclearInteraction: mass: {} GeV", std::to_string(mass / 1_GeV));
-
- // CORSIKA 7 way
- // spectators inherit momentum from original projectile
- const double mass_ratio = mass / ProjMass;
-
- CORSIKA_LOG_DEBUG("NuclearInteraction: mass ratio " + std::to_string(mass_ratio));
-
- auto const Plab = PprojLab * mass_ratio;
-
- CORSIKA_LOG_DEBUG("NuclearInteraction: fragment momentum: {}",
- Plab.getSpaceLikeComponents().getComponents() / 1_GeV);
-
- projectile.addSecondary(
- std::make_tuple(specCode, Plab.getSpaceLikeComponents(), pOrig, tOrig));
- }
-
- // add elastic nucleons to corsika stack
- // TODO: the elastic interaction could be external like the inelastic interaction,
- // e.g. use existing ElasticModel
- CORSIKA_LOG_DEBUG("adding elastically scattered nucleons to particle stack..");
- for (int j = 0; j < nElasticNucleons; ++j) {
- // TODO: sample proton or neutron
- auto const elaNucCode = Code::Proton;
-
- // CORSIKA 7 way
- // elastic nucleons inherit momentum from original projectile
- // neglecting momentum transfer in interaction
- const double mass_ratio = get_mass(elaNucCode) / ProjMass;
- auto const Plab = PprojLab * mass_ratio;
-
- projectile.addSecondary(
- std::make_tuple(elaNucCode, Plab.getSpaceLikeComponents(), pOrig, tOrig));
- }
-
- // add inelastic interactions
- CORSIKA_LOG_DEBUG("calculate inelastic nucleon-nucleon interactions..");
- for (int j = 0; j < nInelNucleons; ++j) {
- // TODO: sample neutron or proton
- auto 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, PprojNucLab.getSpaceLikeComponents(), pOrig, tOrig));
- inelasticNucleon.setNode(projectile.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);
- for (const auto& pSec : nucleon_secondaries) {
- projectile.addSecondary(std::make_tuple(pSec.getPID(), pSec.getMomentum(),
- pSec.getPosition(), pSec.getTime()));
- }
- }
-
- CORSIKA_LOG_DEBUG("NuclearInteraction: DoInteraction: done");
- }
-
-} // namespace corsika::sibyll
diff --git a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
new file mode 100644
index 000000000..8811e0eea
--- /dev/null
+++ b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
@@ -0,0 +1,363 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <nuclib.hpp>
+
+namespace corsika::sibyll {
+
+ template <typename TEnvironment, typename TNucleonModel>
+ inline NuclearInteractionModel<TEnvironment, TNucleonModel>::NuclearInteractionModel(
+ TNucleonModel& hadint, TEnvironment const& env)
+ : environment_(env)
+ , hadronicInteraction_(hadint) {
+
+ // initialize nuclib
+ // TODO: make sure this does not overlap with sibyll
+ nuc_nuc_ini_();
+
+ // initialize cross sections
+ initializeNuclearCrossSections();
+ }
+
+ template <typename TEnvironment, typename TNucleonModel>
+ inline NuclearInteractionModel<TEnvironment,
+ TNucleonModel>::~NuclearInteractionModel() {
+ CORSIKA_LOG_DEBUG("Nuclib::NuclearInteractionModel n={} Nnuc={}", count_, nucCount_);
+ }
+
+ template <typename TEnvironment, typename TNucleonModel>
+ inline void constexpr NuclearInteractionModel<TEnvironment, TNucleonModel>::isValid(
+ Code const projectileId, Code const targetId, HEPEnergyType const sqrtSnn,
+ unsigned int const projectileA, unsigned int const targetA) const {
+
+ // also depends on underlying model, for Proton/Neutron projectile
+ hadronicInteraction_.isValid(Code::Proton, targetId, sqrtSnn, 1, targetA); // throws
+
+ // projectile limits:
+ if (is_nucleus(projectileId)) {
+ throw std::runtime_error("can only handle nuclear projectile");
+ }
+ if (projectileA >= getMaxNucleusAProjectile() || projectileA < 2) {
+ throw std::runtime_error("projectile mass A out of bounds");
+ }
+ } // namespace corsika::sibyll
+
+ template <typename TEnvironment, typename TNucleonModel>
+ inline void
+ NuclearInteractionModel<TEnvironment, TNucleonModel>::printCrossSectionTable(
+ Code const pCode) const {
+ if (pCode == Code::Argon) {
+ CORSIKA_LOG_WARN("SIBYLL cannot handle Argon as target!");
+ return;
+ }
+ int const k = targetComponentsIndex_.at(pCode);
+ Code const pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
+ Code::Neon, Code::Argon, Code::Iron};
+
+ std::ostringstream table;
+ table << "Nuclear CrossSectionTable pCode=" << pCode << " :\n en/A ";
+ for (auto& j : pNuclei) table << std::setw(9) << j;
+ table << "\n";
+
+ // loop over energy bins
+ for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
+ table << " " << i << " ";
+
+ for (auto& n : pNuclei) {
+ auto const j = get_nucleus_A(n);
+ table << " " << std::setprecision(5) << std::setw(8)
+ << cnucsignuc_.sigma[j - 1][k][i];
+ }
+ table << "\n";
+ }
+ CORSIKA_LOG_DEBUG(table.str());
+ }
+
+ template <typename TEnvironment, typename TNucleonModel>
+ inline void
+ NuclearInteractionModel<TEnvironment, TNucleonModel>::initializeNuclearCrossSections() {
+
+ auto& universe = *(environment_.getUniverse());
+ // generate complete list of all nuclei types in universe
+
+ auto const allElementsInUniverse = std::invoke([&]() {
+ std::set<Code> allElementsInUniverse;
+ auto collectElements = [&](auto& vtn) {
+ if (vtn.hasModelProperties()) {
+ auto const& comp =
+ vtn.getModelProperties().getNuclearComposition().getComponents();
+ for (auto const c : comp) allElementsInUniverse.insert(c);
+ }
+ };
+ universe.walk(collectElements);
+ return allElementsInUniverse;
+ });
+
+ CORSIKA_LOG_DEBUG("initializing nuclear cross sections...");
+
+ // loop over target components, at most 4!!
+ int k = -1;
+ for (Code const ptarg : allElementsInUniverse) {
+ 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"
+ hadronicInteraction_.isValid(Code::Proton, ptarg, 100_GeV, 1, ib); // throws
+ targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
+ // loop over energies, fNEnBins log. energy bins
+ for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
+ // hard coded energy grid, has to be aligned to definition in signuc2!!, no
+ // comment..
+ HEPEnergyType const Ecm = pow(10., 1. + 1. * i) * 1_GeV;
+ // get p-p cross sections
+ auto const protonId = Code::Proton;
+ auto const [siginel, sigela] =
+ hadronicInteraction_.getCrossSectionInelEla(protonId, protonId, Ecm);
+ const double dsig = siginel / 1_mb;
+ const double dsigela = sigela / 1_mb;
+ // loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
+ for (unsigned int j = 1; j < gMaxNucleusAProjectile_; ++j) {
+ const int jj = j + 1;
+ double sig_out, dsig_out, sigqe_out, dsigqe_out;
+ sigma_mc_(jj, ib, dsig, dsigela, gNSample_, sig_out, dsig_out, sigqe_out,
+ dsigqe_out);
+ // write to table
+ cnucsignuc_.sigma[j][k][i] = sig_out;
+ cnucsignuc_.sigqe[j][k][i] = sigqe_out;
+ }
+ }
+ }
+ CORSIKA_LOG_DEBUG("cross sections for {} components initialized!",
+ targetComponentsIndex_.size());
+ for (auto& ptarg : allElementsInUniverse) { printCrossSectionTable(ptarg); }
+ }
+
+ template <typename TEnvironment, typename TNucleonModel>
+ inline CrossSectionType
+ NuclearInteractionModel<TEnvironment, TNucleonModel>::readCrossSectionTable(
+ int const ia, Code const pTarget, HEPEnergyType const elabnuc) const {
+
+ int const ib = targetComponentsIndex_.at(pTarget) + 1; // table index in fortran
+ auto const ECoMNuc = sqrt(2. * constants::nucleonMass * elabnuc);
+ if (ECoMNuc < getMinEnergyPerNucleonCoM() || ECoMNuc > getMaxEnergyPerNucleonCoM()) {
+ throw std::runtime_error("energy outside tabulated range!");
+ }
+ double const e0 = elabnuc / 1_GeV;
+ double sig;
+ CORSIKA_LOG_DEBUG("ReadCrossSectionTable: {} {} {}", ia, ib, e0);
+ signuc2_(ia, ib, e0, sig);
+ CORSIKA_LOG_DEBUG("ReadCrossSectionTable: sig={}", sig);
+ return sig * 1_mb;
+ }
+
+ // TODO: remove elastic cross section?
+ template <typename TEnvironment, typename TNucleonModel>
+ CrossSectionType inline NuclearInteractionModel<
+ TEnvironment, TNucleonModel>::getCrossSection(Code const projectileId,
+ Code const targetId,
+ HEPEnergyType const sqrtSnn,
+ unsigned int const projectileA,
+ unsigned int const targetA) const {
+
+ isValid(projectileId, targetId, sqrtSnn, projectileA, targetA); // throws
+ HEPEnergyType const LabEnergyPerNuc =
+ static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
+ auto const sigProd = readCrossSectionTable(projectileA, targetId, LabEnergyPerNuc);
+ CORSIKA_LOG_DEBUG("cross section (mb): {}", sigProd / 1_mb);
+ return sigProd;
+ }
+
+ 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, unsigned int const projectileA,
+ unsigned int const targetA) {
+
+ isValid(projectileId, targetId, sqrtSnn, projectileA, targetA); // throws
+
+ CORSIKA_LOG_DEBUG("pId={} tId={} sqrtSnn={}GeV", projectileId, targetId,
+ sqrtSnn / 1_GeV);
+ count_++;
+
+ HEPEnergyType const ProjMass = projectileA * constants::nucleonMass;
+ // is this really more precise?:
+ // projectile.getNuclearZ() * Proton::mass +
+ // (projectile.getNuclearA() - projectile.getNuclearZ()) * Neutron::mass;
+
+ // 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});
+
+ /*
+ FOR NOW: allow nuclei with A<18 or protons only.
+ when medium composition becomes more complex, approximations will have to be
+ allowed air in atmosphere also contains some Argon.
+ */
+ int kATarget = -1;
+ if (is_nucleus(targetId)) {
+ kATarget = targetA;
+ } else if (targetId == Code::Proton) {
+ kATarget = 1;
+ }
+ CORSIKA_LOG_DEBUG("nuclib target code: {}", kATarget);
+
+ // end of target sampling
+
+ // superposition
+ CORSIKA_LOG_DEBUG("sampling nuc. multiple interaction structure.. ");
+ // get nucleon-nucleon cross section
+ // (needed to determine number of nucleon-nucleon scatterings)
+ auto const protonId = Code::Proton;
+ auto const [prodCrossSection, elaCrossSection] =
+ hadronicInteraction_.getCrossSectionInelEla(protonId, protonId, sqrtSnn);
+ double const sigProd = prodCrossSection / 1_mb;
+ double const sigEla = elaCrossSection / 1_mb;
+ // sample number of interactions (only input variables, output in common cnucms)
+ // nuclear multiple scattering according to glauber (r.i.p.)
+ int_nuc_(kATarget, projectileA, sigProd, sigEla);
+
+ CORSIKA_LOG_DEBUG(
+ "number of nucleons in target : {}\n"
+ "number of wounded nucleons in target : {}\n"
+ "number of nucleons in projectile : {}\n"
+ "number of wounded nucleons in project. : {}\n"
+ "number of inel. nuc.-nuc. interactions : {}\n"
+ "number of elastic nucleons in target : {}\n"
+ "number of elastic nucleons in project. : {}\n"
+ "impact parameter: {}",
+ kATarget, cnucms_.na, projectileA, cnucms_.nb, cnucms_.ni, cnucms_.nael,
+ cnucms_.nbel, cnucms_.b);
+
+ // calculate fragmentation
+ CORSIKA_LOG_DEBUG("calculating nuclear fragments..");
+ // number of interactions
+ // include elastic
+ int const nElasticNucleons = cnucms_.nbel;
+ int const nInelNucleons = cnucms_.nb;
+ int const nIntProj = nInelNucleons + nElasticNucleons;
+ double const impactPar = cnucms_.b; // only needed to avoid passing common var.
+ int nFragments = 0;
+ // number of fragments is limited to 60
+ int AFragments[60];
+ // call fragmentation routine
+ // input: target A, projectile A, number of int. nucleons in projectile, impact
+ // parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
+ // in pf in common fragments, neglected
+ fragm_(kATarget, projectileA, nIntProj, impactPar, nFragments, AFragments);
+
+ // this should not occur but well :) (LCOV_EXCL_START)
+ if (nFragments > (int)getMaxNFragments()) {
+ throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
+ }
+ // (LCOV_EXCL_STOP)
+
+ // position and time of interaction, not used in NUCLIB
+ Point pOrig{boost.getOriginalCS(), {0_m, 0_m, 0_m}}; // = projectile.getPosition();
+ TimeType delay = 0_s; // there is no time in sibyll
+
+ CORSIKA_LOG_DEBUG("Interaction: position of interaction: {} {}",
+ pOrig.getCoordinates(), delay / 1_s);
+ CORSIKA_LOG_DEBUG("number of fragments: {}", nFragments);
+ CORSIKA_LOG_DEBUG("adding nuclear fragments to particle stack..");
+ // put nuclear fragments on corsika stack
+ for (int j = 0; j < nFragments; ++j) {
+ CORSIKA_LOG_DEBUG("fragment {}: A={} px={} py={} pz={}", j, AFragments[j],
+ fragments_.ppp[j][0], fragments_.ppp[j][1], fragments_.ppp[j][2]);
+ Code specCode;
+ auto const nuclA = AFragments[j];
+ // get Z from stability line
+ auto const nuclZ = int(nuclA / 2.15 + 0.7);
+
+ // TODO: do we need to catch single nucleons??
+ if (nuclA == 1) {
+ // TODO: sample neutron or proton
+ specCode = Code::Proton;
+ } else {
+ specCode = Code::Nucleus;
+ }
+ HEPMassType const mass = get_nucleus_mass(nuclA, nuclZ);
+
+ CORSIKA_LOG_DEBUG("adding fragment: {}", get_name(specCode));
+ CORSIKA_LOG_DEBUG("A,Z: {}, {}", nuclA, nuclZ);
+ CORSIKA_LOG_DEBUG("mass: {} GeV", mass / 1_GeV);
+
+ // 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);
+
+ if (nuclA == 1) {
+ // add nucleon
+ view.addSecondary(std::make_tuple(specCode, p3lab, pOrig, delay));
+ } else {
+ // add nucleus
+ view.addSecondary(std::make_tuple(specCode, p3lab, pOrig, delay, nuclA, nuclZ));
+ }
+ }
+
+ // add elastic nucleons to corsika stack
+ // TODO: the elastic interaction could be external like the inelastic interaction,
+ // e.g. use existing ElasticModel
+ CORSIKA_LOG_DEBUG("adding elastically scattered nucleons to particle stack..");
+ for (int j = 0; j < nElasticNucleons; ++j) {
+ // TODO: sample proton or neutron
+ Code const elaNucCode = Code::Proton;
+
+ // 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;
+ view.addSecondary(std::make_tuple(elaNucCode, p3lab, pOrig, delay));
+ }
+
+ // add inelastic interactions
+ CORSIKA_LOG_DEBUG("calculate inelastic nucleon-nucleon interactions..");
+ for (int j = 0; j < nInelNucleons; ++j) {
+ // TODO: sample neutron or proton
+ auto 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));
+ 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, 0, targetA);
+ for (const auto& pSec : nucleon_secondaries) {
+ view.addSecondary(std::make_tuple(pSec.getPID(), pSec.getMomentum(),
+ pSec.getPosition(), pSec.getTime()));
+ }
+ }
+ }
+
+} // namespace corsika::sibyll
diff --git a/corsika/detail/stack/VectorStack.inl b/corsika/detail/stack/VectorStack.inl
index 9d2a9434c..57bc094a4 100644
--- a/corsika/detail/stack/VectorStack.inl
+++ b/corsika/detail/stack/VectorStack.inl
@@ -41,7 +41,7 @@ namespace corsika {
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
- ParticleInterface<StackIteratorInterface> const&,
+ ParticleInterface<StackIteratorInterface> const& parent,
particle_data_momentum_type const& v) {
this->setPID(std::get<0>(v));
MomentumVector const p = std::get<1>(v);
@@ -54,7 +54,7 @@ namespace corsika {
this->setDirection(p / sqrt(P2));
}
this->setPosition(std::get<2>(v));
- this->setTime(std::get<3>(v));
+ this->setTime(std::get<3>(v) + parent.getTime()); // parent time is added
}
template <typename StackIteratorInterface>
@@ -69,12 +69,13 @@ namespace corsika {
template <typename StackIteratorInterface>
inline void ParticleInterface<StackIteratorInterface>::setParticleData(
- ParticleInterface<StackIteratorInterface> const&, particle_data_type const& v) {
+ ParticleInterface<StackIteratorInterface> const& parent,
+ particle_data_type const& v) {
this->setPID(std::get<0>(v));
this->setKineticEnergy(std::get<1>(v));
this->setDirection(std::get<2>(v));
this->setPosition(std::get<3>(v));
- this->setTime(std::get<4>(v));
+ this->setTime(std::get<4>(v) + parent.getTime()); // parent time is added
}
template <typename StackIteratorInterface>
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index 1163950f3..796655a03 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -33,6 +33,8 @@
namespace corsika {
+ class COMBoost; // fwd-decl
+
/**
*
* The Cascade class is constructed from template arguments making
@@ -127,8 +129,10 @@ namespace corsika {
void step(particle_type& vParticle);
ProcessReturn decay(stack_view_type& view, InverseTimeType initial_inv_decay_time);
- ProcessReturn interaction(stack_view_type& view,
- InverseGrammageType initial_inv_int_length);
+ ProcessReturn interaction(stack_view_type& view, COMBoost const& boost,
+ HEPEnergyType const sqrtSnn,
+ NuclearComposition const& composition,
+ CrossSectionType const initial_cross_section);
void setEventType(stack_view_type& view, history::EventType);
// data members
diff --git a/corsika/framework/process/BaseProcess.hpp b/corsika/framework/process/BaseProcess.hpp
index 23e8667b2..a871a065e 100644
--- a/corsika/framework/process/BaseProcess.hpp
+++ b/corsika/framework/process/BaseProcess.hpp
@@ -41,8 +41,8 @@ namespace corsika {
/** @name getRef Return reference to underlying type
@{
*/
- TDerived& ref() { return static_cast<TDerived&>(*this); }
- const TDerived& ref() const { return static_cast<const TDerived&>(*this); }
+ TDerived& getRef() { return static_cast<TDerived&>(*this); }
+ const TDerived& getRef() const { return static_cast<const TDerived&>(*this); }
//! @}
public:
diff --git a/corsika/framework/process/DecayProcess.hpp b/corsika/framework/process/DecayProcess.hpp
index 6fa5a754c..a10434392 100644
--- a/corsika/framework/process/DecayProcess.hpp
+++ b/corsika/framework/process/DecayProcess.hpp
@@ -50,7 +50,7 @@ namespace corsika {
template <typename TDerived>
struct DecayProcess : BaseProcess<TDerived> {
public:
- using BaseProcess<TDerived>::ref;
+ using BaseProcess<TDerived>::getRef;
template <typename TParticle>
InverseTimeType getInverseLifetime(TParticle const& particle) {
@@ -61,7 +61,7 @@ namespace corsika {
"getInteractionLength(TParticle const&)\" required for "
"InteractionProcess<TDerived>. ");
- return 1. / ref().getLifetime(particle);
+ return 1. / getRef().getLifetime(particle);
}
};
diff --git a/corsika/framework/process/InteractionProcess.hpp b/corsika/framework/process/InteractionProcess.hpp
index c446e6f90..469ebbf89 100644
--- a/corsika/framework/process/InteractionProcess.hpp
+++ b/corsika/framework/process/InteractionProcess.hpp
@@ -10,67 +10,54 @@
#include <corsika/framework/process/BaseProcess.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/media/NuclearComposition.hpp>
#include <corsika/detail/framework/process/InteractionProcess.hpp> // for extra traits, method/interface checking
namespace corsika {
/**
- @ingroup Processes
- @{
-
- Process describing the interaction of particles
-
- Create a new InteractionProcess, e.g. for XYModel, via
- @code
- class XYModel : public InteractionProcess<XYModel> {};
- @endcode
-
- and provide the two necessary interface methods
- @code
- template <typename TSecondaryView>
- void XYModel::doInteraction(TSecondaryView&);
-
- template <typename TParticle>
- GrammageType XYModel::getInteractionLength(TParticle const&)
- @endcode
-
- Where, of course, SecondaryView and Particle are the valid
- classes to access particles on the Stack. In user code, those two methods do
- not need to be templated, they could use the types
- e.g. corsika::setup::Stack::particle_type -- but by the cost of
- loosing all flexibility otherwise provided.
-
- SecondaryView allows to retrieve the properties of the projectile
- particles, AND to store new particles (secondaries) which then
- subsequently can be processes by SecondariesProcess. This is how
- the output of interactions can be studied right away.
-
+ * @ingroup Processes
+ * @{
+ *
+ * Process describing the interaction of particles
+ *
+ * Create a new InteractionProcess, e.g. for XYModel, via
+ * @code
+ * class XYModel : public InteractionProcess<XYModel> {};
+ * @endcode
+ *
+ * and provide the two necessary interface methods
+ * @code
+ * template <typename TSecondaryView>
+ * void XYModel::doInteraction(TSecondaryView&);
+ *
+ * template <typename TParticle>
+ * GrammageType XYModel::getInteractionLength(TParticle const&)
+ * @endcode
+ *
+ * Where, of course, SecondaryView and Particle are the valid
+ * classes to access particles on the Stack. In user code, those two methods do
+ * not need to be templated, they could use the types
+ * e.g. corsika::setup::Stack::particle_type -- but by the cost of
+ * loosing all flexibility otherwise provided.
+ *
+ * SecondaryView allows to retrieve the properties of the projectile
+ * particles, AND to store new particles (secondaries) which then
+ * subsequently can be processes by SecondariesProcess. This is how
+ * the output of interactions can be studied right away.
*/
- template <typename TDerived>
- class InteractionProcess : public BaseProcess<TDerived> {
+ template <typename TModel>
+ class InteractionProcess : public BaseProcess<TModel> {
public:
- using BaseProcess<TDerived>::ref;
-
- template <typename TParticle>
- InverseGrammageType getInverseInteractionLength(TParticle const& particle) {
-
- // interface checking on TProcess1
- static_assert(
- has_method_getInteractionLength_v<TDerived, GrammageType, TParticle const&>,
- "TDerived has no method with correct signature \"GrammageType "
- "getInteractionLength(TParticle const&)\" required for "
- "InteractionProcess<TDerived>. ");
-
- return 1. / ref().getInteractionLength(particle);
- }
+ using BaseProcess<TModel>::getRef;
};
/**
- * ProcessTraits specialization to flag InteractionProcess objects
- **/
+ * ProcessTraits specialization to flag InteractionProcess objects.
+ */
template <typename TProcess>
struct is_interaction_process<
TProcess, std::enable_if_t<
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 407f40349..34407d8c4 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -25,14 +25,19 @@
#include <corsika/framework/process/SecondariesProcess.hpp>
#include <corsika/framework/process/StackProcess.hpp>
#include <corsika/framework/process/NullModel.hpp>
+
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
namespace corsika {
+ class COMBoost; // fwd-decl
+ class NuclearComposition; // fwd-decl
+
/**
- count_processes traits specialization to increase process count by
- getNumberOfProcesses(). This is used to statically count processes in the sequence
- **/
+ * count_processes traits specialization to increase process count by
+ * getNumberOfProcesses(). This is used to statically count processes in the sequence.
+ */
template <typename TProcess, int N>
struct count_processes<
TProcess, N,
@@ -43,110 +48,111 @@ namespace corsika {
};
/**
- @defgroup Processes Physics Processes and Modules
-
- Physics processes in CORSIKA 8 are clustered in ProcessSequence and
- SwitchProcessSequence containers. The former is a mere (ordered) collection, while
- the latter has the option to switch between two alternative ProcessSequences.
-
- Depending on the type of data to act on and on the allowed actions of processes there
- are several interface options:
- - InteractionProcess
- - DecayProcess
- - ContinuousProcess
- - StackProcess
- - SecondariesProcess
- - BoundaryCrossingProcess
-
- And all processes (including ProcessSequence and SwitchProcessSequence) are derived
- from BaseProcess.
-
- Processes of any type (e.g. p1, p2, p3,...) can be assembled into a ProcessSequence
- using the `make_sequence` factory function.
-
- @code{.cpp}
- auto sequence1 = make_sequence(p1, p2, p3);
- auto sequence2 = make_sequence(p4, p5, p6, p7);
- auto sequence3 = make_sequence(sequence1, sequemce2, p8, p9);
- @endcode
-
- Note, if the order of processes
- matters, the order of occurence
- in the ProcessSequence determines
- the executiion order.
-
- SecondariesProcess alyways act on
- new secondaries produced (i.e. in
- InteractionProcess and
- DecayProcess) in the scope of
- their ProcessSequence. For
- example if i1 and i2 are
- InteractionProcesses and s1 is a
- SecondariesProcess, then
-
- @code{.cpp}
- auto sequence = make_sequence(i1, make_sequence(i2, s1))
- @endcode
-
- will result in s1 acting only on
- the particles produced by i2 and
- not by i1. This can be very
- useful, e.g. to fine tune thinning.
-
- A special type of ProcessSequence
- is SwitchProcessSequence, which
- has two branches and a functor
- that can select between these two
- branches.
-
- @code{.cpp}
- auto sequence = make_switch(sequence1, sequence2, selector);
- @endcode
-
- where the only requirement to
- `selector` is that it
- provides a `SwitchResult operator()(Particle const& particle) const` method. Thus,
- based on the dynamic properties
- of `particle` the functor
- can make its decision. This is
- clearly important for switching
- between low-energy and
- high-energy models, but not
- limited to this. The selection
- can even be done with a lambda
- function.
-
-
- @class ProcessSequence
- @ingroup Processes
-
- Definition of a static process list/sequence
-
- A compile time static list of processes. The compiler will
- generate a new type based on template logic containing all the
- elements provided by the user.
-
- TProcess1 and TProcess2 must both be derived from BaseProcess,
- and are both references if possible (lvalue), otherwise (rvalue)
- they are just classes. This allows us to handle both, rvalue as
- well as lvalue Processes in the ProcessSequence.
-
- (For your potential interest,
- the static version of the
- ProcessSequence and all Process
- types are based on the CRTP C++
- design pattern)
-
- Template parameters:
- @tparam TProcess1 is of type BaseProcess, either a dedicatd process, or a
- ProcessSequence
- @tparam TProcess2 is of type BaseProcess, either a dedicatd process, or a
- ProcessSequence
- @tparam IndexFirstProcess to count and index each Process in the entire
- process-chain. The offset is the starting value for this ProcessSequence
- @tparam IndexOfProcess1 index of TProcess1 (counting of Process)
- @tparam IndexOfProcess2 index of TProcess2 (counting of Process)
- **/
+ * @defgroup Processes Physics Processes and Modules
+ *
+ * Physics processes in CORSIKA 8 are clustered in ProcessSequence and
+ * SwitchProcessSequence containers. The former is a mere (ordered) collection, while
+ * the latter has the option to switch between two alternative ProcessSequences.
+ *
+ * Depending on the type of data to act on and on the allowed actions of processes
+ there
+ * are several interface options:
+ * - InteractionProcess
+ * - DecayProcess
+ * - ContinuousProcess
+ * - StackProcess
+ * - SecondariesProcess
+ * - BoundaryCrossingProcess
+ *
+ * And all processes (including ProcessSequence and SwitchProcessSequence) are derived
+ * from BaseProcess.
+ *
+ * Processes of any type (e.g. p1, p2, p3,...) can be assembled into a ProcessSequence
+ * using the `make_sequence` factory function.
+ *
+ * @code{.cpp}
+ * auto sequence1 = make_sequence(p1, p2, p3);
+ * auto sequence2 = make_sequence(p4, p5, p6, p7);
+ * auto sequence3 = make_sequence(sequence1, sequemce2, p8, p9);
+ * @endcode
+ *
+ * Note, if the order of processes
+ * matters, the order of occurence
+ * in the ProcessSequence determines
+ * the executiion order.
+ *
+ * SecondariesProcess alyways act on
+ * new secondaries produced (i.e. in
+ * InteractionProcess and
+ * DecayProcess) in the scope of
+ * their ProcessSequence. For
+ * example if i1 and i2 are
+ * InteractionProcesses and s1 is a
+ * SecondariesProcess, then
+ *
+ * @code{.cpp}
+ * auto sequence = make_sequence(i1, make_sequence(i2, s1))
+ * @endcode
+ *
+ * will result in s1 acting only on
+ * the particles produced by i2 and
+ * not by i1. This can be very
+ * useful, e.g. to fine tune thinning.
+ *
+ * A special type of ProcessSequence
+ * is SwitchProcessSequence, which
+ * has two branches and a functor
+ * that can select between these two
+ * branches.
+ *
+ * @code{.cpp}
+ * auto sequence = make_switch(sequence1, sequence2, selector);
+ * @endcode
+ *
+ * where the only requirement to
+ * `selector` is that it
+ * provides a `SwitchResult operator()(Particle const& particle) const` method. Thus,
+ * based on the dynamic properties
+ * of `particle` the functor
+ * can make its decision. This is
+ * clearly important for switching
+ * between low-energy and
+ * high-energy models, but not
+ * limited to this. The selection
+ * can even be done with a lambda
+ * function.
+ *
+ *
+ * @class ProcessSequence
+ * @ingroup Processes
+ *
+ * Definition of a static process list/sequence
+ *
+ * A compile time static list of processes. The compiler will
+ * generate a new type based on template logic containing all the
+ * elements provided by the user.
+ *
+ * TProcess1 and TProcess2 must both be derived from BaseProcess,
+ * and are both references if possible (lvalue), otherwise (rvalue)
+ * they are just classes. This allows us to handle both, rvalue as
+ * well as lvalue Processes in the ProcessSequence.
+ *
+ * (For your potential interest,
+ * the static version of the
+ * ProcessSequence and all Process
+ * types are based on the CRTP C++
+ * design pattern)
+ *
+ * Template parameters:
+ * @tparam TProcess1 is of type BaseProcess, either a dedicatd process, or a
+ * ProcessSequence.
+ * @tparam TProcess2 is of type BaseProcess, either a dedicatd process, or a
+ * ProcessSequence.
+ * @tparam IndexFirstProcess to count and index each Process in the entire
+ * process-chain. The offset is the starting value for this ProcessSequence.
+ * @tparam IndexOfProcess1 index of TProcess1 (counting of Process).
+ * @tparam IndexOfProcess2 index of TProcess2 (counting of Process).
+ */
template <typename TProcess1, typename TProcess2 = NullModel,
int ProcessIndexOffset = 0,
@@ -171,15 +177,15 @@ namespace corsika {
static bool const is_process_sequence = true;
/**
- 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
- construct sequences.
-
- @param in_A BaseProcess or switch/process list
- @param in_B BaseProcess or switch/process list
- **/
+ * 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
+ * construct sequences.
+ *
+ * @param in_A BaseProcess or switch/process list.
+ * @param in_B BaseProcess or switch/process list.
+ */
ProcessSequence(TProcess1 in_A, TProcess2 in_B);
template <typename TParticle>
@@ -195,17 +201,17 @@ namespace corsika {
void doSecondaries(TSecondaries& vS);
/**
- The processes of type StackProcess do have an internal counter,
- so they can be exectuted only each N steps. Often these are
- "maintenacne processes" that do not need to run after each
- single step of the simulations. In the CheckStep function it is
- tested if either A_ or B_ are StackProcess and if they are due
- for execution.
+ * The processes of type StackProcess do have an internal counter,
+ * so they can be exectuted only each N steps. Often these are
+ * "maintenacne processes" that do not need to run after each
+ * single step of the simulations. In the CheckStep function it is
+ * tested if either A_ or B_ are StackProcess and if they are due
+ * for execution.
*/
bool checkStep();
/**
- Execute the StackProcess-es in the ProcessSequence
+ * Execute the StackProcess-es in the ProcessSequence.
*/
template <typename TStack>
void doStack(TStack& stack);
@@ -223,49 +229,45 @@ namespace corsika {
/**
* Calculate the maximum allowed length of the next tracking step, based on all
- * ContinuousProcess-es
+ * ContinuousProcess-es.
*
* The maximum allowed step length is the minimum of the allowed track lenght over all
* ContinuousProcess-es in the ProcessSequence.
*
* @return ContinuousProcessStepLength which contains the step length itself in
* LengthType, and a unique identifier of the related ContinuousProcess.
- **/
+ */
template <typename TParticle, typename TTrack>
- ContinuousProcessStepLength getMaxStepLength(TParticle& particle, TTrack& vTrack);
-
- template <typename TParticle>
- GrammageType getInteractionLength(TParticle&& particle) {
- return 1. / getInverseInteractionLength(particle);
- }
+ ContinuousProcessStepLength getMaxStepLength(TParticle&& particle, TTrack&& vTrack);
template <typename TParticle>
- InverseGrammageType getInverseInteractionLength(TParticle&& particle);
-
- template <typename TSecondaryView>
- ProcessReturn selectInteraction(
- TSecondaryView& view, [[maybe_unused]] InverseGrammageType lambda_inv_select,
- [[maybe_unused]] InverseGrammageType lambda_inv_sum =
- InverseGrammageType::zero());
+ CrossSectionType getCrossSection(TParticle&& projectile, Code const targetId,
+ HEPEnergyType const sqrtSnn) const;
template <typename TParticle>
- TimeType getLifetime(TParticle& particle) {
+ TimeType getLifetime(TParticle&& particle) {
return 1. / getInverseLifetime(particle);
}
+ template <typename TSecondaryView, typename TRNG>
+ inline ProcessReturn selectInteraction(
+ TSecondaryView&& view, COMBoost const& boost, HEPEnergyType const sqrtSnn,
+ NuclearComposition const& composition, TRNG&& rng,
+ CrossSectionType const cx_select,
+ CrossSectionType cx_sum = CrossSectionType::zero());
+
template <typename TParticle>
InverseTimeType getInverseLifetime(TParticle&& particle);
// select decay process
template <typename TSecondaryView>
- ProcessReturn selectDecay(
- TSecondaryView& view, [[maybe_unused]] InverseTimeType decay_inv_select,
- [[maybe_unused]] InverseTimeType decay_inv_sum = InverseTimeType::zero());
+ ProcessReturn selectDecay(TSecondaryView&& view, InverseTimeType decay_inv_select,
+ InverseTimeType decay_inv_sum = InverseTimeType::zero());
/**
* static counter to uniquely index (count) all ContinuousProcess in switch sequence.
- **/
+ */
static unsigned int constexpr getNumberOfProcesses() { return numberOfProcesses_; }
#ifdef CORSIKA_UNIT_TESTING
@@ -274,40 +276,40 @@ namespace corsika {
#endif
private:
- TProcess1 A_; /// process/list A, this is a reference, if possible
- TProcess2 B_; /// process/list B, this is a reference, if possible
+ TProcess1 A_; //! process/list A, this is a reference, if possible
+ TProcess2 B_; //! process/list B, this is a reference, if possible
static unsigned int constexpr numberOfProcesses_ = IndexOfProcess1; // static counter
};
/**
- @fn make_sequence
- @ingroup Processes
-
- Factory function to create a ProcessSequence
-
- to construct ProcessSequences in a flexible and dynamic way the
- `sequence` factory functions are provided
-
- Any objects of type
- - BaseProcess
- - ContinuousProcess and
- - InteractionProcess/DecayProcess
- - StackProcess
- - SecondariesProcess
- - BoundaryCrossingProcess
-
- can be assembled into a ProcessSequence, all
- combinatorics are allowed.
-
- The sequence function checks that all its arguments are all of
- types derived from BaseProcess. Also the ProcessSequence itself
- is derived from type BaseProcess
-
- @tparam TProcesses parameter pack with objects of type BaseProcess
- @tparam TProcess1 another BaseProcess
- @param vA needs to derive from BaseProcess
- @param vB paramter-pack, needs to derive BaseProcess
+ * @fn make_sequence
+ * @ingroup Processes
+ *
+ * Factory function to create a ProcessSequence
+ *
+ * to construct ProcessSequences in a flexible and dynamic way the
+ * `sequence` factory functions are provided
+ *
+ * Any objects of type
+ * - BaseProcess
+ * - ContinuousProcess and
+ * - InteractionProcess/DecayProcess
+ * - StackProcess
+ * - SecondariesProcess
+ * - BoundaryCrossingProcess
+ *
+ * can be assembled into a ProcessSequence, all
+ * combinatorics are allowed.
+ *
+ * The sequence function checks that all its arguments are all of
+ * types derived from BaseProcess. Also the ProcessSequence itself
+ * is derived from type BaseProcess.
+ *
+ * @tparam TProcesses parameter pack with objects of type BaseProcess.
+ * @tparam TProcess1 another BaseProcess.
+ * @param vA needs to derive from BaseProcess.
+ * @param vB paramter-pack, needs to derive BaseProcess.
*/
template <typename... TProcesses, typename TProcess1>
ProcessSequence<TProcess1, decltype(make_sequence(std::declval<TProcesses>()...))>
@@ -318,11 +320,11 @@ namespace corsika {
}
/**
- @fn make_sequence
- @ingroup Processes
-
- Factory function to create ProcessSequence
-
+ * @fn make_sequence
+ * @ingroup Processes
+ *
+ * Factory function to create ProcessSequence
+ *
specialization for two input objects (no paramter pack in vB).
@tparam TProcess1 another BaseProcess
diff --git a/corsika/framework/process/SwitchProcessSequence.hpp b/corsika/framework/process/SwitchProcessSequence.hpp
index 09edeea43..d3f01120c 100644
--- a/corsika/framework/process/SwitchProcessSequence.hpp
+++ b/corsika/framework/process/SwitchProcessSequence.hpp
@@ -146,18 +146,15 @@ namespace corsika {
ContinuousProcessStepLength getMaxStepLength(TParticle& particle, TTrack& vTrack);
template <typename TParticle>
- GrammageType getInteractionLength(TParticle&& particle) {
- return 1. / getInverseInteractionLength(particle);
- }
-
- template <typename TParticle>
- InverseGrammageType getInverseInteractionLength(TParticle&& particle);
-
- template <typename TSecondaryView>
- ProcessReturn selectInteraction(
- TSecondaryView& view, [[maybe_unused]] InverseGrammageType lambda_inv_select,
- [[maybe_unused]] InverseGrammageType lambda_inv_sum =
- InverseGrammageType::zero());
+ CrossSectionType getCrossSection(TParticle const& projectile, Code const targetId,
+ HEPEnergyType const sqrtSnn) const;
+
+ template <typename TSecondaryView, typename TRNG>
+ ProcessReturn selectInteraction(TSecondaryView& view, COMBoost const& boost,
+ HEPEnergyType const sqrtSnn,
+ NuclearComposition const& composition, TRNG& rng,
+ CrossSectionType const cx_select,
+ CrossSectionType cx_sum = CrossSectionType::zero());
template <typename TParticle>
TimeType getLifetime(TParticle&& particle) {
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index eac8dde53..d1f9ee785 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -66,6 +66,9 @@ namespace corsika {
//! returns the rotated coordinate system
CoordinateSystemPtr getRotatedCS() const;
+ //! returns the original coordinate system
+ CoordinateSystemPtr getOriginalCS() const;
+
protected:
//! internal method
void setBoost(double coshEta, double sinhEta);
diff --git a/corsika/media/NuclearComposition.hpp b/corsika/media/NuclearComposition.hpp
index e4ad9513a..184853ef2 100644
--- a/corsika/media/NuclearComposition.hpp
+++ b/corsika/media/NuclearComposition.hpp
@@ -20,48 +20,69 @@
namespace corsika {
- /** Describes the composition of matter
- * Allowes and handles the creation of custom matter compositions
- **/
+ /**
+ * Describes the composition of matter
+ * Allowes and handles the creation of custom matter compositions.
+ */
+
class NuclearComposition {
public:
- /** Constructor
+ /**
+ * Constructor
* The constructore takes a list of elements and a list which describe the relative
* amount. Booth lists need to have the same length and the sum all of fractions
- * should be 1. Otherwise an exception is thrown
- * @param pComponents List of particle types
+ * should be 1. Otherwise an exception is thrown.
+ *
+ * @param pComponents List of particle types.
* @param pFractions List of fractions how much each particle contributes. The sum
- * needs to add up to 1
- **/
+ * needs to add up to 1.
+ */
NuclearComposition(std::vector<Code> const& pComponents,
- std::vector<float> const& pFractions);
+ std::vector<double> const& pFractions);
+
+ /**
+ * Returns a vector of the same length as elements in the material with the weighted
+ * return of "func". The typical default application is for cross section weighted
+ * with fraction in the material.
+ *
+ * @tparam TFunction Type of functions for the weights. The type should be
+ * Code -> CrossSectionType.
+ * @param func Functions for reweighting specific elements.
+ * @retval returns the vector with weighted return types of func.
+ */
+ template <typename TFunction>
+ auto getWeighted(TFunction const& func) const;
- /** Sum all all relative composition weighted by func(element)
+ /**
+ * Sum all all relative composition weighted by func(element)
* This function sums all relative compositions given during this classes
- *construction. Each entry is weighted by the user defined function func given to this
- *function.
+ * construction. Each entry is weighted by the user defined function func given to
+ * this function.
+ *
* @tparam TFunction Type of functions for the weights. The type should be
- * Code -> float
- * @param func Functions for reweighting specific elements
- * @retval returns the weighted sum with the type defined by the return type of func
- **/
+ * Code -> double.
+ * @param func Functions for reweighting specific elements.
+ * @retval returns the weighted sum with the type defined by the return type of func.
+ */
template <typename TFunction>
- auto getWeightedSum(TFunction const& func) const;
+ auto getWeightedSum(TFunction const& func) const
+ -> decltype(func(std::declval<Code>()));
- /** Number of elements in the composition array
- * @retval returns the number of elements in the composition array
- **/
+ /**
+ * Number of elements in the composition array
+ * @retval returns the number of elements in the composition array.
+ */
size_t getSize() const;
//! Returns a const reference to the fraction
- std::vector<float> const& getFractions() const;
+ std::vector<double> const& getFractions() const;
//! Returns a const reference to the fraction
std::vector<Code> const& getComponents() const;
double const getAverageMassNumber() const;
template <class TRNG>
Code sampleTarget(std::vector<CrossSectionType> const& sigma,
- TRNG& randomStream) const;
+ TRNG&& randomStream) const;
// Note: when this class ever modifies its internal data, the hash
// must be updated, too!
@@ -73,8 +94,8 @@ namespace corsika {
private:
void updateHash();
- std::vector<float> const numberFractions_; //!< relative fractions of number density
- std::vector<Code> const components_; //!< particle codes of consitutents
+ std::vector<double> const numberFractions_; //!< relative fractions of number density
+ std::vector<Code> const components_; //!< particle codes of consitutents
double const avgMassNumber_;
diff --git a/corsika/modules/Sibyll.hpp b/corsika/modules/Sibyll.hpp
index 44bebf0ff..401766122 100644
--- a/corsika/modules/Sibyll.hpp
+++ b/corsika/modules/Sibyll.hpp
@@ -9,6 +9,41 @@
#pragma once
#include <corsika/modules/sibyll/ParticleConversion.hpp>
-#include <corsika/modules/sibyll/Interaction.hpp>
+#include <corsika/modules/sibyll/InteractionModel.hpp>
#include <corsika/modules/sibyll/Decay.hpp>
-#include <corsika/modules/sibyll/NuclearInteraction.hpp>
+#include <corsika/modules/sibyll/NuclearInteractionModel.hpp>
+
+#include <corsika/framework/process/InteractionProcess.hpp>
+
+/**
+ * @file Sibyll.hpp
+ *
+ * Includes all the parts of the Sibyll model. Defines the InteractoinProcess<TModel>
+ * classes needed for the ProcessSequence.
+ */
+
+namespace corsika::sibyll {
+ /**
+ * @brief sibyll::Interaction is the process for ProcessSequence.
+ *
+ * The sibyll::Model is wrapped as an InteractionProcess here in order
+ * to provide all the functions for ProcessSequence.
+ */
+ class Interaction : public InteractionModel, public InteractionProcess<Interaction> {};
+
+ /**
+ * @brief sibyll::NuclearInteraction is the process for ProcessSequence.
+ *
+ * The sibyll::NuclearInteractionModel is wrapped as an InteractionProcess here in order
+ * to provide all the functions for ProcessSequence.
+ */
+ template <class TEnvironment, class TNucleonModel>
+ class NuclearInteraction
+ : public NuclearInteractionModel<TEnvironment, TNucleonModel>,
+ public InteractionProcess<NuclearInteraction<TEnvironment, TNucleonModel>> {
+ public:
+ NuclearInteraction(TNucleonModel& model, TEnvironment const& env)
+ : NuclearInteractionModel<TNucleonModel, TEnvironment>(model, env) {}
+ };
+
+} // namespace corsika::sibyll
diff --git a/corsika/modules/sibyll/Interaction.hpp b/corsika/modules/sibyll/Interaction.hpp
deleted file mode 100644
index b53fb6308..000000000
--- a/corsika/modules/sibyll/Interaction.hpp
+++ /dev/null
@@ -1,73 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/process/InteractionProcess.hpp>
-#include <tuple>
-
-namespace corsika::sibyll {
-
- class Interaction : public InteractionProcess<Interaction> {
-
- public:
- Interaction(bool const sibyll_printout_on = false);
- ~Interaction();
-
- bool isValidCoMEnergy(HEPEnergyType const ecm) const {
- return (minEnergyCoM_ <= ecm) && (ecm <= maxEnergyCoM_);
- }
- //! sibyll only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
- //! neutrons (p,n == nucleon)
- bool isValidTarget(Code const TargetId) const {
- return (is_nucleus(TargetId) && (get_nucleus_A(TargetId) >= minNuclearTargetA_) &&
- (get_nucleus_A(TargetId) < maxTargetMassNumber_)) ||
- (TargetId == Code::Proton || TargetId == Code::Hydrogen ||
- TargetId == Code::Neutron);
- }
-
- //! returns production and elastic cross section for hadrons in sibyll. Inputs are:
- //! CorsikaId of beam particle, CorsikaId of target particle and center-of-mass
- //! energy. Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
- std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
- Code const, Code const, HEPEnergyType const) const;
-
- template <typename TParticle>
- GrammageType getInteractionLength(TParticle const&) const;
-
- /**
- In this function SIBYLL is called to produce one event. The
- event is copied (and boosted) into the shower lab frame.
- */
-
- template <typename TSecondaries>
- void doInteraction(TSecondaries&);
-
- private:
- unsigned int constexpr getMaxTargetMassNumber() const { return maxTargetMassNumber_; }
- HEPEnergyType getMinEnergyCoM() const { return minEnergyCoM_; }
- HEPEnergyType getMaxEnergyCoM() const { return maxEnergyCoM_; }
-
- default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sibyll");
- const HEPEnergyType minEnergyCoM_ = 10. * 1e9 * electronvolt;
- const HEPEnergyType maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
- static unsigned int constexpr maxTargetMassNumber_ = 18;
- static unsigned int constexpr minNuclearTargetA_ = 4;
-
- // data members
- int count_ = 0;
- int nucCount_ = 0;
- bool sibyll_listing_;
- };
-
-} // namespace corsika::sibyll
-
-#include <corsika/detail/modules/sibyll/Interaction.inl>
diff --git a/corsika/modules/sibyll/NuclearInteraction.hpp b/corsika/modules/sibyll/NuclearInteraction.hpp
deleted file mode 100644
index 719636313..000000000
--- a/corsika/modules/sibyll/NuclearInteraction.hpp
+++ /dev/null
@@ -1,70 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/process/InteractionProcess.hpp>
-
-namespace corsika::sibyll {
-
- class Interaction; // fwd-decl
-
- /**
- *
- *
- **/
- template <class TEnvironment>
- class NuclearInteraction : public InteractionProcess<NuclearInteraction<TEnvironment>> {
-
- public:
- NuclearInteraction(sibyll::Interaction&, TEnvironment const&);
- ~NuclearInteraction();
-
- void initializeNuclearCrossSections();
- void printCrossSectionTable(Code);
- CrossSectionType readCrossSectionTable(int const, Code const, HEPEnergyType const);
- HEPEnergyType getMinEnergyPerNucleonCoM() { return gMinEnergyPerNucleonCoM_; }
- HEPEnergyType getMaxEnergyPerNucleonCoM() { return gMaxEnergyPerNucleonCoM_; }
- unsigned int constexpr getMaxNucleusAProjectile() { return gMaxNucleusAProjectile_; }
- unsigned int constexpr getMaxNFragments() { return gMaxNFragments_; }
- unsigned int constexpr getNEnergyBins() { return gNEnBins_; }
-
- template <typename Particle>
- std::tuple<CrossSectionType, CrossSectionType> getCrossSection(Particle const& p,
- const Code TargetId);
-
- template <typename Particle>
- GrammageType getInteractionLength(Particle const&);
-
- template <typename TSecondaryView>
- void doInteraction(TSecondaryView&);
-
- private:
- int count_ = 0;
- int nucCount_ = 0;
-
- TEnvironment const& environment_;
- sibyll::Interaction& hadronicInteraction_;
- std::map<Code, int> targetComponentsIndex_;
- default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sibyll");
- static unsigned int constexpr gNSample_ =
- 500; // number of samples in MC estimation of cross section
- static unsigned int constexpr gMaxNucleusAProjectile_ = 56;
- static unsigned int constexpr gNEnBins_ = 6;
- static unsigned int constexpr gMaxNFragments_ = 60;
- // energy limits defined by table used for cross section in signuc.f
- // 10**1 GeV to 10**6 GeV
- static HEPEnergyType constexpr gMinEnergyPerNucleonCoM_ = 10. * 1e9 * electronvolt;
- static HEPEnergyType constexpr gMaxEnergyPerNucleonCoM_ = 1.e6 * 1e9 * electronvolt;
- };
-
-} // namespace corsika::sibyll
-
-#include <corsika/detail/modules/sibyll/NuclearInteraction.inl>
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index 2a96ef3df..df1b020a2 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -53,13 +53,13 @@ namespace corsika {
/**
* Set data of new particle.
*
- * @param p parent particle
+ * @param parent parent particle
* @param v tuple containing: PID, Momentum Vector, Position, Time
*
* MomentumVector is only used to determine the DirectionVector, the normalization
* is lost.
*/
- void setParticleData(ParticleInterface<TStackIterator> const& p,
+ void setParticleData(ParticleInterface<TStackIterator> const& parent,
particle_data_type const& v);
/**
@@ -73,11 +73,11 @@ namespace corsika {
/**
* Set data of new particle.
*
- * @param p parent particle
+ * @param parent parent particle
* @param v tuple containing: PID, kinetic Energy, Direction Vector, Position, Time
*
*/
- void setParticleData(ParticleInterface<TStackIterator> const& p,
+ void setParticleData(ParticleInterface<TStackIterator> const& parent,
particle_data_momentum_type const& v);
///! Set particle corsika::Code
diff --git a/tests/common/SetupTestEnvironment.hpp b/tests/common/SetupTestEnvironment.hpp
index fa48ec023..358ed6df2 100644
--- a/tests/common/SetupTestEnvironment.hpp
+++ b/tests/common/SetupTestEnvironment.hpp
@@ -49,7 +49,7 @@ namespace corsika::setup::testing {
world->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, Vector(cs, 0_T, 0_T, BfieldZ), 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{vTargetCode}, std::vector<float>{1.}));
+ NuclearComposition(std::vector<Code>{vTargetCode}, std::vector<double>{1.}));
setup::Environment::BaseNodeType* nodePtr = world.get();
universe.addChild(std::move(world));
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index 9c9918729..b97f95497 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -17,6 +17,8 @@
#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>
@@ -55,8 +57,8 @@ auto make_dummy_env() {
Point{env.getCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
- using MyEmptyModel = Empty<IEmpty>;
- world->setModelProperties<MyEmptyModel>();
+ NuclearComposition const composition({Code::Proton}, {1.});
+ world->setModelProperties<TestEnvironmentInterface>(19.2_g / cube(1_cm), composition);
universe.addChild(std::move(world));
return env;
@@ -86,16 +88,15 @@ public:
class ProcessSplit : public InteractionProcess<ProcessSplit> {
- int calls_ = 0;
-
public:
- template <typename Particle>
- GrammageType getInteractionLength(Particle const&) const {
- return 0_g / square(1_cm);
+ CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const,
+ unsigned int const, unsigned int const) const {
+ return 1_mb;
}
template <typename TView>
- void doInteraction(TView& view) {
+ void doInteraction(TView& view, COMBoost const&, Code, Code, HEPEnergyType,
+ unsigned int, unsigned int) {
++calls_;
auto vP = view.getProjectile();
const HEPEnergyType Ekin = vP.getKineticEnergy();
@@ -106,14 +107,13 @@ public:
}
int getCalls() const { return calls_; }
+
+private:
+ int calls_ = 0;
};
class ProcessCut : public SecondariesProcess<ProcessCut> {
- int count_ = 0;
- int calls_ = 0;
- HEPEnergyType fEcrit;
-
public:
ProcessCut(HEPEnergyType e)
: fEcrit(e) {}
@@ -136,6 +136,11 @@ public:
int getCount() const { return count_; }
int getCalls() const { return calls_; }
+
+private:
+ int count_ = 0;
+ int calls_ = 0;
+ HEPEnergyType fEcrit;
};
TEST_CASE("Cascade", "[Cascade]") {
@@ -153,7 +158,7 @@ TEST_CASE("Cascade", "[Cascade]") {
StackInspector<TestCascadeStack> stackInspect(100, true, E0);
NullModel nullModel;
- const HEPEnergyType Ecrit = 85_MeV;
+ HEPEnergyType const Ecrit = 85_MeV;
ProcessSplit split;
ProcessCut cut(Ecrit);
auto sequence = make_sequence(nullModel, stackInspect, split, cut);
@@ -172,7 +177,6 @@ TEST_CASE("Cascade", "[Cascade]") {
SECTION("full cascade") {
EAS.run();
-
CHECK(cut.getCount() == 2048);
CHECK(cut.getCalls() == 2047); // final particle is still on stack and not yet deleted
CHECK(split.getCalls() == 2047);
diff --git a/tests/framework/testCascade.hpp b/tests/framework/testCascade.hpp
index cb6e0b5a0..8c37a22d2 100644
--- a/tests/framework/testCascade.hpp
+++ b/tests/framework/testCascade.hpp
@@ -9,14 +9,15 @@
#pragma once
#include <corsika/media/Environment.hpp>
-#include <corsika/media/IEmpty.hpp>
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/framework/stack/CombinedStack.hpp>
#include <corsika/framework/stack/SecondaryView.hpp>
#include <corsika/stack/GeometryNodeStackExtension.hpp>
#include <corsika/stack/VectorStack.hpp>
-using TestEnvironmentInterface = corsika::IEmpty;
+using TestEnvironmentInterface = corsika::HomogeneousMedium<corsika::IMediumModel>;
using TestEnvironmentType = corsika::Environment<TestEnvironmentInterface>;
template <typename T>
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index cf9dc523f..39b230e14 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -9,10 +9,15 @@
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/ProcessTraits.hpp>
#include <corsika/framework/process/ContinuousProcessStepLength.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/utility/COMBoost.hpp>
+
+#include <corsika/media/NuclearComposition.hpp>
+
#include <catch2/catch.hpp>
#include <array>
@@ -30,6 +35,12 @@
using namespace corsika;
using namespace std;
+struct DummyRNG {
+ int max() const { return 10; }
+ int min() const { return 0; }
+ double operator()() const { return 0.5; }
+};
+
static int const nData = 10;
// DummyNode is only needed for BoundaryCrossingProcess
@@ -46,6 +57,10 @@ struct DummyStack {};
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 {}
+ HEPEnergyType getEnergy() const { return 10_GeV; }
+ unsigned int getNuclearA() const { return 1; }
};
// there is no real trajectory/track
@@ -193,14 +208,15 @@ public:
}
template <typename TView>
- void doInteraction(TView& v) const {
+ void doInteraction(TView& v, COMBoost const&, Code const, Code const,
+ HEPEnergyType const, unsigned int const, unsigned int const) const {
checkInteract |= 1;
for (int i = 0; i < nData; ++i) v.parent().data_[i] += 1 + i;
}
- template <typename TParticle>
- GrammageType getInteractionLength(TParticle&) const {
- return 10_g / square(1_cm);
+ CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const,
+ unsigned int const, unsigned int const) const {
+ return 10_mb;
}
private:
@@ -219,15 +235,17 @@ public:
}
template <typename TView>
- void doInteraction(TView& v) const {
+ void doInteraction(TView& v, COMBoost const&, Code const, Code const,
+ HEPEnergyType const, unsigned int const, unsigned int const) const {
checkInteract |= 2;
for (int i = 0; i < nData; ++i) v.parent().data_[i] /= 1.1;
CORSIKA_LOG_DEBUG("Process2::doInteraction");
}
- template <typename Particle>
- GrammageType getInteractionLength(Particle&) const {
- CORSIKA_LOG_DEBUG("Process2::GetInteractionLength");
- return 20_g / (1_cm * 1_cm);
+
+ CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const,
+ unsigned int const, unsigned int const) const {
+ CORSIKA_LOG_DEBUG("Process2::getCrossSection");
+ return 20_mb;
}
private:
@@ -246,15 +264,17 @@ public:
}
template <typename TView>
- void doInteraction(TView& v) const {
+ void doInteraction(TView& v, COMBoost const&, Code const, Code const,
+ HEPEnergyType const, unsigned int const, unsigned int const) const {
checkInteract |= 4;
for (int i = 0; i < nData; ++i) v.parent().data_[i] *= 1.01;
CORSIKA_LOG_DEBUG("Process3::doInteraction");
}
- template <typename Particle>
- GrammageType getInteractionLength(Particle&) const {
- CORSIKA_LOG_DEBUG("Process3::GetInteractionLength");
- return 30_g / (1_cm * 1_cm);
+
+ CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const,
+ unsigned int const, unsigned int const) const {
+ CORSIKA_LOG_DEBUG("Process3::getCrossSection");
+ return 30_mb;
}
private:
@@ -280,7 +300,8 @@ public:
return ProcessReturn::Ok;
}
template <typename TView>
- void doInteraction(TView&) const {
+ void doInteraction(TView&, COMBoost const&, Code const, Code const, HEPEnergyType const,
+ unsigned int const, unsigned int const) const {
checkInteract |= 8;
}
@@ -439,7 +460,7 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
DummyData particle;
auto sequence2 = make_sequence(cp1, m2, m3);
- GrammageType const tot = sequence2.getInteractionLength(particle);
+ /*GrammageType const tot = sequence2.getInteractionLength(particle);
InverseGrammageType const tot_inv = sequence2.getInverseInteractionLength(particle);
CORSIKA_LOG_DEBUG(
"lambda_tot={}"
@@ -447,7 +468,7 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
tot, tot_inv);
CHECK(tot / 1_g * square(1_cm) == 12);
- CHECK(tot_inv * 1_g / square(1_cm) == 1. / 12);
+ CHECK(tot_inv * 1_g / square(1_cm) == 1. / 12);*/
globalCount = 0;
}
@@ -595,7 +616,9 @@ TEST_CASE("ProcessSequence General", "ProcessSequence") {
TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::trace);
+
+ CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
/**
* In this example switching is done only by "data_[0]>0", where
@@ -686,24 +709,28 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
CHECK(checkCont == 0b101);
CHECK(checkSec == 0);
- // 1/(30g/cm2) is Process3
- InverseGrammageType lambda_select = .9 / 30. * square(1_cm) / 1_g;
- InverseTimeType time_select = 0.1 / second;
+ // 30_mb is Process3
+ CrossSectionType cx_select = .9 * 30_mb;
+ InverseTimeType time_select = 0.1 / second; // for decay
checkDecay = 0;
checkInteract = 0;
checkSec = 0;
checkCont = 0;
particle.data_[0] = 100; // data positive --> sequence1
- sequence3.selectInteraction(view, lambda_select);
+
+ DummyRNG rng;
+ COMBoost const noBoost({10_GeV, {rootCS, {0_eV, 0_eV, 0_eV}}}, 0_GeV);
+ NuclearComposition const noComposition({Code::Nitrogen}, {1});
+ sequence3.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
sequence3.selectDecay(view, time_select);
CHECK(checkInteract == 0b100); // this is Process3
CHECK(checkDecay == 0b001); // this is Decay1
CHECK(checkCont == 0);
CHECK(checkSec == 0);
- lambda_select = 1.01 / 30. * square(1_cm) / 1_g;
+ cx_select = 1.01 * 30_mb;
checkInteract = 0;
- sequence3.selectInteraction(view, lambda_select);
+ sequence3.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
CHECK(checkInteract == 0b001); // this is Process1
checkDecay = 0;
@@ -711,7 +738,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkSec = 0;
checkCont = 0;
particle.data_[0] = -100; // data negative --> sequence2
- sequence3.selectInteraction(view, lambda_select);
+ sequence3.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
sequence3.selectDecay(view, time_select);
CHECK(checkInteract == 0b010); // this is Process2
CHECK(checkDecay == 0b010); // this is Decay2
@@ -739,7 +766,7 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
checkSec = 0;
checkCont = 0;
particle.data_[0] = -100; // data negative --> sequence1
- sequence4.selectInteraction(view, lambda_select);
+ sequence4.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
sequence4.doSecondaries(view);
sequence4.selectDecay(view, time_select);
sequence4.doSecondaries(view);
@@ -749,11 +776,11 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
CHECK(checkSec == 0);
// check that large "select" value will correctly ignore the call
- lambda_select = 1e5 * square(1_cm) / 1_g;
+ cx_select = 1e5_mb;
time_select = 1e5 / second;
checkDecay = 0;
checkInteract = 0;
- sequence3.selectInteraction(view, lambda_select);
+ sequence3.selectInteraction(view, noBoost, 10_GeV, noComposition, rng, cx_select);
sequence3.selectDecay(view, time_select);
CHECK(checkInteract == 0);
CHECK(checkDecay == 0);
diff --git a/tests/media/CMakeLists.txt b/tests/media/CMakeLists.txt
index ffcbe11d1..c1265193c 100644
--- a/tests/media/CMakeLists.txt
+++ b/tests/media/CMakeLists.txt
@@ -1,5 +1,6 @@
set (test_media_sources
TestMain.cpp
+ testNuclearComposition.cpp
testEnvironment.cpp
testShowerAxis.cpp
testMedium.cpp
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index 43476644e..81c1ac757 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -82,9 +82,6 @@ TEST_CASE("HomogeneousMedium") {
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton}, {1.});
HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
- CHECK(protonComposition.getFractions() == std::vector<float>{1.});
- CHECK(protonComposition.getComponents() == std::vector<Code>{Code::Proton});
-
CHECK_THROWS(NuclearComposition({Code::Proton}, {1.1}));
CHECK_THROWS(NuclearComposition({Code::Proton}, {0.99}));
}
@@ -104,7 +101,7 @@ TEST_CASE("FlatExponential") {
LengthType const length = 2_m;
TimeType const tEnd = length / speed;
- CHECK(medium.getNuclearComposition().getFractions() == std::vector<float>{1.});
+ CHECK(medium.getNuclearComposition().getFractions() == std::vector<double>{1.});
CHECK(medium.getNuclearComposition().getComponents() ==
std::vector<Code>{Code::Proton});
@@ -514,7 +511,7 @@ TEST_CASE("InhomogeneousMedium") {
LengthType const length = tEnd * speed;
- NuclearComposition const composition{{Code::Proton}, {1.f}};
+ NuclearComposition const composition{{Code::Proton}, {1.}};
InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
CORSIKA_LOG_INFO("test={} l={} {} {}", rho.getIntegrateGrammage(trajectory), length,
diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp
index 6e9f347b8..083ccf568 100644
--- a/tests/media/testMagneticField.cpp
+++ b/tests/media/testMagneticField.cpp
@@ -34,8 +34,7 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>;
// the composition we use for the homogenous medium
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
// create a magnetic field vector
Vector B0(gCS, 0_T, 0_T, 0_T);
diff --git a/tests/media/testMedium.cpp b/tests/media/testMedium.cpp
index fb55aaa59..2cbf6a7bb 100644
--- a/tests/media/testMedium.cpp
+++ b/tests/media/testMedium.cpp
@@ -58,8 +58,7 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
const auto density{19.2_g / cube(1_cm)};
// the composition we use for the homogenous medium
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
// the refrative index that we use
const Medium type = corsika::Medium::AirDry1Atm;
diff --git a/tests/media/testNuclearComposition.cpp b/tests/media/testNuclearComposition.cpp
new file mode 100644
index 000000000..7af1f34a2
--- /dev/null
+++ b/tests/media/testNuclearComposition.cpp
@@ -0,0 +1,68 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+
+struct DummyRNG {
+ double v_;
+ DummyRNG(double v)
+ : v_(v) {}
+ int max() const { return 10; }
+ int min() const { return 0; }
+ double operator()() const { return v_; }
+};
+
+TEST_CASE("NuclearComposition") {
+
+ logging::set_level(logging::level::info);
+
+ // incompatible input: wrong vectors
+ CHECK_THROWS(
+ NuclearComposition({Code::Oxygen, Code::Carbon}, {0.20, 0.05, 1 - 0.20 - 0.05}));
+ // incompatible input: wrong fractions
+ CHECK_THROWS(
+ NuclearComposition({Code::Oxygen, Code::Carbon}, {0.21, 0.05, 1 - 0.20 - 0.05}));
+ // incompatible input: wrong fractions
+ CHECK_THROWS(
+ NuclearComposition({Code::Oxygen, Code::Carbon}, {0.19, 0.05, 1 - 0.20 - 0.05}));
+
+ NuclearComposition const testComposition({Code::Oxygen, Code::Carbon, Code::Nitrogen},
+ {0.20, 0.05, 1 - 0.20 - 0.05});
+
+ CHECK(testComposition.getSize() == 3);
+ CHECK(testComposition.getFractions() == std::vector<double>{0.2, 0.05, 1 - 0.2 - 0.05});
+ CHECK(testComposition.getComponents() ==
+ std::vector<Code>{Code::Oxygen, Code::Carbon, Code::Nitrogen});
+
+ CHECK(testComposition.getHash() == 18183071370253166150U);
+ CHECK(testComposition.getAverageMassNumber() == 14.3);
+
+ CHECK(testComposition.getWeighted([](Code) -> double { return 1; }) ==
+ std::vector<double>{0.2, 0.05, 1 - 0.2 - 0.05});
+
+ std::vector<CrossSectionType> const testCX =
+ testComposition.getWeighted([](Code) -> CrossSectionType { return 1_mb; });
+ std::vector<CrossSectionType> const checkCX{0.2_mb, 0.05_mb, 1_mb - 0.2_mb - 0.05_mb};
+ for (auto i1 = testCX.begin(), i2 = checkCX.begin(); i1 != testCX.end(); ++i1, ++i2) {
+ CHECK(*i1 / 1_mb == Approx(*i2 / 1_mb));
+ }
+
+ CHECK(testComposition.getWeightedSum([](Code) -> double { return 1; }) == 1);
+
+ CHECK(testComposition.getWeightedSum([](Code) -> CrossSectionType { return 1_mb; }) ==
+ 1_mb);
+
+ CHECK(testComposition.sampleTarget(testCX, DummyRNG(0.1)) == Code::Oxygen);
+}
diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
index 62422f7a8..53c9a694c 100644
--- a/tests/media/testRefractiveIndex.cpp
+++ b/tests/media/testRefractiveIndex.cpp
@@ -42,8 +42,7 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
const auto density{19.2_g / cube(1_cm)};
// the composition we use for the homogenous medium
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
// the refrative index that we use
const double n{1.000327};
@@ -113,8 +112,7 @@ TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
const auto density{19.2_g / cube(1_cm)};
// the composition we use for the homogenous medium
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
// a new refractive index
const double n0{2};
diff --git a/tests/media/testShowerAxis.cpp b/tests/media/testShowerAxis.cpp
index 9378d9b41..d8285e322 100644
--- a/tests/media/testShowerAxis.cpp
+++ b/tests/media/testShowerAxis.cpp
@@ -36,8 +36,7 @@ auto setupEnvironment(Code vTargetCode) {
using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
theMedium->setModelProperties<MyHomogeneousModel>(
- density,
- NuclearComposition(std::vector<Code>{vTargetCode}, std::vector<float>{1.}));
+ density, NuclearComposition({vTargetCode}, {1.}));
auto const* nodePtr = theMedium.get();
universe.addChild(std::move(theMedium));
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index 16539ba56..d04f3ac80 100644
--- a/tests/modules/testCONEX.cpp
+++ b/tests/modules/testCONEX.cpp
@@ -100,7 +100,7 @@ TEST_CASE("CONEXSourceCut") {
// need to initialize Sibyll, done in constructor:
corsika::sibyll::Interaction sibyll;
- [[maybe_unused]] corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+ [[maybe_unused]] corsika::sibyll::NuclearInteractionModel sibyllNuc(sibyll, env);
CONEXhybrid conex(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton));
conex.initCascadeEquations();
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 140515382..8bc4879f9 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -13,6 +13,7 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
#include <catch2/catch.hpp>
#include <tuple>
@@ -100,68 +101,89 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
return sum;
}
-TEST_CASE("SibyllInteractionInterface", "modules") {
+/*
+calculate COM boost object assuming fixed target collision with projectile
+*/
+COMBoost getCOMboost(HEPEnergyType const& eProjectileLab,
+ MomentumVector const& pProjectileLab,
+ CoordinateSystemPtr const& cs) {
+ // define target
+ // for Sibyll is always a single nucleon
+ // FOR NOW: target is always at rest
+ auto const pTargetLab = MomentumVector(cs, 0_GeV, 0_GeV, 0_GeV);
+ FourVector const P4projLab(eProjectileLab, pProjectileLab);
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Sibyll projectile: +z
+ COMBoost const boost(P4projLab, constants::nucleonMass);
+ return boost;
+}
+
+TEST_CASE("SibyllInterface", "modules") {
logging::set_level(logging::level::info);
+ // the environment and stack should eventually disappear from here
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
auto const& cs = *csPtr;
{ [[maybe_unused]] auto const& env_dummy = env; }
+ auto [stack, viewPtr] = setup::testing::setup_stack(
+ Code::Proton, 0, 0, 10_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ setup::StackView& view = *viewPtr;
+
RNGManager<>::getInstance().registerRandomStream("sibyll");
SECTION("InteractionInterface - valid targets") {
- Interaction model;
+ corsika::sibyll::InteractionModel model;
// sibyll only accepts protons or nuclei with 4<=A<=18 as targets
- CHECK_FALSE(model.isValidTarget(Code::Electron));
- CHECK(model.isValidTarget(Code::Hydrogen));
- CHECK_FALSE(model.isValidTarget(Code::Deuterium));
- CHECK(model.isValidTarget(Code::Helium));
- CHECK_FALSE(model.isValidTarget(Code::Helium3));
- CHECK_FALSE(model.isValidTarget(Code::Iron));
- CHECK(model.isValidTarget(Code::Oxygen));
+ CHECK_THROWS(model.isValid(Code::Proton, Code::Electron, 100_GeV, 1, 1));
+ CHECK_NOTHROW(model.isValid(Code::Proton, Code::Hydrogen, 100_GeV, 1, 1));
+ CHECK_THROWS(model.isValid(Code::Proton, Code::Deuterium, 100_GeV, 1, 2));
+ CHECK_NOTHROW(model.isValid(Code::Proton, Code::Helium, 100_GeV, 1, 4));
+ CHECK_THROWS(model.isValid(Code::Proton, Code::Helium3, 100_GeV, 1, 3));
+ CHECK_THROWS(model.isValid(Code::Proton, Code::Iron, 100_GeV, 1, 56));
+ CHECK_NOTHROW(model.isValid(Code::Proton, Code::Oxygen, 100_GeV, 1, 16));
+ // beam particles
+ CHECK_NOTHROW(model.isValid(Code::Electron, Code::Oxygen, 100_GeV, 1, 1));
+ CHECK_NOTHROW(model.isValid(Code::Nucleus, Code::Oxygen, 100_GeV, 1, 20));
+ // energy too low
+ CHECK_THROWS(model.isValid(Code::Proton, Code::Proton, 9_GeV, 1, 1));
+ CHECK_NOTHROW(model.isValid(Code::Proton, Code::Proton, 11_GeV, 1, 1));
+ // energy too high
+ CHECK_THROWS(model.isValid(Code::Proton, Code::Proton, 1000001_GeV, 1, 1));
+ CHECK_NOTHROW(model.isValid(Code::Proton, Code::Proton, 999999_GeV, 1, 1));
// hydrogen target == proton target == neutron target
auto const [xs_prod_pp, xs_ela_pp] =
- model.getCrossSection(Code::Proton, Code::Proton, 100_GeV);
+ model.getCrossSectionInelEla(Code::Proton, Code::Proton, 100_GeV);
auto const [xs_prod_pn, xs_ela_pn] =
- model.getCrossSection(Code::Proton, Code::Neutron, 100_GeV);
+ model.getCrossSectionInelEla(Code::Proton, Code::Neutron, 100_GeV);
auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] =
- model.getCrossSection(Code::Proton, Code::Hydrogen, 100_GeV);
+ model.getCrossSectionInelEla(Code::Proton, Code::Hydrogen, 100_GeV);
CHECK(xs_prod_pp == xs_prod_pHydrogen);
CHECK(xs_prod_pp == xs_prod_pn);
CHECK(xs_ela_pp == xs_ela_pHydrogen);
CHECK(xs_ela_pn == xs_ela_pHydrogen);
CHECK_THROWS(convertFromSibyll(corsika::sibyll::SibyllCode::Unknown));
-
- // out of range
- // beam particle
- CHECK_THROWS(
- std::get<0>(model.getCrossSection(Code::Electron, Code::Hydrogen, 100_GeV)));
- // target particle
- CHECK(std::get<0>(model.getCrossSection(Code::Proton, Code::Electron, 100_GeV)) ==
- std::numeric_limits<double>::infinity() * 1_mb);
- // energy out of range
- CHECK_THROWS(std::get<0>(model.getCrossSection(Code::Proton, Code::Hydrogen, 5_GeV)));
}
SECTION("InteractionInterface - low energy") {
const HEPEnergyType P0 = 60_GeV;
- auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Proton, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
MomentumVector plab =
MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
- setup::StackView& view = *viewPtr;
-
- auto particle = stack->first();
-
// also print particles after sibyll was called
- Interaction model(true);
-
- model.doInteraction(view);
+ 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);
+ view.clear();
+ COMBoost boost = getCOMboost(Elab, plab, cs);
+ model.doInteraction(view, boost, Code::Proton, Code::Oxygen, sqrtSnn, 0,
+ get_nucleus_A(Code::Oxygen));
auto const pSum = sumMomentum(view, cs);
/*
@@ -227,82 +249,34 @@ TEST_CASE("SibyllInteractionInterface", "modules") {
CHECK((pSum - plab).getNorm() / 1_GeV ==
Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
- [[maybe_unused]] GrammageType const length = model.getInteractionLength(particle);
- CHECK(length / 1_g * 1_cm * 1_cm == Approx(88.7).margin(0.1));
+ [[maybe_unused]] CrossSectionType const cx = model.getCrossSection(
+ Code::Proton, Code::Oxygen, sqrtSnn, 0, get_nucleus_A(Code::Oxygen));
+ 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
// changes"
}
- SECTION("InteractionInterface - energy too low") {
-
- const HEPEnergyType P0 = 5_GeV;
- auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Proton, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
- MomentumVector plab =
- MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
- setup::StackView& view = *viewPtr;
-
- auto particle = stack->first();
-
- Interaction model;
- CHECK_THROWS(model.doInteraction(view));
-
- [[maybe_unused]] GrammageType const length = model.getInteractionLength(particle);
- CHECK(model.getInteractionLength(particle) / 1_g * 1_cm * 1_cm ==
- std::numeric_limits<double>::infinity());
- }
-
- SECTION("InteractionInterface - energy too high") {
-
- const HEPEnergyType P0 = 1000_EeV;
- auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Proton, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
- { [[maybe_unused]] auto const& dummy1 = stack; }
- MomentumVector plab =
- MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
- setup::StackView& view = *viewPtr;
-
- Interaction model;
- CHECK_THROWS(model.doInteraction(view));
- }
-
- SECTION("InteractionInterface - target nucleus out of range") {
- auto [env1, csPtr1, nodePtr1] = setup::testing::setup_environment(Code::Argon);
- { [[maybe_unused]] auto const& dummy1 = env1; }
- auto const& cs1 = *csPtr1;
- const HEPEnergyType P0 = 150_GeV;
- auto [stack, viewPtr] = setup::testing::setup_stack(
- Code::Electron, P0, (setup::Environment::BaseNodeType* const)nodePtr1, cs1);
- { [[maybe_unused]] auto const& dummy1 = stack; }
- MomentumVector plab = MomentumVector(
- cs1, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
- setup::StackView& view = *viewPtr;
-
- Interaction model;
- CHECK_THROWS(model.doInteraction(view));
- }
-
SECTION("NuclearInteractionInterface") {
- auto [stack, viewPtr] =
- setup::testing::setup_stack(get_nucleus_code(8, 4), 900_GeV,
- (setup::Environment::BaseNodeType* const)nodePtr, cs);
- setup::StackView& view = *viewPtr;
- auto particle = stack->first();
-
- Interaction hmodel;
- NuclearInteraction model(hmodel, *env);
-
- model.doInteraction(view);
- [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
+ HEPMomentumType const P0 = 2500_GeV; // per nucleon
+ MomentumVector plab = MomentumVector(cs, {P0, 0_eV, 0_eV}); // per nucleon
+ corsika::sibyll::InteractionModel hmodel;
+ NuclearInteractionModel model(hmodel, *env);
+ HEPEnergyType const ElabNuc =
+ sqrt(static_pow<2>(P0 * 8) + static_pow<2>(get_nucleus_mass(8, 4))) / 8;
+ HEPEnergyType const sqrtSnn = sqrt((ElabNuc + constants::nucleonMass + P0) *
+ (ElabNuc + constants::nucleonMass - P0));
+ model.doInteraction(view, getCOMboost(ElabNuc, plab, cs), Code::Nucleus, Code::Oxygen,
+ sqrtSnn, 8, get_nucleus_A(Code::Oxygen));
+ CrossSectionType const cx = model.getCrossSection(
+ Code::Nucleus, Code::Oxygen, sqrtSnn, 8, get_nucleus_A(Code::Oxygen));
// 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(length / 1_g * 1_cm * 1_cm ==
- Approx(31).margin(5)); // this is not physics validation
+ CHECK(cx / 1_mb == Approx(870).margin(60)); // 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(90).margin(90)); // this is not physics validation
}
}
@@ -341,7 +315,7 @@ TEST_CASE("SibyllDecayInterface", "modules") {
Decay model;
model.printDecayConfig();
- [[maybe_unused]] const TimeType time = model.getLifetime(particle);
+ [[maybe_unused]] TimeType const time = model.getLifetime(particle);
auto const gamma = particle.getEnergy() / particle.getMass();
CHECK(time == get_lifetime(Code::Lambda0) * gamma);
model.doDecay(view);
@@ -372,7 +346,7 @@ TEST_CASE("SibyllDecayInterface", "modules") {
CHECK(model.isDecayHandled(Code::PiMinus));
CHECK_FALSE(model.isDecayHandled(Code::KPlus));
- const std::vector<Code> particleTestList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
+ std::vector<Code> const particleTestList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
Code::Lambda0Bar, Code::D0Bar};
// setup decays
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index d1d967d6f..828397c1b 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -141,19 +141,19 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
MagneticFieldVector magneticfield(cs, 0_T, 0_T, Bfield);
target->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<double>{1.}));
target_neutral->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<double>{1.}));
target_2->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<double>{1.}));
target_2_behind->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<double>{1.}));
target_2_partly_behind->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<double>{1.}));
auto* targetPtr = target.get();
auto* targetPtr_2 = target_2.get();
auto* targetPtr_neutral = target_neutral.get();
@@ -287,7 +287,7 @@ TEST_CASE("TrackingLeapFrogCurved") {
MagneticFieldVector magneticfield(cs, 100_T, 0_T, 0_uT);
target->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, magneticfield, 1_g / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<float>{1.}));
+ NuclearComposition(std::vector<Code>{Code::Oxygen}, std::vector<double>{1.}));
auto* targetPtr = target.get();
worldPtr->addChild(std::move(target));
--
GitLab