IAP GITLAB

Skip to content
Snippets Groups Projects
Commit d82b8ca0 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan: Committed by Maximilian Reininghaus
Browse files

updated to latest coding style

parent 59d25dd5
No related branches found
No related tags found
2 merge requests!234WIP: Initial example of python as script language from C++,!205UrQMD improvements
...@@ -36,18 +36,18 @@ UrQMD::UrQMD(std::string const& xs_file) { ...@@ -36,18 +36,18 @@ UrQMD::UrQMD(std::string const& xs_file) {
iniurqmd_(); iniurqmd_();
} }
CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code vProjectileCode, CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode,
corsika::particles::Code vTargetCode, corsika::particles::Code targetCode,
HEPEnergyType vLabEnergy) const { HEPEnergyType labEnergy) const {
// translated to C++ from CORSIKA 7 subroutine cxtot_u // translated to C++ from CORSIKA 7 subroutine cxtot_u
auto const kinEnergy = vLabEnergy - particles::GetMass(vProjectileCode); auto const kinEnergy = labEnergy - particles::GetMass(projectileCode);
assert(kinEnergy >= HEPEnergyType::zero()); assert(kinEnergy >= HEPEnergyType::zero());
double const logKinEnergy = std::log10(kinEnergy * (1 / 1_GeV)); double const logKinEnergy = std::log10(kinEnergy * (1 / 1_GeV));
double const ye = std::max(10 * logKinEnergy + 10.5, 1.); 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)); int const je = std::min(int(ye), int(xs_interp_support_table_.shape()[2] - 2));
std::array<double, 3> w; std::array<double, 3> w;
w[2 - 1] = ye - je; w[2 - 1] = ye - je;
w[3 - 1] = w[2 - 1] * (w[2 - 1] - 1.) * .5; w[3 - 1] = w[2 - 1] * (w[2 - 1] - 1.) * .5;
...@@ -55,7 +55,7 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code vProjectileCode ...@@ -55,7 +55,7 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code vProjectileCode
w[2 - 1] = w[2 - 1] - 2 * w[3 - 1]; w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
int projectileIndex; int projectileIndex;
switch (vProjectileCode) { switch (projectileCode) {
case particles::Code::Proton: case particles::Code::Proton:
projectileIndex = 0; projectileIndex = 0;
break; break;
...@@ -89,13 +89,13 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code vProjectileCode ...@@ -89,13 +89,13 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code vProjectileCode
projectileIndex = 8; projectileIndex = 8;
break; break;
default: default:
std::cout << "WARNING: UrQMD cross-section not tabulated for " << vProjectileCode std::cout << "WARNING: UrQMD cross-section not tabulated for " << projectileCode
<< std::endl; << std::endl;
return CrossSectionType::zero(); return CrossSectionType::zero();
} }
int targetIndex; int targetIndex;
switch (vTargetCode) { switch (targetCode) {
case particles::Code::Nitrogen: case particles::Code::Nitrogen:
targetIndex = 0; targetIndex = 0;
break; break;
...@@ -107,38 +107,37 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code vProjectileCode ...@@ -107,38 +107,37 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code vProjectileCode
break; break;
default: default:
std::stringstream ss; std::stringstream ss;
ss << "UrQMD cross-section not tabluated for target " << vTargetCode; ss << "UrQMD cross-section not tabluated for target " << targetCode;
throw std::runtime_error(ss.str().data()); throw std::runtime_error(ss.str().data());
} }
auto result = CrossSectionType::zero(); auto result = CrossSectionType::zero();
for (int i = 0; i < 3; ++i) { for (int i = 0; i < 3; ++i) {
result += result +=
xs_interp_support_table[projectileIndex][targetIndex][je + i - 1 - 1] * w[i]; xs_interp_support_table_[projectileIndex][targetIndex][je + i - 1 - 1] * w[i];
} }
return result; return result;
} }
CrossSectionType UrQMD::GetCrossSection(particles::Code vProjectileCode, CrossSectionType UrQMD::GetCrossSection(particles::Code projectileCode,
corsika::particles::Code vTargetCode, corsika::particles::Code targetCode,
HEPEnergyType vLabEnergy, HEPEnergyType labEnergy, int projectileA) const {
int vAProjectile) const {
// the following is a (incomplete!) translation of ptsigtot() into C++ // the following is a (incomplete!) translation of ptsigtot() into C++
if (vProjectileCode != particles::Code::Nucleus && if (projectileCode != particles::Code::Nucleus &&
!IsNucleus(vTargetCode)) { // both particles are "special" !IsNucleus(targetCode)) { // both particles are "special"
auto const mProj = particles::GetMass(vProjectileCode); auto const mProj = particles::GetMass(projectileCode);
auto const mTar = particles::GetMass(vTargetCode); auto const mTar = particles::GetMass(targetCode);
double sqrtS = sqrt(units::si::detail::static_pow<2>(mProj) + double sqrtS = sqrt(units::si::detail::static_pow<2>(mProj) +
units::si::detail::static_pow<2>(mTar) + 2 * vLabEnergy * mTar) * units::si::detail::static_pow<2>(mTar) + 2 * labEnergy * mTar) *
(1 / 1_GeV); (1 / 1_GeV);
// we must set some UrQMD globals first... // we must set some UrQMD globals first...
auto const [ityp, iso3] = ConvertToUrQMD(vProjectileCode); auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
inputs_.spityp[0] = ityp; inputs_.spityp[0] = ityp;
inputs_.spiso3[0] = iso3; inputs_.spiso3[0] = iso3;
auto const [itypTar, iso3Tar] = ConvertToUrQMD(vTargetCode); auto const [itypTar, iso3Tar] = ConvertToUrQMD(targetCode);
inputs_.spityp[1] = itypTar; inputs_.spityp[1] = itypTar;
inputs_.spiso3[1] = iso3Tar; inputs_.spiso3[1] = iso3Tar;
...@@ -178,8 +177,8 @@ CrossSectionType UrQMD::GetCrossSection(particles::Code vProjectileCode, ...@@ -178,8 +177,8 @@ CrossSectionType UrQMD::GetCrossSection(particles::Code vProjectileCode,
return sigEl * 0_mb; return sigEl * 0_mb;
} }
} else { } else {
int const Ap = vAProjectile; int const Ap = projectileA;
int const At = IsNucleus(vTargetCode) ? particles::GetNucleusA(vTargetCode) : 1; int const At = IsNucleus(targetCode) ? particles::GetNucleusA(targetCode) : 1;
double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1]; double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1];
return 10_mb * M_PI * units::si::detail::static_pow<2>(maxImpact); return 10_mb * M_PI * units::si::detail::static_pow<2>(maxImpact);
...@@ -206,7 +205,7 @@ CrossSectionType UrQMD::GetCrossSection(TParticle const& projectile, ...@@ -206,7 +205,7 @@ CrossSectionType UrQMD::GetCrossSection(TParticle const& projectile,
return GetTabulatedCrossSection(projectileCode, targetCode, projectileEnergyLab); return GetTabulatedCrossSection(projectileCode, targetCode, projectileEnergyLab);
} }
bool UrQMD::CanInteract(particles::Code vCode) const { bool UrQMD::CanInteract(particles::Code code) const {
// According to the manual, UrQMD can use all mesons, baryons and nucleons // According to the manual, UrQMD can use all mesons, baryons and nucleons
// which are modeled also as input particles. I think it is safer to accept // which are modeled also as input particles. I think it is safer to accept
// only the usual long-lived species as input. // only the usual long-lived species as input.
...@@ -222,52 +221,52 @@ bool UrQMD::CanInteract(particles::Code vCode) const { ...@@ -222,52 +221,52 @@ bool UrQMD::CanInteract(particles::Code vCode) const {
particles::Code::K0Long}; particles::Code::K0Long};
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes), return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes); code) != std::cend(validProjectileCodes);
} }
GrammageType UrQMD::GetInteractionLength(SetupParticle const& vParticle) const { GrammageType UrQMD::GetInteractionLength(SetupParticle const& particle) const {
if (!CanInteract(vParticle.GetPID())) { if (!CanInteract(particle.GetPID())) {
// we could do the canInteract check in GetCrossSection, too but if // we could do the canInteract check in GetCrossSection, too but if
// we do it here we have the advantage of avoiding the loop // we do it here we have the advantage of avoiding the loop
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm); return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
} }
auto const& mediumComposition = auto const& mediumComposition =
vParticle.GetNode()->GetModelProperties().GetNuclearComposition(); particle.GetNode()->GetModelProperties().GetNuclearComposition();
using namespace std::placeholders; using namespace std::placeholders;
CrossSectionType const weightedProdCrossSection = mediumComposition.WeightedSum( CrossSectionType const weightedProdCrossSection = mediumComposition.WeightedSum(
std::bind(&UrQMD::GetCrossSection<decltype(vParticle)>, this, vParticle, _1)); std::bind(&UrQMD::GetCrossSection<decltype(particle)>, this, particle, _1));
return mediumComposition.GetAverageMassNumber() * units::constants::u / return mediumComposition.GetAverageMassNumber() * units::constants::u /
weightedProdCrossSection; weightedProdCrossSection;
} }
corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjectile) { corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectile) {
using namespace units::si; using namespace units::si;
auto projectileCode = vProjectile.GetPID(); auto projectileCode = projectile.GetPID();
auto const projectileEnergyLab = vProjectile.GetEnergy(); auto const projectileEnergyLab = projectile.GetEnergy();
auto const& projectileMomentumLab = vProjectile.GetMomentum(); auto const& projectileMomentumLab = projectile.GetMomentum();
auto const& projectilePosition = vProjectile.GetPosition(); auto const& projectilePosition = projectile.GetPosition();
auto const projectileTime = vProjectile.GetTime(); auto const projectileTime = projectile.GetTime();
// sample target particle // sample target particle
auto const& mediumComposition = auto const& mediumComposition =
vProjectile.GetNode()->GetModelProperties().GetNuclearComposition(); projectile.GetNode()->GetModelProperties().GetNuclearComposition();
auto const componentCrossSections = std::invoke([&]() { auto const componentCrossSections = std::invoke([&]() {
auto const& components = mediumComposition.GetComponents(); auto const& components = mediumComposition.GetComponents();
std::vector<CrossSectionType> crossSections; std::vector<CrossSectionType> crossSections;
crossSections.reserve(components.size()); crossSections.reserve(components.size());
for (auto const c : components) { for (auto const c : components) {
crossSections.push_back(GetCrossSection(vProjectile, c)); crossSections.push_back(GetCrossSection(projectile, c));
} }
return crossSections; return crossSections;
}); });
auto const targetCode = mediumComposition.SampleTarget(componentCrossSections, fRNG); auto const targetCode = mediumComposition.SampleTarget(componentCrossSections, rng_);
auto const targetA = particles::GetNucleusA(targetCode); auto const targetA = particles::GetNucleusA(targetCode);
auto const targetZ = particles::GetNucleusZ(targetCode); auto const targetZ = particles::GetNucleusZ(targetCode);
...@@ -281,10 +280,10 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti ...@@ -281,10 +280,10 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti
// is this everything? // is this everything?
inputs_.prspflg = 0; inputs_.prspflg = 0;
sys_.Ap = vProjectile.GetNuclearA(); sys_.Ap = projectile.GetNuclearA();
sys_.Zp = vProjectile.GetNuclearZ(); sys_.Zp = projectile.GetNuclearZ();
rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV) / rsys_.ebeam = (projectileEnergyLab - projectile.GetMass()) * (1 / 1_GeV) /
vProjectile.GetNuclearA(); projectile.GetNuclearA();
rsys_.bdist = nucrad_(targetA) + nucrad_(sys_.Ap) + 2 * options_.CTParam[30 - 1]; rsys_.bdist = nucrad_(targetA) + nucrad_(sys_.Ap) + 2 * options_.CTParam[30 - 1];
...@@ -294,11 +293,11 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti ...@@ -294,11 +293,11 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti
inputs_.prspflg = 1; inputs_.prspflg = 1;
sys_.Ap = 1; // even for non-baryons this has to be set, see vanilla UrQMD.f sys_.Ap = 1; // even for non-baryons this has to be set, see vanilla UrQMD.f
rsys_.bdist = nucrad_(targetA) + nucrad_(1) + 2 * options_.CTParam[30 - 1]; rsys_.bdist = nucrad_(targetA) + nucrad_(1) + 2 * options_.CTParam[30 - 1];
rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV); rsys_.ebeam = (projectileEnergyLab - projectile.GetMass()) * (1 / 1_GeV);
if (projectileCode == particles::Code::K0Long || if (projectileCode == particles::Code::K0Long ||
projectileCode == particles::Code::K0Short) { projectileCode == particles::Code::K0Short) {
projectileCode = fBooleanDist(fRNG) ? particles::Code::K0 : particles::Code::K0Bar; projectileCode = booleanDist_(rng_) ? particles::Code::K0 : particles::Code::K0Bar;
} }
auto const [ityp, iso3] = ConvertToUrQMD(projectileCode); auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
...@@ -332,7 +331,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti ...@@ -332,7 +331,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti
for (int i = 0; i < sys_.npart; ++i) { for (int i = 0; i < sys_.npart; ++i) {
auto code = ConvertFromUrQMD(isys_.ityp[i], isys_.iso3[i]); auto code = ConvertFromUrQMD(isys_.ityp[i], isys_.iso3[i]);
if (code == particles::Code::K0 || code == particles::Code::K0Bar) { if (code == particles::Code::K0 || code == particles::Code::K0Bar) {
code = fBooleanDist(fRNG) ? particles::Code::K0Short : particles::Code::K0Long; code = booleanDist_(rng_) ? particles::Code::K0Short : particles::Code::K0Long;
} }
// "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well // "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well
...@@ -346,7 +345,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti ...@@ -346,7 +345,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti
momentum.rebase(originalCS); // transform back into standard lab frame momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << i << " " << code << " " << momentum.GetComponents() << std::endl; std::cout << i << " " << code << " " << momentum.GetComponents() << std::endl;
vProjectile.AddSecondary( projectile.AddSecondary(
std::tuple<particles::Code, HEPEnergyType, stack::MomentumVector, geometry::Point, std::tuple<particles::Code, HEPEnergyType, stack::MomentumVector, geometry::Point,
TimeType>{code, energy, momentum, projectilePosition, projectileTime}); TimeType>{code, energy, momentum, projectilePosition, projectileTime});
} }
...@@ -452,8 +451,8 @@ void UrQMD::readXSFile(std::string const& filename) { ...@@ -452,8 +451,8 @@ void UrQMD::readXSFile(std::string const& filename) {
int nTargets, nProjectiles, nSupports; int nTargets, nProjectiles, nSupports;
ss >> dummy >> nTargets >> nProjectiles >> nSupports; ss >> dummy >> nTargets >> nProjectiles >> nSupports;
decltype(xs_interp_support_table)::extent_gen extents; decltype(xs_interp_support_table_)::extent_gen extents;
xs_interp_support_table.resize(extents[nProjectiles][nTargets][nSupports]); xs_interp_support_table_.resize(extents[nProjectiles][nTargets][nSupports]);
for (int i = 0; i < nTargets; ++i) { for (int i = 0; i < nTargets; ++i) {
for (int j = 0; j < nProjectiles; ++j) { for (int j = 0; j < nProjectiles; ++j) {
...@@ -462,7 +461,7 @@ void UrQMD::readXSFile(std::string const& filename) { ...@@ -462,7 +461,7 @@ void UrQMD::readXSFile(std::string const& filename) {
std::stringstream s(line); std::stringstream s(line);
double energy, sigma; double energy, sigma;
s >> energy >> sigma; s >> energy >> sigma;
xs_interp_support_table[j][i][k] = sigma * 1_mb; xs_interp_support_table_[j][i][k] = sigma * 1_mb;
} }
std::getline(file, line); std::getline(file, line);
......
...@@ -52,14 +52,14 @@ namespace corsika::process::UrQMD { ...@@ -52,14 +52,14 @@ namespace corsika::process::UrQMD {
private: private:
void readXSFile(std::string const&); void readXSFile(std::string const&);
corsika::random::RNG& fRNG = corsika::random::RNG& rng_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD"); corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");
std::uniform_int_distribution<int> fBooleanDist{0, 1}; std::uniform_int_distribution<int> booleanDist_{0, 1};
boost::multi_array<corsika::units::si::CrossSectionType, 3> xs_interp_support_table; boost::multi_array<corsika::units::si::CrossSectionType, 3> xs_interp_support_table_;
}; };
namespace constants { namespace details::constants {
// from coms.f // from coms.f
int constexpr nmax = 500; int constexpr nmax = 500;
...@@ -70,10 +70,10 @@ namespace corsika::process::UrQMD { ...@@ -70,10 +70,10 @@ namespace corsika::process::UrQMD {
// from inputs.f // from inputs.f
int constexpr aamax = 300; int constexpr aamax = 300;
} // namespace constants } // namespace details::constants
template <typename T> template <typename T>
using nmaxArray = std::array<T, constants::nmax>; using nmaxArray = std::array<T, details::constants::nmax>;
using nmaxIntArray = nmaxArray<int>; using nmaxIntArray = nmaxArray<int>;
using nmaxDoubleArray = nmaxArray<double>; using nmaxDoubleArray = nmaxArray<double>;
...@@ -134,8 +134,8 @@ namespace corsika::process::UrQMD { ...@@ -134,8 +134,8 @@ namespace corsika::process::UrQMD {
// defined in options.f // defined in options.f
extern struct { extern struct {
std::array<double, constants::numcto> CTOption; std::array<double, details::constants::numcto> CTOption;
std::array<double, constants::numctp> CTParam; std::array<double, details::constants::numctp> CTParam;
} options_; } options_;
extern struct { extern struct {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment