From 1b90cc27db8a3b3304e020c9f2505eae5ecdebb8 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 6 Apr 2021 15:50:09 +0200
Subject: [PATCH 1/8] UrQMD cross-section read-out from table (copied from
!205)
---
.../detail/framework/utility/CorsikaData.inl | 7 +-
corsika/detail/modules/urqmd/UrQMD.inl | 164 ++++++++++++++++--
corsika/framework/utility/CorsikaData.hpp | 4 +-
corsika/modules/urqmd/UrQMD.hpp | 11 +-
4 files changed, 162 insertions(+), 24 deletions(-)
diff --git a/corsika/detail/framework/utility/CorsikaData.inl b/corsika/detail/framework/utility/CorsikaData.inl
index e2ea1e99b..fc638b81a 100644
--- a/corsika/detail/framework/utility/CorsikaData.inl
+++ b/corsika/detail/framework/utility/CorsikaData.inl
@@ -7,14 +7,15 @@
*/
#pragma once
+#include <boost/filesystem/path.hpp>
+
#include <cstdlib>
#include <stdexcept>
#include <string>
-inline std::string corsika::corsika_data(std::string const& key) {
+inline boost::filesystem::path corsika::corsika_data(boost::filesystem::path const& key) {
if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
- auto const path = std::string(p) + "/" + key;
- return path;
+ return boost::filesystem::path(p) / key;
} else {
throw std::runtime_error("CORSIKA_DATA not set");
}
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index a1abaa6ab..1fc6d285e 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -15,15 +15,112 @@
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <boost/filesystem/path.hpp>
+#include <boost/multi_array.hpp>
+
#include <algorithm>
+#include <cassert>
#include <functional>
#include <iostream>
+#include <fstream>
+#include <sstream>
#include <urqmd.hpp>
namespace corsika::urqmd {
- inline UrQMD::UrQMD() { ::urqmd::iniurqmdc8_(); }
+ inline UrQMD::UrQMD(boost::filesystem::path const& xs_file) {
+ readXSFile(xs_file);
+ ::urqmd::iniurqmdc8_();
+ }
+
+ inline CrossSectionType UrQMD::getTabulatedCrossSection(Code projectileCode,
+ Code targetCode,
+ HEPEnergyType labEnergy) const {
+ // translated to C++ from CORSIKA 7 subroutine cxtot_u
+
+ auto const kinEnergy = labEnergy - get_mass(projectileCode);
+
+ assert(kinEnergy >= HEPEnergyType::zero());
+
+ double const logKinEnergy = std::log10(kinEnergy * (1 / 1_GeV));
+ double const ye = std::max(10 * logKinEnergy + 10.5, 1.);
+ int const je = std::min(int(ye), int(xs_interp_support_table_.shape()[2] - 2));
+ std::array<double, 3> w;
+ w[2 - 1] = ye - je;
+ w[3 - 1] = w[2 - 1] * (w[2 - 1] - 1.) * .5;
+ w[1 - 1] = 1 - w[2 - 1] + w[3 - 1];
+ w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
+
+ int projectileIndex;
+ switch (projectileCode) {
+ case Code::Proton:
+ projectileIndex = 0;
+ break;
+ case Code::AntiProton:
+ projectileIndex = 1;
+ break;
+ case Code::Neutron:
+ projectileIndex = 2;
+ break;
+ case Code::AntiNeutron:
+ projectileIndex = 3;
+ break;
+ case Code::PiPlus:
+ projectileIndex = 4;
+ break;
+ case Code::PiMinus:
+ projectileIndex = 5;
+ break;
+ case Code::KPlus:
+ projectileIndex = 6;
+ break;
+ case Code::KMinus:
+ projectileIndex = 7;
+ break;
+ case Code::K0Short:
+ case Code::K0Long:
+ /* since K0Short and K0Long are treated the same, we can also add K0 and K0Bar
+ * to the list. This is a deviation from CORSIKA 7. */
+ case Code::K0:
+ case Code::K0Bar:
+ projectileIndex = 8;
+ break;
+ default: {
+ CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode);
+ return CrossSectionType::zero();
+ }
+ }
+
+ int targetIndex;
+ switch (targetCode) {
+ case Code::Nitrogen:
+ targetIndex = 0;
+ break;
+ case Code::Oxygen:
+ targetIndex = 1;
+ break;
+ case Code::Argon:
+ targetIndex = 2;
+ break;
+ default:
+ std::stringstream ss;
+ ss << "UrQMD cross-section not tabluated for target " << targetCode;
+ throw std::runtime_error(ss.str().data());
+ }
+
+ auto result = CrossSectionType::zero();
+ for (int i = 0; i < 3; ++i) {
+ result +=
+ xs_interp_support_table_[projectileIndex][targetIndex][je + i - 1 - 1] * w[i];
+ }
+
+ CORSIKA_LOG_TRACE(
+ "UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={} GeV, sigma={}",
+ get_name(projectileCode), get_name(targetCode), labEnergy / 1_GeV, result);
+
+ return result;
+ }
inline CrossSectionType UrQMD::getCrossSection(Code vProjectileCode, Code vTargetCode,
HEPEnergyType vLabEnergy,
@@ -63,20 +160,19 @@ namespace corsika::urqmd {
template <typename TParticle> // need template here, as this is called both with
// SetupParticle as well as SetupProjectile
- inline CrossSectionType UrQMD::getCrossSection(TParticle const& vProjectile,
- Code vTargetCode) const {
- // TODO: return 0 for non-hadrons?
-
- auto const projectileCode = vProjectile.getPID();
- auto const projectileEnergyLab = vProjectile.getEnergy();
-
- if (projectileCode == Code::K0Long) {
- return 0.5 * (getCrossSection(Code::K0, vTargetCode, projectileEnergyLab) +
- getCrossSection(Code::K0Bar, vTargetCode, projectileEnergyLab));
+ inline CrossSectionType UrQMD::getCrossSection(TParticle const& projectile,
+ Code targetCode) const {
+ auto const projectileCode = projectile.getPID();
+
+ if (projectileCode == Code::Nucleus) {
+ /*
+ * unfortunately unavoidable at the moment until we have tools to get the actual
+ * inealstic cross-section from UrQMD
+ */
+ return CrossSectionType::zero();
}
- int const Ap = (projectileCode == Code::Nucleus) ? vProjectile.getNuclearA() : 1;
- return getCrossSection(projectileCode, vTargetCode, projectileEnergyLab, Ap);
+ return getTabulatedCrossSection(projectileCode, targetCode, projectile.getEnergy());
}
inline bool UrQMD::canInteract(Code vCode) const {
@@ -85,10 +181,9 @@ namespace corsika::urqmd {
// only the usual long-lived species as input.
// TODO: Charmed mesons should be added to the list, too
- static Code const validProjectileCodes[] = {
- Code::Nucleus, Code::Proton, Code::AntiProton, Code::Neutron,
- Code::AntiNeutron, Code::PiPlus, Code::PiMinus, Code::KPlus,
- Code::KMinus, Code::K0, Code::K0Bar, Code::K0Long};
+ static std::array constexpr validProjectileCodes{
+ Code::Proton, Code::AntiProton, Code::Neutron, Code::AntiNeutron, Code::PiPlus,
+ Code::PiMinus, Code::KPlus, Code::KMinus, Code::K0Short, Code::K0Long};
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes);
@@ -301,4 +396,39 @@ namespace corsika::urqmd {
return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
}
+ inline void UrQMD::readXSFile(boost::filesystem::path const& filename) {
+ std::ifstream file(filename, std::ios::in);
+
+ if (!file.is_open()) {
+ throw std::runtime_error(filename.native() + " could not be opened.");
+ }
+
+ std::string line;
+
+ std::getline(file, line);
+ std::stringstream ss(line);
+
+ char dummy;
+ int nTargets, nProjectiles, nSupports;
+ ss >> dummy >> nTargets >> nProjectiles >> nSupports;
+
+ decltype(xs_interp_support_table_)::extent_gen extents;
+ xs_interp_support_table_.resize(extents[nProjectiles][nTargets][nSupports]);
+
+ for (int i = 0; i < nTargets; ++i) {
+ for (int j = 0; j < nProjectiles; ++j) {
+ for (int k = 0; k < nSupports; ++k) {
+ std::getline(file, line);
+ std::stringstream s(line);
+ double energy, sigma;
+ s >> energy >> sigma;
+ xs_interp_support_table_[j][i][k] = sigma * 1_mb;
+ }
+
+ std::getline(file, line);
+ std::getline(file, line);
+ }
+ }
+ }
+
} // namespace corsika::urqmd
diff --git a/corsika/framework/utility/CorsikaData.hpp b/corsika/framework/utility/CorsikaData.hpp
index d117db8da..9f24d5036 100644
--- a/corsika/framework/utility/CorsikaData.hpp
+++ b/corsika/framework/utility/CorsikaData.hpp
@@ -8,13 +8,13 @@
#pragma once
-#include <string>
+#include <boost/filesystem/path.hpp>
namespace corsika {
/**
* returns the full path of the file \p filename within the CORSIKA_DATA directory
*/
- std::string corsika_data(std::string const& filename);
+ boost::filesystem::path corsika_data(boost::filesystem::path const& filename);
} // namespace corsika
#include <corsika/detail/framework/utility/CorsikaData.inl>
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
index 8f2fd8ede..9492f49f7 100644
--- a/corsika/modules/urqmd/UrQMD.hpp
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -12,21 +12,27 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
#include <corsika/setup/SetupStack.hpp>
+#include <boost/multi_array.hpp>
+
#include <array>
#include <utility>
+#include <string>
namespace corsika::urqmd {
class UrQMD : public InteractionProcess<UrQMD> {
public:
- UrQMD();
+ UrQMD(boost::filesystem::path const& path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat"));
template <typename TParticle>
GrammageType getInteractionLength(TParticle const&) const;
+ CrossSectionType getTabulatedCrossSection(Code, Code, HEPEnergyType) const;
+
template <typename TParticle>
CrossSectionType getCrossSection(TParticle const&, Code) const;
@@ -39,11 +45,12 @@ namespace corsika::urqmd {
private:
static CrossSectionType getCrossSection(Code, Code, HEPEnergyType, int);
+ void readXSFile(boost::filesystem::path const&);
// data members
default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("urqmd");
-
std::uniform_int_distribution<int> booleanDist_{0, 1};
+ boost::multi_array<CrossSectionType, 3> xs_interp_support_table_;
};
/**
--
GitLab
From 2f82efd94059bce12ada7682061b17aaf197302c Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 6 Apr 2021 15:55:29 +0200
Subject: [PATCH 2/8] missing include path
---
corsika/modules/urqmd/UrQMD.hpp | 1 +
1 file changed, 1 insertion(+)
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
index 9492f49f7..3f6d3e993 100644
--- a/corsika/modules/urqmd/UrQMD.hpp
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -16,6 +16,7 @@
#include <corsika/setup/SetupStack.hpp>
+#include <boost/filesystem/path.hpp>
#include <boost/multi_array.hpp>
#include <array>
--
GitLab
From 4815e08dbc6846d7ff6499967ddc0da5bb78986c Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 6 Apr 2021 16:33:47 +0200
Subject: [PATCH 3/8] inlcude complete boost/filesystem.hpp
---
corsika/detail/modules/urqmd/UrQMD.inl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index 1fc6d285e..a23098a04 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -15,7 +15,7 @@
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
-#include <boost/filesystem/path.hpp>
+#include <boost/filesystem.hpp>
#include <boost/multi_array.hpp>
#include <algorithm>
--
GitLab
From afac3fd7b083470ab59e1cda67ada99b4edad04d Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 6 Apr 2021 16:43:02 +0200
Subject: [PATCH 4/8] boost::filesystem::ifstream
---
corsika/detail/modules/urqmd/UrQMD.inl | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index a23098a04..e85a6abc3 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -397,7 +397,7 @@ namespace corsika::urqmd {
}
inline void UrQMD::readXSFile(boost::filesystem::path const& filename) {
- std::ifstream file(filename, std::ios::in);
+ boost::filesystem::ifstream file(filename, std::ios::in);
if (!file.is_open()) {
throw std::runtime_error(filename.native() + " could not be opened.");
--
GitLab
From 9a85572f6fe8414ff96579fbb7bcf286798958ed Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 6 Apr 2021 17:17:51 +0200
Subject: [PATCH 5/8] recovered testUrQMD from !205
---
tests/modules/testUrQMD.cpp | 20 +++++++++-----------
1 file changed, 9 insertions(+), 11 deletions(-)
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 28ed39a15..92a544932 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -77,21 +77,19 @@ TEST_CASE("UrQMD") {
auto const& cs = *csPtr;
{ [[maybe_unused]] auto const& env_dummy = env; }
- Code validProjectileCodes[] = {Code::PiPlus, Code::PiMinus, Code::Proton,
- Code::Neutron, Code::KPlus, Code::KMinus,
- Code::K0, Code::K0Bar, Code::K0Long};
+ Code validProjectileCodes[] = {Code::PiPlus, Code::PiMinus, Code::Proton,
+ Code::AntiProton, Code::AntiNeutron, Code::Neutron,
+ Code::KPlus, Code::KMinus, Code::K0,
+ Code::K0Bar, Code::K0Long};
for (auto code : validProjectileCodes) {
- auto [stack, view] = setup::testing::setup_stack(
- code, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
- CHECK(stack->getEntries() == 1);
- CHECK(view->getEntries() == 0);
+ auto [stack, view] = setup::testing::setup_stack(code, 0, 0, 100_GeV, nodePtr, cs);
+ REQUIRE(stack->getEntries() == 1);
+ REQUIRE(view->getEntries() == 0);
// simple check whether the cross-section is non-vanishing
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Proton) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Nitrogen) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Oxygen) / 1_mb > 0);
- CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Argon) / 1_mb > 0);
+ // only nuclei with available tabluated data so far
+ REQUIRE(urqmd.getInteractionLength(stack->getNextParticle()) > 1_g / square(1_cm));
}
}
--
GitLab
From 0e1f0607e3878f68932b48861ac50bddf75ba3be Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 6 Apr 2021 17:25:29 +0200
Subject: [PATCH 6/8] removed now unused function
---
modules/urqmd/urqmdInterface.F | 13 -------------
1 file changed, 13 deletions(-)
diff --git a/modules/urqmd/urqmdInterface.F b/modules/urqmd/urqmdInterface.F
index ec2bd092a..446221dfc 100644
--- a/modules/urqmd/urqmdInterface.F
+++ b/modules/urqmd/urqmdInterface.F
@@ -432,16 +432,3 @@ c~ xsegymin=0.25d0
end
-c
-c M. Reininghaus, 2020-04-08
-c
- integer function ReadSigmaLn(ia, ib, ic)
- implicit none
-
- include 'comres.f'
-
- integer :: ia, ib, ic
-
- ReadSigmaLn = SigmaLn(ia, ib, ic)
-
- end function ReadSigmaLn
--
GitLab
From a0534868f30f1a511b6c4a0d613bb1f6b272e03e Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 13 Apr 2021 09:25:28 +0200
Subject: [PATCH 7/8] try to increase test coverage, some documentation
---
.../detail/modules/qgsjetII/Interaction.inl | 14 ++--
corsika/detail/modules/urqmd/UrQMD.inl | 8 ++-
corsika/framework/utility/COMBoost.hpp | 26 ++++++-
corsika/framework/utility/CorsikaData.hpp | 5 ++
corsika/framework/utility/QuarticSolver.hpp | 5 ++
.../framework/utility/SaveBoostHistogram.hpp | 2 +
documentation/index.rst | 1 +
documentation/utilities.rst | 9 +++
tests/modules/testQGSJetII.cpp | 68 ++++++++++++++----
tests/modules/testUrQMD.cpp | 69 ++++++++++++++++++-
10 files changed, 181 insertions(+), 26 deletions(-)
create mode 100644 documentation/utilities.rst
diff --git a/corsika/detail/modules/qgsjetII/Interaction.inl b/corsika/detail/modules/qgsjetII/Interaction.inl
index 828b9fcfe..217fe30b7 100644
--- a/corsika/detail/modules/qgsjetII/Interaction.inl
+++ b/corsika/detail/modules/qgsjetII/Interaction.inl
@@ -239,10 +239,14 @@ namespace corsika::qgsjetII {
if (is_nucleus(targetCode)) { // nucleus
targetMassNumber = get_nucleus_A(targetCode);
if (targetMassNumber > maxMassNumber_)
- throw std::runtime_error("QgsjetII target mass outside range.");
+ throw std::runtime_error(
+ "QgsjetII target mass outside range."); // LCOV_EXCL_LINE there is no
+ // allowed path here
} else {
- if (targetCode != Proton::code)
- throw std::runtime_error("QgsjetII Taget not possible.");
+ if (targetCode != Proton::code) // LCOV_EXCL_LINE there is no allowed path here
+ throw std::runtime_error(
+ "QgsjetII Taget not possible."); // LCOV_EXCL_LINE there is no allowed path
+ // here
}
CORSIKA_LOG_DEBUG("Interaction: target qgsjetII code/A: {}", targetMassNumber);
@@ -252,7 +256,9 @@ namespace corsika::qgsjetII {
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
projectileMassNumber = projectile.getNuclearA();
if (projectileMassNumber > maxMassNumber_)
- throw std::runtime_error("QgsjetII projectile mass outside range.");
+ throw std::runtime_error(
+ "QgsjetII projectile mass outside range."); // LCOV_EXCL_LINE there is no
+ // allowed path here
std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
std::uniform_int_distribution select(0, 1);
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index e85a6abc3..ee29642e0 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -86,9 +86,9 @@ namespace corsika::urqmd {
case Code::K0Bar:
projectileIndex = 8;
break;
- default: {
+ default: { // LCOV_EXCL_START since this can never happen due to canInteract
CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode);
- return CrossSectionType::zero();
+ return CrossSectionType::zero(); // LCOV_EXCL_STOP
}
}
@@ -400,7 +400,9 @@ namespace corsika::urqmd {
boost::filesystem::ifstream file(filename, std::ios::in);
if (!file.is_open()) {
- throw std::runtime_error(filename.native() + " could not be opened.");
+ throw std::runtime_error(
+ filename.native() +
+ " could not be opened."); // LCOV_EXCL_LINE since this is pointless to test
}
std::string line;
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index befd30250..43a123cf1 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -19,14 +19,36 @@
namespace corsika {
/**
+ @defgroup Utilities
+
+ Collection of classes and methods to perform recurring tasks.
+ **/
+
+ /**
+ @class
+ @ingroup Utilities
+
This utility class handles Lorentz boost between different
referenence frames, using FourVector.
+
+ The class is initialized with projectile and optionally target
+ energy/momentum data. During initialization, a rotation matrix is
+ calculated to represent the projectile movement (and thus the
+ boost) along the z-axis. Also the inverse of this rotation is
+ calculated. The Lorentz boost matrix and its inverse are
+ determined as 2x2 matrices considering the energy and
+ pz-momentum.
+
+ Different constructors are offered with different specialization
+ for the cases of collisions (projectile-target) or just decays
+ (projectile only).
*/
class COMBoost {
public:
- //! construct a COMBoost given four-vector of projectile and mass of target
+ //! construct a COMBoost given four-vector of projectile and mass of target (target at
+ //! rest)
COMBoost(FourVector<HEPEnergyType, MomentumVector> const& Pprojectile,
HEPEnergyType const massTarget);
@@ -41,9 +63,11 @@ namespace corsika {
template <typename FourVector>
FourVector fromCoM(FourVector const& p) const;
+ //! returns the rotated coordinate system
CoordinateSystemPtr getRotatedCS() const;
protected:
+ //! internal method
void setBoost(double coshEta, double sinhEta);
private:
diff --git a/corsika/framework/utility/CorsikaData.hpp b/corsika/framework/utility/CorsikaData.hpp
index 9f24d5036..2ab2dec7d 100644
--- a/corsika/framework/utility/CorsikaData.hpp
+++ b/corsika/framework/utility/CorsikaData.hpp
@@ -12,9 +12,14 @@
namespace corsika {
/**
+ * @file CorsikaData.hpp
+ * @ingroup Utilities
+ * @{
* returns the full path of the file \p filename within the CORSIKA_DATA directory
*/
boost::filesystem::path corsika_data(boost::filesystem::path const& filename);
+
+ //! @}
} // namespace corsika
#include <corsika/detail/framework/utility/CorsikaData.inl>
diff --git a/corsika/framework/utility/QuarticSolver.hpp b/corsika/framework/utility/QuarticSolver.hpp
index 0bcdfcedf..ad9484939 100644
--- a/corsika/framework/utility/QuarticSolver.hpp
+++ b/corsika/framework/utility/QuarticSolver.hpp
@@ -11,6 +11,9 @@
#include <complex>
/**
+ * @file QuarticSolver.hpp
+ * @ingroup Utilities
+ * @{
* \todo convert to class
*/
@@ -56,6 +59,8 @@ namespace corsika::quartic_solver {
// afterwards.
DComplex* solve_quartic(double a, double b, double c, double d);
+ //! @}
+
} // namespace corsika::quartic_solver
#include <corsika/detail/framework/utility/QuarticSolver.inl>
diff --git a/corsika/framework/utility/SaveBoostHistogram.hpp b/corsika/framework/utility/SaveBoostHistogram.hpp
index 14e037996..30cbdda2a 100644
--- a/corsika/framework/utility/SaveBoostHistogram.hpp
+++ b/corsika/framework/utility/SaveBoostHistogram.hpp
@@ -13,6 +13,8 @@
namespace corsika {
/**
+ * @ingroup Utilities
+ *
* This functions saves a boost::histogram into a numpy file. Only rather basic axis
* types are supported: regular, variable, integer, category<int>. Only "ordinary" bin
* counts (i.e. a double or int) are supported, nothing fancy like profiles.
diff --git a/documentation/index.rst b/documentation/index.rst
index c1be80a37..79440c503 100644
--- a/documentation/index.rst
+++ b/documentation/index.rst
@@ -13,6 +13,7 @@ Welcome to the CORSIKA 8 air shower simulation framework.
units
environment
stack
+ utilities
api
diff --git a/documentation/utilities.rst b/documentation/utilities.rst
new file mode 100644
index 000000000..37c289fde
--- /dev/null
+++ b/documentation/utilities.rst
@@ -0,0 +1,9 @@
+Utilities
+==========
+
+.. doxygengroup:: Utilities
+ :project: CORSIKA8
+ :members:
+
+
+
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 477dedc6c..4a604fd8b 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -175,7 +175,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
SECTION("InteractionInterface Nuclei") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Nucleus, 20, 10, 10100_GeV,
+ Code::Nucleus, 60, 30, 20100_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
@@ -183,10 +183,15 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto const projectileMomentum = projectile.getMomentum();
corsika::qgsjetII::Interaction model;
- model.doInteraction(view);
+ model.doInteraction(view); // this also should produce some fragments
+ CHECK(view.getSize() == Approx(188).margin(2)); // this is not physics validation
+ int countFragments = 0;
+ for (auto const& sec : view) { countFragments += (sec.getPID() == Code::Nucleus); }
+ CHECK(countFragments == Approx(2).margin(1)); // this is not physics validation
[[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
- CHECK(length / (1_g / square(1_cm)) == Approx(20.13).margin(0.1));
+ CHECK(length / (1_g / square(1_cm)) ==
+ Approx(12).margin(2)); // this is not physics validation
}
SECTION("Heavy nuclei") {
@@ -210,17 +215,50 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
}
SECTION("Allowed Particles") {
-
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Electron, 0, 0, 1100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
- auto particle = stackPtr->first();
- auto projectile = secViewPtr->getProjectile();
- auto const projectileMomentum = projectile.getMomentum();
-
- corsika::qgsjetII::Interaction model;
-
- GrammageType const length = model.getInteractionLength(particle);
- CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
+ { // electron
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Electron, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ auto particle = stackPtr->first();
+ corsika::qgsjetII::Interaction model;
+ GrammageType const length = model.getInteractionLength(particle);
+ CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
+ }
+ { // pi0 is internally converted into pi+/pi-
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Pi0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(18).margin(2)); // this is not physics validation
+ }
+ { // rho0 is internally converted into pi-/pi+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Rho0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(7).margin(2)); // this is not physics validation
+ }
+ { // Lambda is internally converted into neutron
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Lambda0, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+ }
+ { // AntiLambda is internally converted into anti neutron
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Lambda0Bar, 0, 0, 100_GeV,
+ (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ setup::StackView& view = *(secViewPtr.get());
+ corsika::qgsjetII::Interaction model;
+ model.doInteraction(view);
+ CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+ }
}
}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 92a544932..da788b797 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -84,15 +84,34 @@ TEST_CASE("UrQMD") {
for (auto code : validProjectileCodes) {
auto [stack, view] = setup::testing::setup_stack(code, 0, 0, 100_GeV, nodePtr, cs);
- REQUIRE(stack->getEntries() == 1);
- REQUIRE(view->getEntries() == 0);
+ CHECK(stack->getEntries() == 1);
+ CHECK(view->getEntries() == 0);
// simple check whether the cross-section is non-vanishing
// only nuclei with available tabluated data so far
- REQUIRE(urqmd.getInteractionLength(stack->getNextParticle()) > 1_g / square(1_cm));
+ CHECK(urqmd.getInteractionLength(stack->getNextParticle()) > 1_g / square(1_cm));
}
}
+ SECTION("targets options") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
+ auto const& cs = *csPtr;
+ { [[maybe_unused]] auto const& env_dummy = env; }
+ auto [stack, view] =
+ setup::testing::setup_stack(Code::Proton, 0, 0, 100_GeV, nodePtr, cs);
+ CHECK(urqmd.getInteractionLength(stack->getNextParticle()) / 1_g * square(1_cm) ==
+ Approx(105).margin(5));
+ }
+
+ SECTION("invalid targets options") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Omega);
+ auto const& cs = *csPtr;
+ { [[maybe_unused]] auto const& env_dummy = env; }
+ auto [stack, view] =
+ setup::testing::setup_stack(Code::Neutron, 0, 0, 100_GeV, nodePtr, cs);
+ CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle()));
+ }
+
SECTION("nucleus projectile") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
@@ -144,6 +163,32 @@ TEST_CASE("UrQMD") {
Approx(0).margin(1e-2));
}
+ SECTION("\"special\" projectile and target") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ CHECK(stackPtr->getEntries() == 1);
+ CHECK(secViewPtr->getEntries() == 0);
+
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
+
+ urqmd.doInteraction(*secViewPtr);
+
+ CHECK(sumCharge(*secViewPtr) ==
+ get_charge_number(Code::PiPlus) + get_charge_number(Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
+ CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
+ Approx(0).margin(1e-2));
+ }
+
SECTION("K0Long projectile") {
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
@@ -169,4 +214,22 @@ TEST_CASE("UrQMD") {
CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
Approx(0).margin(1e-2));
}
+
+ SECTION("K0Short projectile") {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::K0Short, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ CHECK(stackPtr->getEntries() == 1);
+ CHECK(secViewPtr->getEntries() == 0);
+
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
+
+ CHECK_THROWS(urqmd.doInteraction(*secViewPtr));
+ }
}
--
GitLab
From 4f7e2bb824539390cd2ed82191bd4726734a6d78 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 13 Apr 2021 09:48:46 +0200
Subject: [PATCH 8/8] fixed one unit test
---
tests/modules/testUrQMD.cpp | 49 +++++++++++++++++++++++--------------
1 file changed, 31 insertions(+), 18 deletions(-)
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index da788b797..0f9f4a862 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -164,29 +164,42 @@ TEST_CASE("UrQMD") {
}
SECTION("\"special\" projectile and target") {
- auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
- [[maybe_unused]] auto const& env_dummy = env; // against warnings
- [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+ {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); // Code::Proton not a valid target
+ }
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
- CHECK(stackPtr->getEntries() == 1);
- CHECK(secViewPtr->getEntries() == 0);
+ {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
- // must be assigned to variable, cannot be used as rvalue?!
- auto projectile = secViewPtr->getProjectile();
- auto const projectileMomentum = projectile.getMomentum();
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
+ CHECK(stackPtr->getEntries() == 1);
+ CHECK(secViewPtr->getEntries() == 0);
- urqmd.doInteraction(*secViewPtr);
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
- CHECK(sumCharge(*secViewPtr) ==
- get_charge_number(Code::PiPlus) + get_charge_number(Code::Oxygen));
+ urqmd.doInteraction(*secViewPtr);
- auto const secMomSum =
- sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
- CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
- Approx(0).margin(1e-2));
+ CHECK(sumCharge(*secViewPtr) ==
+ get_charge_number(Code::PiPlus) + get_charge_number(Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
+ CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
+ Approx(0).margin(1e-2));
+ }
}
SECTION("K0Long projectile") {
--
GitLab