From 27967bbeb5bd393e0da02d42cceee203922d624b Mon Sep 17 00:00:00 2001
From: Felix Riehn <friehn@lip.pt>
Date: Fri, 9 Aug 2024 19:14:44 +0200
Subject: [PATCH] switch cross section application to getCrossSectionInelEla,
add corr. routines for Sibyll and QgsJet (INCOMPLETE! returning zero for
elastic cross section)
---
applications/cross_section.cpp | 4 ++--
.../modules/qgsjetII/InteractionModel.inl | 6 ++++++
.../detail/modules/sibyll/InteractionModel.inl | 11 +++++++++++
corsika/modules/qgsjetII/InteractionModel.hpp | 18 ++++++++++++++++++
corsika/modules/sibyll/InteractionModel.hpp | 3 +++
5 files changed, 40 insertions(+), 2 deletions(-)
diff --git a/applications/cross_section.cpp b/applications/cross_section.cpp
index 1d9c55298..19661f1f9 100644
--- a/applications/cross_section.cpp
+++ b/applications/cross_section.cpp
@@ -99,8 +99,8 @@ void calculate_cross_sections(TModel& model, std::string const& model_name,
}
// p-p
- CrossSectionType const xs_prod_pp =
- model->getCrossSection(Code::Proton, Code::Proton, p4protonProj, p4protonTarg);
+ auto const [xs_prod_pp, xs_ela_pp ]=
+ model->getCrossSectionInelEla(Code::Proton, Code::Proton, p4protonProj, p4protonTarg);
// p-O
CrossSectionType const xs_prod_pO =
model->getCrossSection(Code::Proton, Code::Oxygen, p4protonProj, p4oxygenTarg);
diff --git a/corsika/detail/modules/qgsjetII/InteractionModel.inl b/corsika/detail/modules/qgsjetII/InteractionModel.inl
index bb5585af3..32e8398b7 100644
--- a/corsika/detail/modules/qgsjetII/InteractionModel.inl
+++ b/corsika/detail/modules/qgsjetII/InteractionModel.inl
@@ -102,6 +102,12 @@ namespace corsika::qgsjetII {
return sigProd * 1_mb;
}
+inline std::tuple<CrossSectionType, CrossSectionType> InteractionModel::getCrossSectionInelEla(
+ Code projCode, Code targetCode, FourMomentum const& proj4mom,
+ FourMomentum const& target4mom) const {
+ return { getCrossSection(projCode, targetCode, proj4mom, target4mom), CrossSectionType::zero()};
+ }
+
template <typename TSecondaries>
inline void InteractionModel::doInteraction(TSecondaries& view, Code const projectileId,
Code const targetId,
diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
index 5c2be73da..8a3b58f7b 100644
--- a/corsika/detail/modules/sibyll/InteractionModel.inl
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -50,6 +50,17 @@ namespace corsika::sibyll {
target4mom);
}
+inline std::tuple<CrossSectionType, CrossSectionType> InteractionModel::getCrossSectionInelEla(
+ Code projCode, Code targetCode, FourMomentum const& proj4mom,
+ FourMomentum const& target4mom) const {
+ if (is_nucleus(projCode))
+ return { getNuclearInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
+ target4mom), CrossSectionType::zero()};
+ else
+ return getHadronInteractionModel().getCrossSectionInelEla(projCode, targetCode, proj4mom,
+ target4mom);
+ }
+
template <typename TSecondaries>
inline void InteractionModel::doInteraction(TSecondaries& view, Code projCode,
Code targetCode,
diff --git a/corsika/modules/qgsjetII/InteractionModel.hpp b/corsika/modules/qgsjetII/InteractionModel.hpp
index c6fd881f9..bbd3b24e7 100644
--- a/corsika/modules/qgsjetII/InteractionModel.hpp
+++ b/corsika/modules/qgsjetII/InteractionModel.hpp
@@ -55,6 +55,24 @@ namespace corsika::qgsjetII {
FourMomentum const& projectileP4,
FourMomentum const& targetP4) const;
+
+/**
+ * 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).
+ *
+ * @param projectile is the Code of the projectile
+ * @param target is the Code of the target
+ * @param projectileP4: four-momentum of projectile
+ * @param targetP4: four-momentum of target
+ *
+ * @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;
+
/**
* In this function QGSJETII is called to produce one event.
*
diff --git a/corsika/modules/sibyll/InteractionModel.hpp b/corsika/modules/sibyll/InteractionModel.hpp
index c5fc50087..05135bdd6 100644
--- a/corsika/modules/sibyll/InteractionModel.hpp
+++ b/corsika/modules/sibyll/InteractionModel.hpp
@@ -31,6 +31,9 @@ namespace corsika::sibyll {
CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
FourMomentum const&) const;
+ std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(Code, Code, FourMomentum const&,
+ FourMomentum const&) const;
+
template <typename TSecondaries>
void doInteraction(TSecondaries&, Code, Code, FourMomentum const&,
FourMomentum const&);
--
GitLab