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] 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