diff --git a/corsika/detail/framework/process/InteractionProcess.hpp b/corsika/detail/framework/process/InteractionProcess.hpp
index 462556675b8526acafd06dfed68861afbaf045af..7319027194015b5caaf23946a1c12d2b4fd79c84 100644
--- a/corsika/detail/framework/process/InteractionProcess.hpp
+++ b/corsika/detail/framework/process/InteractionProcess.hpp
@@ -40,9 +40,9 @@ namespace corsika {
public:
/**
- @name traits results
- @{
- */
+ * @name traits results
+ * @{
+ */
using type = decltype(test<std::decay_t<TProcess>>(nullptr));
static const bool value = type::value;
//! @}
@@ -54,11 +54,11 @@ namespace corsika {
has_method_doInteract<TProcess, TReturn, TTemplate, TArgs...>::value;
/**
- * traits test for InteractionProcess::getInteractionLength method.
+ * traits test for TEMPLATED InteractionProcess::getCrossSection method (PROPOSAL).
*/
- template <class TProcess, typename TReturn, typename... TArgs>
- struct has_method_getInteractionLength
+ template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+ struct has_method_getCrossSectionTemplate
: public detail::has_method_signature<TReturn, TArgs...> {
///! method signature
@@ -70,12 +70,12 @@ namespace corsika {
//! templated parameter option
template <class T>
- static decltype(testSignature(&T::template getInteractionLength<TArgs...>)) test(
+ static decltype(testSignature(&T::template getCrossSection<TTemplate>)) test(
std::nullptr_t);
//! non templated parameter option
template <class T>
- static decltype(testSignature(&T::getInteractionLength)) test(std::nullptr_t);
+ static decltype(testSignature(&T::getCrossSection)) test(std::nullptr_t);
public:
/**
@@ -88,9 +88,9 @@ namespace corsika {
};
//! value traits type shortcut
- template <class TProcess, typename TReturn, typename... TArgs>
- bool constexpr has_method_getInteractionLength_v =
- has_method_getInteractionLength<TProcess, TReturn, TArgs...>::value;
+ template <class TProcess, typename TReturn, typename TTemplate, typename... TArgs>
+ bool constexpr has_method_getCrossSectionTemplate_v =
+ has_method_getCrossSectionTemplate<TProcess, TReturn, TTemplate, TArgs...>::value;
/**
* traits test for InteractionProcess::getCrossSection method.
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 5711ff24991ce6113d16b74cf487c354f3a205d6..c8bd3d0a4b9b9a2cb065a21cf45c732775a9d2a1 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -322,9 +322,22 @@ namespace corsika {
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>) {
- tot += A_.getCrossSection(projectile.getPID(), targetId,
- {projectile.getEnergy(), projectile.getMomentum()},
- targetP4);
+
+ bool constexpr has_signature_cx1 =
+ has_method_getCrossSection_v<TProcess1, // process object
+ CrossSectionType, // return type
+ Code, Code, // parameters
+ FourMomentum const&, FourMomentum const&>;
+
+ if constexpr (has_signature_cx1) {
+ tot += A_.getCrossSection(projectile.getPID(), targetId,
+ {projectile.getEnergy(), projectile.getMomentum()},
+ targetP4);
+ } else { // for PROPOSAL
+ tot += A_.getCrossSection(projectile, projectile.getPID(),
+ {projectile.getEnergy(), projectile.getMomentum()});
+ }
+
} else if constexpr (process1_type::is_process_sequence) {
tot += A_.getCrossSection(projectile, targetId, targetP4);
}
@@ -332,9 +345,22 @@ namespace corsika {
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>) {
- tot += B_.getCrossSection(projectile.getPID(), targetId,
- {projectile.getEnergy(), projectile.getMomentum()},
- targetP4);
+
+ bool constexpr has_signature_cx1 =
+ has_method_getCrossSection_v<TProcess2, // process object
+ CrossSectionType, // return type
+ Code, Code, // parameters
+ FourMomentum const&, FourMomentum const&>;
+
+ if constexpr (has_signature_cx1) {
+ tot += B_.getCrossSection(projectile.getPID(), targetId,
+ {projectile.getEnergy(), projectile.getMomentum()},
+ targetP4);
+ } else { // for PROPOSAL
+ tot += B_.getCrossSection(projectile, projectile.getPID(),
+ {projectile.getEnergy(), projectile.getMomentum()});
+ }
+
} else if constexpr (process2_type::is_process_sequence) {
tot += B_.getCrossSection(projectile, targetId, targetP4);
}
@@ -443,52 +469,74 @@ namespace corsika {
// get cross section vector for all material components
// for selected process A
- static_assert(
+
+ bool constexpr has_signature_cx1 =
has_method_getCrossSection_v<TProcess1, // process object
CrossSectionType, // return type
Code, Code, // parameters
- FourMomentum const&, FourMomentum const&>,
- "TProcess1 has no method with correct signature \"CrossSectionType "
- "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
- "const&)\" required by "
- "InteractionProcess<TProcess1>. ");
-
- std::vector<CrossSectionType> const weightedCrossSections =
- composition.getWeighted([=](Code const targetId) -> CrossSectionType {
- FourMomentum const targetP4(
- get_mass(targetId),
- MomentumVector(projectile.getMomentum().getCoordinateSystem(),
- {0_GeV, 0_GeV, 0_GeV}));
- return A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
- });
-
- cx_sum += std::accumulate(weightedCrossSections.cbegin(),
- weightedCrossSections.cend(), CrossSectionType::zero());
+ FourMomentum const&, FourMomentum const&>;
+ bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+ has_method_getCrossSectionTemplate_v<
+ TProcess1, // process object
+ CrossSectionType, // return type
+ decltype(projectile) const&, // template argument
+ decltype(projectile) const&, // parameters
+ Code, FourMomentum const&>;
+ static_assert((has_signature_cx1 || has_signature_cx2),
+ "TProcess1 has no method with correct signature \"CrossSectionType "
+ "getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
+ "const&)\" required by "
+ "InteractionProcess<TProcess1>. ");
+
+ std::vector<CrossSectionType> weightedCrossSections;
+ if constexpr (has_signature_cx1) {
+ /*std::vector<CrossSectionType> const*/ weightedCrossSections =
+ composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+ return A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+ });
+
+ cx_sum +=
+ std::accumulate(weightedCrossSections.cbegin(),
+ weightedCrossSections.cend(), CrossSectionType::zero());
+
+ } else { // this is for PROPOSAL
+ cx_sum += A_.template getCrossSection(projectile, projectileId, projectileP4);
+ }
// check if we should execute THIS process and then EXIT
if (cx_select <= cx_sum) {
- // now also sample targetId from weighted cross sections
- Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
- FourMomentum const targetP4(
- get_mass(targetId),
- MomentumVector(projectile.getMomentum().getCoordinateSystem(),
- {0_GeV, 0_GeV, 0_GeV}));
-
- // interface checking on TProcess1
- static_assert(
- has_method_doInteract_v<TProcess1, // process object
- void, // return type
- TSecondaryView, // template argument
- TSecondaryView&, // method parameters
- Code, Code, FourMomentum const&,
- FourMomentum const&>,
- "TProcess1 has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, "
- "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
- "InteractionProcess<TProcess1>. ");
-
- A_.template doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+ if constexpr (has_signature_cx1) {
+ // now also sample targetId from weighted cross sections
+ Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+
+ // interface checking on TProcess1
+ static_assert(
+ has_method_doInteract_v<TProcess1, // process object
+ void, // return type
+ TSecondaryView, // template argument
+ TSecondaryView&, // method parameters
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>,
+ "TProcess1 has no method with correct signature \"void "
+ "doInteraction<TSecondaryView>(TSecondaryView&, "
+ "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+ "InteractionProcess<TProcess1>. ");
+
+ A_.template doInteraction(view, projectileId, targetId, projectileP4,
+ targetP4);
+
+ } else { // this is for PROPOSAL
+ A_.template doInteraction(view, projectileId, projectileP4);
+ }
return ProcessReturn::Interacted;
}
@@ -508,52 +556,72 @@ namespace corsika {
Code const projectileId = projectile.getPID();
// get cross section vector for all material components, for selected process B
- static_assert(has_method_getCrossSection_v<TProcess2, // process object
- CrossSectionType, // return type
- Code, Code, FourMomentum const&,
- FourMomentum const&>, // parameters
+ bool constexpr has_signature_cx1 =
+ has_method_getCrossSection_v<TProcess2, // process object
+ CrossSectionType, // return type
+ Code, Code, // parameters
+ FourMomentum const&, FourMomentum const&>;
+ bool constexpr has_signature_cx2 = // needed for PROPOSAL interface
+ has_method_getCrossSectionTemplate_v<
+ TProcess2, // process object
+ CrossSectionType, // return type
+ decltype(*projectile) const&, // template argument
+ decltype(*projectile) const&, // parameters
+ Code, // parameters
+ FourMomentum const&>;
+ static_assert((has_signature_cx1 || has_signature_cx2),
"TProcess2 has no method with correct signature \"CrossSectionType "
"getCrossSection(Code, Code, FourMomentum const&, FourMomentum "
"const&)\" required by "
"InteractionProcess<TProcess1>. ");
- std::vector<CrossSectionType> const weightedCrossSections =
- composition.getWeighted([=](Code const targetId) -> CrossSectionType {
- FourMomentum const targetP4(
- get_mass(targetId),
- MomentumVector(projectile.getMomentum().getCoordinateSystem(),
- {0_GeV, 0_GeV, 0_GeV}));
- return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
- });
-
- cx_sum += std::accumulate(weightedCrossSections.begin(),
- weightedCrossSections.end(), CrossSectionType::zero());
+ std::vector<CrossSectionType> weightedCrossSections;
+ if constexpr (has_signature_cx1) {
+ /* std::vector<CrossSectionType> const*/ weightedCrossSections =
+ composition.getWeighted([=](Code const targetId) -> CrossSectionType {
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+ return B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+ });
+
+ cx_sum +=
+ std::accumulate(weightedCrossSections.begin(), weightedCrossSections.end(),
+ CrossSectionType::zero());
+ } else { // this is for PROPOSAL
+ cx_sum += B_.template getCrossSection(projectile, projectileId, projectileP4);
+ }
// check if we should execute THIS process and then EXIT
if (cx_select <= cx_sum) {
- // now also sample targetId from weighted cross sections
- Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
- FourMomentum const targetP4(
- get_mass(targetId),
- MomentumVector(projectile.getMomentum().getCoordinateSystem(),
- {0_GeV, 0_GeV, 0_GeV}));
-
- // interface checking on TProcess2
- static_assert(
- has_method_doInteract_v<TProcess2, // process object
- void, // return type
- TSecondaryView, // template argument
- TSecondaryView&, // method parameters
- Code, Code, FourMomentum const&,
- FourMomentum const&>,
- "TProcess1 has no method with correct signature \"void "
- "doInteraction<TSecondaryView>(TSecondaryView&, "
- "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
- "InteractionProcess<TProcess2>. ");
-
- B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
-
+ if constexpr (has_signature_cx1) {
+
+ // now also sample targetId from weighted cross sections
+ Code const targetId = composition.sampleTarget(weightedCrossSections, rng);
+ FourMomentum const targetP4(
+ get_mass(targetId),
+ MomentumVector(projectile.getMomentum().getCoordinateSystem(),
+ {0_GeV, 0_GeV, 0_GeV}));
+
+ // interface checking on TProcess2
+ static_assert(
+ has_method_doInteract_v<TProcess2, // process object
+ void, // return type
+ TSecondaryView, // template argument
+ TSecondaryView&, // method parameters
+ Code, Code, FourMomentum const&,
+ FourMomentum const&>,
+ "TProcess1 has no method with correct signature \"void "
+ "doInteraction<TSecondaryView>(TSecondaryView&, "
+ "Code, Code, FourMomentum const&, FourMomentum const&)\" required for "
+ "InteractionProcess<TProcess2>. ");
+
+ B_.doInteraction(view, projectileId, targetId, projectileP4, targetP4);
+ } else { // this is for PROPOSAL
+ B_.doInteraction(view, projectileId, projectileP4);
+ }
return ProcessReturn::Interacted;
}
}
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index c450594146f13a7afaa809c99d0f3354d9427b02..e81f9984b7b880ddd208556b087df0d0302b3d59 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -39,24 +39,26 @@ namespace corsika::proposal {
auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Look which interactions take place and build the corresponding
- // interaction and secondarie builder. The interaction integral will
+ // interaction and secondary builder. The interaction integral will
// interpolated too and saved in the calc map by a key build out of a hash
// of composed of the component and particle code.
auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
- calc[std::make_pair(comp.getHash(), code)] = std::make_tuple(
+ calc_[std::make_pair(comp.getHash(), code)] = std::make_tuple(
PROPOSAL::make_secondaries(inter_types, particle[code], media.at(comp.getHash())),
PROPOSAL::make_interaction(c, true));
}
template <typename TStackView>
- inline ProcessReturn Interaction::doInteraction(TStackView& view) {
+ inline ProcessReturn Interaction::doInteraction(TStackView& view,
+ Code const projectileId,
+ FourMomentum const& projectileP4) {
auto const projectile = view.getProjectile();
- if (canInteract(projectile.getPID())) {
+ if (canInteract(projectileId)) {
// get or build corresponding calculators
- auto c = getCalculator(projectile, calc);
+ auto c = getCalculator(projectile, calc_);
// get the rates of the interaction types for every component.
std::uniform_real_distribution<double> distr(0., 1.);
@@ -103,14 +105,27 @@ namespace corsika::proposal {
}
template <typename TParticle>
- inline GrammageType Interaction::getInteractionLength(TParticle const& projectile) {
-
- if (canInteract(projectile.getPID())) {
- auto c = getCalculator(projectile, calc);
- return std::get<eINTERACTION>(c->second)->MeanFreePath(projectile.getEnergy() /
- 1_MeV) *
- 1_g / (1_cm * 1_cm);
+ inline CrossSectionType Interaction::getCrossSection(TParticle const& projectile,
+ Code const projectileId,
+ FourMomentum const& projectileP4) {
+
+ // ==============================================
+ // this block better diappears. RU 26.10.2021
+ //
+ // determine the volume where the particle is (last) known to be
+ auto const* currentLogicalNode = projectile.getNode();
+ NuclearComposition const& composition =
+ currentLogicalNode->getModelProperties().getNuclearComposition();
+ auto const meanMass = composition.getAverageMassNumber() * constants::u;
+ // ==============================================
+
+ if (canInteract(projectileId)) {
+ auto c = getCalculator(projectile, calc_);
+ return meanMass / (std::get<eINTERACTION>(c->second)->MeanFreePath(
+ projectileP4.getTimeLikeComponent() / 1_MeV) *
+ 1_g / (1_cm * 1_cm));
}
- return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+
+ return CrossSectionType::zero();
}
} // namespace corsika::proposal
diff --git a/corsika/detail/modules/proposal/ProposalProcessBase.inl b/corsika/detail/modules/proposal/ProposalProcessBase.inl
index c7c38c10423ab6fe270c0851f7b77ca21d8ab73c..936798ea1cd05340ba2f0c4b397d5af9c47a068e 100644
--- a/corsika/detail/modules/proposal/ProposalProcessBase.inl
+++ b/corsika/detail/modules/proposal/ProposalProcessBase.inl
@@ -56,8 +56,8 @@ namespace corsika::proposal {
PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str();
}
- inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const
- noexcept {
+ inline size_t ProposalProcessBase::hash::operator()(
+ const calc_key_t& p) const noexcept {
return p.first ^ std::hash<Code>{}(p.second);
}
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 9e13f4de82834de41ad14f411583169ec811a79f..5c495e69375dd6b7b767a8501c07611ee5e649f1 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -80,22 +80,22 @@ namespace corsika::pythia8 {
Pythia8::Pythia::particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
}
- inline void Interaction::configureLabFrameCollision(Code const BeamId,
- Code const TargetId,
+ inline void Interaction::configureLabFrameCollision(Code const projectileId,
+ Code const targetId,
HEPEnergyType const BeamEnergy) {
// Pythia configuration of the current event
// very clumsy. I am sure this can be done better..
// set beam
// beam id for pythia
- auto const pdgBeam = static_cast<int>(get_PDG(BeamId));
+ auto const pdgBeam = static_cast<int>(get_PDG(projectileId));
std::stringstream stBeam;
stBeam << "Beams:idA = " << pdgBeam;
Pythia8::Pythia::readString(stBeam.str());
// set target
- auto pdgTarget = static_cast<int>(get_PDG(TargetId));
+ auto pdgTarget = static_cast<int>(get_PDG(targetId));
// replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
- if (TargetId == Code::Hydrogen) pdgTarget = static_cast<int>(get_PDG(Code::Proton));
+ if (targetId == Code::Hydrogen) pdgTarget = static_cast<int>(get_PDG(Code::Proton));
std::stringstream stTarget;
stTarget << "Beams:idB = " << pdgTarget;
Pythia8::Pythia::readString(stTarget.str());
@@ -116,26 +116,34 @@ namespace corsika::pythia8 {
// LCOV_EXCL_STOP
}
- inline bool Interaction::canInteract(Code const pCode) {
+ inline bool Interaction::canInteract(Code const pCode) const {
return pCode == Code::Proton || pCode == Code::Neutron || pCode == Code::AntiProton ||
pCode == Code::AntiNeutron || pCode == Code::PiMinus || pCode == Code::PiPlus;
}
- inline std::tuple<CrossSectionType, CrossSectionType> Interaction::getCrossSection(
- Code const BeamId, Code const TargetId, HEPEnergyType const CoMenergy) {
+ inline std::tuple<CrossSectionType, CrossSectionType>
+ Interaction::getCrossSectionInelEla(Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+
+ HEPEnergyType const CoMenergy = (projectileP4 + targetP4).getNorm();
+
// interaction possible in pythia?
- if (TargetId == Code::Proton || TargetId == Code::Hydrogen) {
- if (canInteract(BeamId) && isValidCoMEnergy(CoMenergy)) {
+ if (targetId == Code::Proton || targetId == Code::Hydrogen) {
+ if (canInteract(projectileId) && isValidCoMEnergy(CoMenergy)) {
// input particle PDG
- auto const pdgCodeBeam = static_cast<int>(get_PDG(BeamId));
- auto const pdgCodeTarget = static_cast<int>(get_PDG(TargetId));
+ auto const pdgCodeBeam = static_cast<int>(get_PDG(projectileId));
+ auto const pdgCodeTarget = static_cast<int>(get_PDG(targetId));
double const ecm = CoMenergy / 1_GeV;
+ //! @todo: remove this const_cast, when Pythia8 becomes const-correct! CHECK!
+ Pythia8::SigmaTotal& sigma = *const_cast<Pythia8::SigmaTotal*>(&sigma_);
+
// calculate cross section
- sigma_.calc(pdgCodeBeam, pdgCodeTarget, ecm);
- if (sigma_.hasSigmaTot()) {
- double const sigEla = sigma_.sigmaEl();
- double const sigProd = sigma_.sigmaTot() - sigEla;
+ sigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
+ if (sigma.hasSigmaTot()) {
+ double const sigEla = sigma.sigmaEl();
+ double const sigProd = sigma.sigmaTot() - sigEla;
return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
@@ -153,92 +161,24 @@ namespace corsika::pythia8 {
}
}
- template <typename TParticle>
- inline GrammageType Interaction::getInteractionLength(TParticle const& particle) {
-
- // coordinate system, get global frame of reference
- MomentumVector const& pMomentum = particle.getMomentum();
- CoordinateSystemPtr const& labCS = pMomentum.getCoordinateSystem();
-
- Code corsikaBeamId = particle.getPID();
-
- // beam particles for pythia : 1, 2, 3 for p, pi, k
- // read from cross section code table
- bool const kInteraction = canInteract(corsikaBeamId);
-
- // FOR NOW: assume target is at rest
- MomentumVector pTarget(labCS, {0_GeV, 0_GeV, 0_GeV});
-
- // total momentum and energy
- HEPEnergyType Elab = particle.getEnergy() + constants::nucleonMass;
- MomentumVector pTotLab(labCS, {0_GeV, 0_GeV, 0_GeV});
- pTotLab += pMomentum;
- 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: {}",
- particle.getEnergy() / 1_GeV, kInteraction, particle.getPID());
-
- // TODO: move limits into variables
- if (kInteraction && Elab >= 8.5_GeV && 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 = particle.getNode();
- auto const mediumComposition =
- currentNode->getModelProperties().getNuclearComposition();
- // determine average interaction length
-
- auto const weightedProdCrossSection =
- mediumComposition.getWeightedSum([=](auto vTargetID) {
- return std::get<0>(this->getCrossSection(corsikaBeamId, vTargetID, ECoM));
- });
-
- CORSIKA_LOG_DEBUG(
- "Interaction: IntLength: weighted CrossSection (mb): {} "
- "Interaction: IntLength: average mass number: {} ",
- weightedProdCrossSection / 1_mb, mediumComposition.getAverageMassNumber());
-
- // 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);
- }
-
template <class TView>
- inline void Interaction::doInteraction(TView& view) {
+ inline void Interaction::doInteraction(TView& view, Code const projectileId,
+ Code const targetId, FourMomentum const&,
+ FourMomentum const&) {
auto projectile = view.getProjectile();
- const auto corsikaBeamId = projectile.getPID();
CORSIKA_LOG_DEBUG(
"Pythia::Interaction: "
"DoInteraction: {} interaction? ",
- corsikaBeamId, corsika::pythia8::Interaction::canInteract(corsikaBeamId));
+ projectileId, corsika::pythia8::Interaction::canInteract(projectileId));
- if (is_nucleus(corsikaBeamId)) {
+ if (is_nucleus(projectileId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
}
- if (corsika::pythia8::Interaction::canInteract(corsikaBeamId)) {
+ if (corsika::pythia8::Interaction::canInteract(projectileId)) {
// define projectile
HEPEnergyType const eProjectileLab = projectile.getEnergy();
@@ -299,32 +239,8 @@ namespace corsika::pythia8 {
// 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<si::CrossSectionType> cross_section_of_components(compVec.size());
-
- for (size_t i = 0; i < compVec.size(); ++i) {
- auto const targetId = compVec[i];
- auto const [sigProd, sigEla] = getCrossSection(corsikaBeamId, targetId, Ecm);
- [[maybe_unused]] auto const& dummy_sigEla = sigEla;
- cross_section_of_components[i] = sigProd;
- }
-
- auto const corsikaTargetId =
- mediumComposition.sampleTarget(cross_section_of_components, RNG_);
- CORSIKA_LOG_DEBUG("Interaction: target selected: {}", corsikaTargetId);
-
- if (corsikaTargetId != Code::Hydrogen && corsikaTargetId != Code::Neutron &&
- corsikaTargetId != Code::Proton)
+ if (targetId != Code::Hydrogen && targetId != Code::Neutron &&
+ targetId != Code::Proton)
throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
CORSIKA_LOG_DEBUG(
@@ -343,7 +259,7 @@ namespace corsika::pythia8 {
} else {
count_++;
- configureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
+ configureLabFrameCollision(projectileId, targetId, eProjectileLab);
// create event in pytia. LCOV_EXCL_START: we don't validate pythia8 internals
if (!Pythia8::Pythia::next())
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 06111acf8c0f534adbddf1aabfa91f1a0aac8338..3b647ec8238585a4633d1135277981d197a6d8d7 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -50,111 +50,109 @@ 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,
@@ -190,6 +188,15 @@ namespace corsika {
*/
ProcessSequence(TProcess1 in_A, TProcess2 in_B);
+ /**
+ * List of all BoundaryProcess.
+ *
+ * @tparam TParticle
+ * @param particle The particle.
+ * @param from Volume the particle is exiting.
+ * @param to Volume the particle is entering.
+ * @return ProcessReturn
+ */
template <typename TParticle>
ProcessReturn doBoundaryCrossing(TParticle& particle,
typename TParticle::node_type const& from,
@@ -199,6 +206,15 @@ namespace corsika {
ProcessReturn doContinuous(TParticle& particle, TTrack& vT,
ContinuousProcessIndex const limitID);
+ /**
+ * Process all secondaries in TSecondaries.
+ *
+ * The seondaries produced by other processes and accessible via TSecondaries
+ * are processed by all SecondariesProcesse via a call here.
+ *
+ * @tparam TSecondaries
+ * @param vS
+ */
template <typename TSecondaries>
void doSecondaries(TSecondaries& vS);
@@ -243,6 +259,15 @@ namespace corsika {
template <typename TParticle, typename TTrack>
ContinuousProcessStepLength getMaxStepLength(TParticle&& particle, TTrack&& vTrack);
+ /**
+ * @brief Calculates the cross section of a projectile with a target.
+ *
+ * @tparam TParticle
+ * @param projectile
+ * @param targetId
+ * @param targetP4
+ * @return CrossSectionType
+ */
template <typename TParticle>
CrossSectionType getCrossSection(TParticle const& projectile, Code const targetId,
FourMomentum const& targetP4) const;
@@ -253,7 +278,7 @@ namespace corsika {
}
/**
- * @brief Selects one concrete InteractionProcess and samples a target nucleus from
+ * Selects one concrete InteractionProcess and samples a target nucleus from
* the material.
*
* The selectInteraction method statically loops over all active InteractionProcess
diff --git a/corsika/modules/proposal/Interaction.hpp b/corsika/modules/proposal/Interaction.hpp
index b6f2fabbab60f772f69eb83bd9f78697a08fb8f8..3b769bc34281487885aab4f3e022e21a54fe564b 100644
--- a/corsika/modules/proposal/Interaction.hpp
+++ b/corsika/modules/proposal/Interaction.hpp
@@ -13,6 +13,7 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/framework/process/ProcessReturn.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/random/UniformRealDistribution.hpp>
@@ -32,7 +33,7 @@ namespace corsika::proposal {
std::unique_ptr<PROPOSAL::Interaction>>;
std::unordered_map<calc_key_t, calculator_t, hash>
- calc; //!< Stores the secondaries and interaction calculators.
+ calc_; //!< Stores the secondaries and interaction calculators.
//!
//! Build the secondaries and interaction calculators and add it to calc.
@@ -53,13 +54,15 @@ namespace corsika::proposal {
//! produce the corresponding secondaries and store them on the particle stack.
//!
template <typename TSecondaryView>
- ProcessReturn doInteraction(TSecondaryView&);
+ ProcessReturn doInteraction(TSecondaryView&, Code const projectileId,
+ FourMomentum const& projectileP4);
//!
- //! Calculates the mean free path length
+ //! Calculates and returns the cross section.
//!
template <typename TParticle>
- GrammageType getInteractionLength(TParticle const& p);
+ CrossSectionType getCrossSection(TParticle const& p, Code const projectileId,
+ FourMomentum const& projectileP4);
};
} // namespace corsika::proposal
diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp
index c099997a4019de54deb00980ef7349669b744e08..63376ab17e69943c9287cce169a33603fd152c39 100644
--- a/corsika/modules/pythia8/Interaction.hpp
+++ b/corsika/modules/pythia8/Interaction.hpp
@@ -28,23 +28,62 @@ namespace corsika::pythia8 {
void setUnstable(const Code);
void setStable(const Code);
- bool isValidCoMEnergy(HEPEnergyType ecm) { return (10_GeV < ecm) && (ecm < 1_PeV); }
+ bool isValidCoMEnergy(HEPEnergyType const ecm) const {
+ return (10_GeV < ecm) && (ecm < 1_PeV);
+ }
- bool canInteract(const Code);
+ bool canInteract(const Code) const;
void configureLabFrameCollision(const Code, const Code, const HEPEnergyType);
- std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
- const Code BeamId, const Code TargetId, const HEPEnergyType CoMenergy);
+ /**
+ * Returns inelastic AND elastic cross sections.
+ *
+ * These cross sections must correspond to the process described in doInteraction
+ * AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
+ * nuclei or single nucleons (p,n,hydrogen). This "InelEla" method is used since
+ * Sibyll must be useful inside the NuclearInteraction model, which requires that.
+ *
+ * @param projectile is the Code of the projectile
+ * @param target is the Code of the target
+ * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
+ * @param Aprojectil is the mass number of the projectils, if it is a nucleus
+ * @param Atarget is the mass number of the target, if it is a nucleus
+ *
+ * @return a tuple of: inelastic cross section, elastic cross section
+ */
+ std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+ Code const projectile, Code const target, FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const;
- template <typename TParticle>
- GrammageType getInteractionLength(TParticle const&);
+ /**
+ * Returns inelastic (production) cross section.
+ *
+ * This cross section must correspond to the process described in doInteraction.
+ * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+ *
+ * @param projectile is the Code of the projectile
+ * @param target is the Code of the target
+ * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
+ * @param Aprojectil is the mass number of the projectils, if it is a nucleus
+ * @param Atarget is the mass number of the target, if it is a nucleus
+ *
+ * @return inelastic cross section
+ * elastic cross section
+ */
+ CrossSectionType getCrossSection(Code const projectile, Code const target,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+ return std::get<0>(
+ getCrossSectionInelEla(projectile, target, projectileP4, targetP4));
+ }
/**
- In this function PYTHIA is called to produce one event. The
- event is copied (and boosted) into the shower lab frame.
+ * In this function PYTHIA is called to produce one event. The
+ * event is copied (and boosted) into the shower lab frame.
*/
template <typename TView>
- void doInteraction(TView&);
+ void doInteraction(TView& output, Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
private:
default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("pythia");
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index 86734017907bd60e4b1c8013bb3185781aa5b19d..dbf99a5cba85272435cfd0827862981698a23e39 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -317,8 +317,10 @@ int main(int argc, char** argv) {
HEPEnergyType const emcut = 1_GeV;
HEPEnergyType const hadcut = 1_GeV;
ParticleCut cut(emcut, emcut, hadcut, hadcut, true);
+
corsika::proposal::Interaction emCascade(env);
- InteractionCounter emCascadeCounted(emCascade);
+ // NOT available for PROPOSAL due to interface trouble:
+ // InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess emContinuous(env);
BetheBlochPDG emContinuous(showerAxis);
@@ -335,7 +337,7 @@ int main(int argc, char** argv) {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
- bool operator()(const Particle& p) { return (p.getKineticEnergy() < cutE_); }
+ bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
};
auto hadronSequence = make_select(EnergySwitch(63.1_GeV), urqmdCounted, heModelCounted);
auto decaySequence = make_sequence(decayPythia, decaySibyll);
@@ -352,8 +354,8 @@ int main(int argc, char** argv) {
output.add("particles", observationLevel);
// assemble the final process sequence
- auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence,
- emCascadeCounted, cut, emContinuous, // trackWriter,
+ auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, cut,
+ emCascade, emContinuous, // trackWriter,
observationLevel, longprof);
/* === END: SETUP PROCESS LIST === */
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 62aa91e564010d12a0fbcb3ab48062917ca42da4..b01d7b608b85d661d7fdce6004a1ef9efcd3cd4b 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -144,7 +144,9 @@ int main(int argc, char** argv) {
ParticleCut cut(60_GeV, 60_GeV, 100_PeV, 100_PeV, true);
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess emContinuous(env);
- InteractionCounter emCascadeCounted(emCascade);
+
+ // NOT possible right now, due to interface differenc in PROPOSAL
+ // InteractionCounter emCascadeCounted(emCascade);
TrackWriter trackWriter;
output.add("tracks", trackWriter); // register TrackWriter
@@ -157,8 +159,8 @@ int main(int argc, char** argv) {
obsPlane, DirectionVector(rootCS, {1., 0., 0.}), "particles.dat");
output.add("obsplane", observationLevel);
- auto sequence = make_sequence(emCascadeCounted, emContinuous, longprof, cut,
- observationLevel, trackWriter);
+ auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, observationLevel,
+ trackWriter);
// define air shower object, run simulation
setup::Tracking tracking;
Cascade EAS(env, tracking, sequence, output, stack);
@@ -183,9 +185,6 @@ int main(int argc, char** argv) {
cut.reset();
emContinuous.reset();
- auto const hists = emCascadeCounted.getHistogram();
- save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
- save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
longprof.save("longprof_emShower.txt");
output.endOfLibrary();
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 4a85351b0e9e5ebdb4c24910a9fa77655c14215a..0ff10af164ff51fc3db6cde5bbb7324dea410d24 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -248,7 +248,7 @@ int main(int argc, char** argv) {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
- bool operator()(const setup::Stack::particle_type& p) {
+ bool operator()(const setup::Stack::particle_type& p) const {
return (p.getEnergy() < cutE_);
}
};
diff --git a/examples/mars.cpp b/examples/mars.cpp
index 60b4f87d43d29004285c6fafd59168f693d5ca43..a284f1ecd85580ace4ec82867bf5ad8573ea5b05 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -345,7 +345,8 @@ int main(int argc, char** argv) {
HEPEnergyType const hadcut = 1_GeV;
ParticleCut cut(emcut, emcut, hadcut, hadcut, true);
corsika::proposal::Interaction emCascade(env);
- InteractionCounter emCascadeCounted(emCascade);
+ // NOT possible right now, due to interface difference for PROPOSAL:
+ // InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess emContinuous(env);
BetheBlochPDG emContinuous(showerAxis);
@@ -378,8 +379,8 @@ int main(int argc, char** argv) {
// assemble the final process sequence
auto sequence =
- make_sequence(stackInspect, hadronSequence, decaySequence, emCascadeCounted,
- emContinuous, cut, trackWriter, observationLevel, longprof);
+ make_sequence(stackInspect, hadronSequence, decaySequence, emCascade, emContinuous,
+ cut, trackWriter, observationLevel, longprof);
/* === END: SETUP PROCESS LIST === */
// create the cascade object using the default stack and tracking implementation
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 911f45c528f6f440c45b7a94048c2e3478b5a8a3..d6c1af1dd4dada9e68f4adb4d2cec4f92facb42b 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -236,7 +236,7 @@ int main(int argc, char** argv) {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
- bool operator()(const Particle& p) { return (p.getEnergy() < cutE_); }
+ bool operator()(const Particle& p) const { return (p.getEnergy() < cutE_); }
};
auto hadronSequence =
make_select(EnergySwitch(55_GeV), urqmdCounted,
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index 7cf2c55c48fdcc32cb81daa0b65a59e4a4600bc3..9b6dce5056a78abe10529b96f75a32e521b31055 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -50,9 +50,6 @@ struct DummyNode {
int data_ = 0;
};
-// The stack is non-existent for this example
-struct DummyStack {};
-
// our data object (particle) is a simple arrary of doubles
struct DummyData {
double data_[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
@@ -65,6 +62,9 @@ struct DummyData {
HEPEnergyType getEnergy() const { return 10_GeV; }
};
+// The stack is non-existent for this example
+struct DummyStack {};
+
// there is no real trajectory/track
struct DummyTrajectory {};
@@ -75,6 +75,7 @@ struct DummyView {
: p_(p) {}
DummyData& p_;
DummyData& parent() { return p_; }
+ // this is only needed because of PROPOSAL interface right now:
};
int globalCount = 0; // simple counter
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 1f36ee5f4a6d29d3640039cae80206d35979220d..a4e1c1b63c1b86df66e347d3ab28ca294985e97e 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -5,13 +5,13 @@ set (test_modules_sources
## testExecTime.cpp --> needs to be fixed, see #326
testObservationPlane.cpp
testQGSJetII.cpp
- #testPythia8.cpp
+ testPythia8.cpp
testUrQMD.cpp
- #testCONEX.cpp
+ testCONEX.cpp
## testOnShellCheck.cpp
testParticleCut.cpp
testSibyll.cpp
- #testEpos.cpp
+ testEpos.cpp
)
CORSIKA_ADD_TEST (testModules SOURCES ${test_modules_sources})
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index 0e449f16b6a38eaeac30b95a3475c7b302532f76..2ccc5407d7c26e42df552c942b62a8c0f6affd7f 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -94,6 +94,8 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("Pythia8Interface", "modules") {
logging::set_level(logging::level::info);
+
+ auto const rootCS = get_root_CoordinateSystem();
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
auto const& cs = *csPtr;
{
@@ -166,7 +168,6 @@ TEST_CASE("Pythia8Interface", "modules") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Proton, 7_TeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
auto& view = *secViewPtr;
- auto const particle = stackPtr->getNextParticle();
corsika::pythia8::Interaction collision;
@@ -179,50 +180,36 @@ TEST_CASE("Pythia8Interface", "modules") {
CHECK_FALSE(collision.canInteract(Code::Electron));
// nuclei not supported
- CHECK_THROWS(collision.getCrossSection(Code::Proton, Code::Helium, 1_TeV));
std::tuple<CrossSectionType, CrossSectionType> xs_test =
- collision.getCrossSection(Code::Iron, Code::Hydrogen, 1_GeV);
- CHECK(std::get<0>(xs_test) == std::numeric_limits<double>::infinity() * 1_mb);
- CHECK(std::get<1>(xs_test) == std::numeric_limits<double>::infinity() * 1_mb);
-
- collision.getInteractionLength(particle);
-
- collision.doInteraction(view);
- [[maybe_unused]] const GrammageType length = collision.getInteractionLength(particle);
- CHECK(length / 1_kg * square(1_m) == Approx(43.04).margin(5e-1));
+ collision.getCrossSectionInelEla(
+ Code::Proton, Code::Hydrogen,
+ {sqrt(static_pow<2>(Iron::mass) + static_pow<2>(100_GeV)),
+ {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
+ CHECK(std::get<0>(xs_test) / 1_mb == Approx(353).margin(2));
+ CHECK(std::get<1>(xs_test) / 1_mb == Approx(75).margin(2));
+
+ collision.doInteraction(view, Code::Proton, Code::Hydrogen,
+ {sqrt(static_pow<2>(Iron::mass) + static_pow<2>(100_GeV)),
+ {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}});
CHECK(view.getSize() == 38);
}
- SECTION("pythia nucleus projectile") {
-
- // this is a projectile nucleus with very little energy
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Oxygen, 17_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
- auto& view = *secViewPtr;
- auto particle = stackPtr->first();
-
- corsika::pythia8::Interaction collision;
-
- GrammageType lambda_test = collision.getInteractionLength(particle);
- CHECK(lambda_test == std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm));
-
- CHECK_THROWS(collision.doInteraction(view));
- }
-
SECTION("pythia too low energy") {
// this is a projectile neutron with very little energy
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Neutron, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
auto& view = *secViewPtr;
- auto particle = stackPtr->first();
corsika::pythia8::Interaction collision;
- GrammageType lambda_test = collision.getInteractionLength(particle);
- CHECK(lambda_test == std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm));
-
- CHECK_THROWS(collision.doInteraction(view));
+ CHECK_THROWS(collision.doInteraction(
+ view, Code::Neutron, Code::Hydrogen,
+ {sqrt(static_pow<2>(Neutron::mass) + static_pow<2>(100_GeV)),
+ {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
}
SECTION("pythia wrong target") {
@@ -244,6 +231,32 @@ TEST_CASE("Pythia8Interface", "modules") {
corsika::pythia8::Interaction collision;
- CHECK_THROWS(collision.doInteraction(view));
+ CHECK_THROWS(collision.getCrossSection(
+ Code::Proton, Code::Iron,
+ {sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
+ {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {Iron::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
+
+ CHECK_THROWS(collision.doInteraction(
+ view, Code::Proton, Code::Iron,
+ {sqrt(static_pow<2>(Proton::mass) + static_pow<2>(100_GeV)),
+ {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {Iron::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
+ }
+
+ SECTION("pythia wrong projectile") {
+
+ // resonable projectile, but tool low energy
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Iron, 1_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ { [[maybe_unused]] auto const& dummy_StackPtr = stackPtr; }
+
+ corsika::pythia8::Interaction collision;
+
+ CHECK_THROWS(
+ collision.doInteraction(*secViewPtr, Code::Iron, Code::Hydrogen,
+ {sqrt(static_pow<2>(Iron::mass) + static_pow<2>(100_GeV)),
+ {rootCS, {0_eV, 0_eV, 100_GeV}}},
+ {Hydrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}}));
}
}
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 94fb90792bfb3c7184f0db69051b7ad3f7d411bb..105f5b983905137804bc1ec5eabcfc157702406d 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -28,10 +28,6 @@
*/
#include <corsika/modules/sibyll/Random.hpp>
-// TODODODODODODODODOODOD THIS MUST BE REMOVED
-#include <corsika/modules/epos/Random.hpp>
-// TODODODODODODODODOODOD THIS MUST BE REMOVED
-
using namespace corsika;
using namespace corsika::sibyll;
@@ -278,7 +274,7 @@ TEST_CASE("SibyllInterface", "modules") {
// CHECK(view.getSize() == 11);
CHECK(cx / 1_mb == Approx(1100).margin(100)); // this is not physics validation
// CHECK(view.getSize() == 20); // also sibyll not stable wrt. to compiler changes
- CHECK(view.getSize() == Approx(90).margin(10)); // this is not physics validation
+ CHECK(view.getSize() == Approx(25).margin(10)); // this is not physics validation
}
}