From 8a9ba7673fba68877a8c4a0c4cf5ab2367e9af4c Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Mon, 20 Apr 2020 15:58:17 +0200
Subject: [PATCH] read UrQMD cross-section file
---
Processes/UrQMD/UrQMD.cc | 127 ++++++++++++++++++++++++++++++++++++++-
Processes/UrQMD/UrQMD.h | 13 +++-
2 files changed, 137 insertions(+), 3 deletions(-)
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index c91dca3e5..7225a9551 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -15,19 +15,107 @@
#include <corsika/units/PhysicalUnits.h>
#include <algorithm>
+#include <array>
+#include <cassert>
+#include <cmath>
+#include <fstream>
#include <functional>
#include <iostream>
#include <random>
+#include <sstream>
using namespace corsika::process::UrQMD;
using namespace corsika::units::si;
-UrQMD::UrQMD() { iniurqmd_(); }
-
using SetupStack = corsika::setup::Stack;
using SetupParticle = corsika::setup::Stack::StackIterator;
using SetupProjectile = corsika::setup::StackView::StackIterator;
+UrQMD::UrQMD(std::filesystem::path const& xs_file) {
+ readXSFile(xs_file);
+ iniurqmd_();
+}
+
+CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code vProjectileCode,
+ corsika::particles::Code vTargetCode,
+ HEPEnergyType vLabEnergy) const {
+ // translated to C++ from CORSIKA 7 subroutine cxtot_u
+
+ auto const kinEnergy = vLabEnergy - particles::GetMass(vProjectileCode);
+
+ 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 (vProjectileCode) {
+ case particles::Code::Proton:
+ projectileIndex = 0;
+ break;
+ case particles::Code::AntiProton:
+ projectileIndex = 1;
+ break;
+ case particles::Code::Neutron:
+ projectileIndex = 2;
+ break;
+ case particles::Code::AntiNeutron:
+ projectileIndex = 3;
+ break;
+ case particles::Code::PiPlus:
+ projectileIndex = 4;
+ break;
+ case particles::Code::PiMinus:
+ projectileIndex = 5;
+ break;
+ case particles::Code::KPlus:
+ projectileIndex = 6;
+ break;
+ case particles::Code::KMinus:
+ projectileIndex = 7;
+ break;
+ case particles::Code::K0Short:
+ case particles::Code::K0Long:
+ projectileIndex = 8;
+ break;
+ default:
+ std::cout << "WARNING: UrQMD cross-section not tabulated for " << vProjectileCode
+ << std::endl;
+ return CrossSectionType::zero();
+ }
+
+ int targetIndex;
+ switch (vTargetCode) {
+ case particles::Code::Nitrogen:
+ targetIndex = 0;
+ break;
+ case particles::Code::Oxygen:
+ targetIndex = 1;
+ break;
+ case particles::Code::Argon:
+ targetIndex = 2;
+ break;
+ default:
+ std::stringstream ss;
+ ss << "UrQMD cross-section not tabluated for target " << vTargetCode;
+ 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];
+ }
+
+ return result;
+}
+
CrossSectionType UrQMD::GetCrossSection(particles::Code vProjectileCode,
corsika::particles::Code vTargetCode,
HEPEnergyType vLabEnergy,
@@ -344,3 +432,38 @@ std::pair<int, int> corsika::process::UrQMD::ConvertToUrQMD(
return mapPDGToUrQMD.at(static_cast<int>(GetPDG(code)));
}
+
+void UrQMD::readXSFile(std::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);
+ }
+ }
+}
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index d74fed7b8..893fe3697 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -16,15 +16,20 @@
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaData.h>
+
+#include <boost/multi_array.hpp>
#include <array>
+#include <filesystem>
#include <random>
#include <utility>
namespace corsika::process::UrQMD {
class UrQMD : public corsika::process::InteractionProcess<UrQMD> {
public:
- UrQMD();
+ UrQMD(
+ std::filesystem::path const& path = utl::CorsikaData("UrQMD/UrQMD-1.3.1-xs.dat"));
void Init() {}
corsika::units::si::GrammageType GetInteractionLength(
corsika::setup::Stack::StackIterator const&) const;
@@ -37,16 +42,22 @@ namespace corsika::process::UrQMD {
particles::Code, particles::Code, corsika::units::si::HEPEnergyType,
int Ap = 1) const;
+ corsika::units::si::CrossSectionType GetTabulatedCrossSection(
+ particles::Code, particles::Code, corsika::units::si::HEPEnergyType) const;
+
corsika::process::EProcessReturn DoInteraction(
corsika::setup::StackView::StackIterator&);
bool CanInteract(particles::Code) const;
private:
+ void readXSFile(std::filesystem::path const&);
+
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");
std::uniform_int_distribution<int> fBooleanDist{0, 1};
+ boost::multi_array<corsika::units::si::CrossSectionType, 3> xs_interp_support_table;
};
namespace constants {
--
GitLab