From 3c70ce61d834850d17eaece36c58fba24e069bae Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 21 Dec 2020 22:55:50 +0100
Subject: [PATCH] onshellcheck
---
.gitlab-ci.yml | 17 +-
.../framework/utility/CorsikaFenvFallback.inl | 4 +-
corsika/detail/modules/ObservationPlane.inl | 2 +-
corsika/detail/modules/pythia8/Decay.inl | 249 +++++++++++++-----
.../detail/modules/pythia8/Interaction.inl | 243 ++++++++---------
corsika/detail/modules/pythia8/Random.inl | 2 +-
corsika/detail/modules/urqmd/UrQMD.inl | 156 ++++++-----
corsika/modules/pythia8/Decay.hpp | 51 +++-
corsika/modules/pythia8/Interaction.hpp | 40 ++-
corsika/modules/pythia8/Random.hpp | 5 +-
corsika/modules/urqmd/UrQMD.hpp | 28 +-
tests/modules/CMakeLists.txt | 4 +-
tests/modules/testOnShellCheck.cpp | 63 ++---
tests/modules/testPythia8.cpp | 108 ++++----
tests/modules/testUrQMD.cpp | 138 +++++-----
15 files changed, 606 insertions(+), 504 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 4d1503e37..d483397b8 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -67,8 +67,18 @@ check-clang-format:
- if: $CI_COMMIT_BRANCH
allow_failure: true
+### CodeQuality tool ####
+include:
+ - template: Code-Quality.gitlab-ci.yml
-
+code_quality:
+ stage: quality
+ rules:
+ - if: '$CODE_QUALITY_DISABLED'
+ when: never
+ - if: '$CI_PIPELINE_SOURCE == "merge_request_event"' # Run code quality job in merge request pipelines
+ - if: '$CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH' # Run code quality job in pipelines on the master branch (but not in other branch pipelines)
+ - if: '$CI_COMMIT_TAG' # Run code quality job in pipelines for tags
####### CONFIG ##############
@@ -87,6 +97,11 @@ check-clang-format:
- if: $CI_MERGE_REQUEST_ID
- if: $CI_COMMIT_TAG
- if: $CI_COMMIT_BRANCH
+ artifacts:
+ when: on_failure
+ expire_in: 3 days
+ paths:
+ - ${CI_PROJECT_DIR}/build/CMakeFiles/CMakeOutput.log
cache:
paths:
- ${CI_PROJECT_DIR}/build/
diff --git a/corsika/detail/framework/utility/CorsikaFenvFallback.inl b/corsika/detail/framework/utility/CorsikaFenvFallback.inl
index ba0dd888a..40a1079ec 100644
--- a/corsika/detail/framework/utility/CorsikaFenvFallback.inl
+++ b/corsika/detail/framework/utility/CorsikaFenvFallback.inl
@@ -12,7 +12,7 @@
extern "C" {
#warning No enabling/disabling of floating point exceptions - platform needs better implementation
-int feenableexcept(int excepts) { return -1; }
+inline int feenableexcept(int excepts) { return -1; }
-int fedisableexcept(int excepts) { return -1; }
+inline int fedisableexcept(int excepts) { return -1; }
}
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index 42dd6ad5f..68ef5b49a 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -66,7 +66,7 @@ namespace corsika::observation_plane {
auto const pointOfIntersection = trajectory.getPosition(timeOfIntersection);
auto dist = (trajectory.getStartPoint() - pointOfIntersection).getNorm() * 1.0001;
- CORSIKA_LOG_TRACE("ObservationPlane::MaxStepLength l={} m", dist / 1_m);
+ CORSIKA_LOG_TRACE("ObservationPlane::getMaxStepLength l={} m", dist / 1_m);
return dist;
}
diff --git a/corsika/detail/modules/pythia8/Decay.inl b/corsika/detail/modules/pythia8/Decay.inl
index 2e639d2b8..224b48025 100644
--- a/corsika/detail/modules/pythia8/Decay.inl
+++ b/corsika/detail/modules/pythia8/Decay.inl
@@ -10,127 +10,238 @@
#include <corsika/modules/pythia8/Decay.hpp>
#include <corsika/modules/pythia8/Random.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
namespace corsika::pythia8 {
- Decay::Decay(std::vector<corsika::Code> pParticles)
- : fTrackedParticles(pParticles) {}
-
- Decay::~Decay() { std::cout << "Pythia::Decay n=" << fCount << std::endl; }
-
- void Decay::Init() {
-
- Decay::SetParticleListStable(fTrackedParticles);
+ Decay::Decay(const bool print_listing)
+ : print_listing_(print_listing) {
// set random number generator in pythia
Pythia8::RndmEngine* rndm = new corsika::pythia8::Random();
- fPythia.setRndmEnginePtr(rndm);
-
- fPythia.readString("Next:numberShowInfo = 0");
- fPythia.readString("Next:numberShowProcess = 0");
- fPythia.readString("Next:numberShowEvent = 0");
+ pythia_.setRndmEnginePtr(rndm);
+
+ /*
+ issue xyz: definition of particles and decay channels use the same mechanism in
+ corsika and pythia we should force pythia to use the file in corsika.
+ */
+ // bool ParticleData::reInit(string startFile, bool xmlFormat = true)
+ // read in particle data from Corsika 8
+ // pythia_.particleData.reInit("/home/felix/ngcorsika/corsika-build/include/corsika/particles/ParticleData.xml");
+ // pythia_.particleData.checkTable();
+
+ pythia_.readString("Next:numberShowInfo = 0");
+ pythia_.readString("Next:numberShowProcess = 0");
+ pythia_.readString("Next:numberShowEvent = 0");
+
+ pythia_.readString("Print:quiet = on");
+ pythia_.readString("Check:particleData = 0");
+
+ /*
+ switching off event check in pythia is needed to allow decays that are off-shell
+ according to the mass definition in pythia.
+ the consistency of particle masses between event generators is an unsolved issues
+ */
+ std::cout << "Pythia::Init: switching off event checking in pythia.." << std::endl;
+ pythia_.readString("Check:event = 1");
+
+ pythia_.readString("ProcessLevel:all = off");
+ pythia_.readString("ProcessLevel:resonanceDecays = off");
+
+ // making sure
+ setStable(Code::Pi0);
+
+ // pythia_.particleData.readString("59:m0 = 101.00");
+
+ if (!pythia_.init())
+ throw std::runtime_error("Pythia::Decay: Initialization failed!");
+ }
- fPythia.readString("Print:quiet = on");
+ Decay::Decay(std::set<Code> const& those)
+ : handleAllDecays_(false)
+ , handledDecays_(those) {}
+
+ Decay::~Decay() { std::cout << "Pythia::Decay n=" << count_ << std::endl; }
+
+ bool Decay::canHandleDecay(Code const vParticleCode) {
+ // if known to pythia and not proton, electron or neutrino it can decay
+ if (vParticleCode == Code::Proton || vParticleCode == Code::AntiProton ||
+ vParticleCode == Code::NuE || vParticleCode == Code::NuMu ||
+ vParticleCode == Code::NuTau || vParticleCode == Code::NuEBar ||
+ vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
+ vParticleCode == Code::Electron || vParticleCode == Code::Positron)
+ return false;
+ else if (canDecay(vParticleCode)) // non-zero for particles known to sibyll
+ return true;
+ else
+ return false;
+ }
- fPythia.readString("ProcessLevel:all = off");
- fPythia.readString("ProcessLevel:resonanceDecays = off");
+ void Decay::setHandleDecay(Code const vParticleCode) {
+ handleAllDecays_ = false;
+ std::cout << "Pythia::Decay: set to handle decay of " << vParticleCode << std::endl;
+ if (Decay::canHandleDecay(vParticleCode))
+ handledDecays_.insert(vParticleCode);
+ else
+ throw std::runtime_error("this decay can not be handled by pythia!");
+ }
- fPythia.particleData.readString("59:m0 = 101.00");
+ void Decay::setHandleDecay(std::vector<Code> const& vParticleList) {
+ handleAllDecays_ = false;
+ for (auto p : vParticleList) setHandleDecay(p);
+ }
- fPythia.init();
+ bool Decay::isDecayHandled(Code const vParticleCode) {
+ if (handleAllDecays_ && canHandleDecay(vParticleCode))
+ return true;
+ else
+ return handledDecays_.find(vParticleCode) != Decay::handledDecays_.end();
}
- void Decay::SetParticleListStable(const std::vector<corsika::Code> particleList) {
- for (auto p : particleList) Decay::SetStable(p);
+ void Decay::setStable(std::vector<Code> const& particleList) {
+ for (auto p : particleList) Decay::setStable(p);
}
- void Decay::SetUnstable(const corsika::Code pCode) {
+ void Decay::setUnstable(Code const pCode) {
std::cout << "Pythia::Decay: setting " << pCode << " unstable.." << std::endl;
- fPythia.particleData.mayDecay(static_cast<int>(corsika::get_PDG(pCode)), true);
+ pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
}
- void Decay::SetStable(const corsika::Code pCode) {
+ void Decay::setStable(Code const pCode) {
std::cout << "Pythia::Decay: setting " << pCode << " stable.." << std::endl;
- fPythia.particleData.mayDecay(static_cast<int>(corsika::get_PDG(pCode)), false);
+ pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
}
- template <typename TParticle>
- TimeType Decay::GetLifetime(TParticle const& p) {
+ bool Decay::isStable(Code const vCode) {
+ return pythia_.particleData.canDecay(static_cast<int>(get_PDG(vCode)));
+ }
- HEPEnergyType E = p.GetEnergy();
- HEPMassType m = p.GetMass();
+ bool Decay::canDecay(Code const pCode) {
+ std::cout << "Pythia::Decay: checking if particle: " << pCode
+ << " can decay in PYTHIA? ";
+ const bool ans = pythia_.particleData.canDecay(static_cast<int>(get_PDG(pCode)));
+ std::cout << ans << std::endl;
+ return ans;
+ }
- const double gamma = E / m;
+ void Decay::printDecayConfig(const Code vCode) {
+ std::cout << "Decay: Pythia decay configuration:" << std::endl;
+ std::cout << vCode << " is ";
+ if (isStable(vCode))
+ std::cout << "stable" << std::endl;
+ else
+ std::cout << "unstable" << std::endl;
+ }
- const TimeType t0 = corsika::get_lifetime(p.GetPID());
- auto const lifetime = gamma * t0;
+ void Decay::printDecayConfig() {
+ std::cout << "Pythia::Decay: decay configuration:" << std::endl;
+ if (handleAllDecays_)
+ std::cout << " all particles known to Pythia are handled by Pythia::Decay!"
+ << std::endl;
+ else
+ for (auto& pCode : handledDecays_)
+ std::cout << "Decay of " << pCode << " is handled by Pythia!" << std::endl;
+ }
- return lifetime;
+ template <typename TParticle>
+ TimeType Decay::getLifetime(TParticle const& particle) {
+
+ const auto pid = particle.getPID();
+ if (canDecay(pid)) {
+ HEPEnergyType E = particle.getEnergy();
+ HEPMassType m = particle.getMass();
+
+ const double gamma = E / m;
+
+ const TimeType t0 = get_lifetime(pid);
+ auto const lifetime = gamma * t0;
+ std::cout << "Pythia::Decay: code: " << particle.getPID() << std::endl;
+ std::cout << "Pythia::Decay: MinStep: t0: " << t0 << std::endl;
+ std::cout << "Pythia::Decay: MinStep: energy: " << E / 1_GeV << " GeV" << std::endl;
+ std::cout << "Pythia::Decay: momentum: "
+ << particle.getMomentum().getComponents() / 1_GeV << " GeV" << std::endl;
+ std::cout << "Pythia::Decay: MinStep: gamma: " << gamma << std::endl;
+ std::cout << "Pythia::Decay: MinStep: tau: " << lifetime << std::endl;
+
+ return lifetime;
+ } else
+ return std::numeric_limits<double>::infinity() * 1_s;
}
- template <typename TProjectile>
- void Decay::DoDecay(TProjectile& vP) {
- using corsika::Point;
+ template <typename TView>
+ void Decay::doDecay(TView& view) {
+
+ auto projectile = view.getProjectile();
- auto const decayPoint = vP.GetPosition();
- auto const t0 = vP.GetTime();
+ auto const& decayPoint = projectile.getPosition();
+ auto const t0 = projectile.getTime();
- // coordinate system, get global frame of reference
- corsika::CoordinateSystem& rootCS =
- corsika::RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ auto const& labMomentum = projectile.getMomentum();
+ CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
- fCount++;
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Pythia projectile: +z
+ COMBoost const boost(labMomentum, projectile.getMass());
+ auto const& rotatedCS = boost.getRotatedCS();
+
+ count_++;
// pythia stack
- Pythia8::Event& event = fPythia.event;
+ Pythia8::Event& event = pythia_.event;
event.reset();
+ auto const particleId = projectile.getPID();
+
// set particle unstable
- Decay::SetUnstable(vP.GetPID());
+ Decay::setUnstable(particleId);
// input particle PDG
- auto const pdgCode = static_cast<int>(corsika::get_PDG(vP.GetPID()));
+ auto const pdgCode = static_cast<int>(get_PDG(particleId));
- auto const pcomp = vP.GetMomentum().GetComponents();
- double px = pcomp[0] / 1_GeV;
- double py = pcomp[1] / 1_GeV;
- double pz = pcomp[2] / 1_GeV;
- double en = vP.GetEnergy() / 1_GeV;
- double m = corsika::get_mass(vP.GetPID()) / 1_GeV;
+ double constexpr px = 0;
+ double constexpr py = 0;
+ double constexpr pz = 0;
+ double const en = projectile.getMass() / 1_GeV;
+ double const m = en;
// add particle to pythia stack
event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
- if (!fPythia.next())
- std::cout << "Pythia::Decay: decay failed!" << std::endl;
+ if (!pythia_.next())
+ throw std::runtime_error("Pythia::Decay: decay failed!");
else
std::cout << "Pythia::Decay: particles after decay: " << event.size() << std::endl;
- // list final state
- event.list();
+ if (print_listing_) {
+ // list final state
+ event.list();
+ }
// loop over final state
for (int i = 0; i < event.size(); ++i)
if (event[i].isFinal()) {
- auto const pyId =
- corsika::convert_from_PDG(static_cast<corsika::PDGCode>(event[i].id()));
- HEPEnergyType pyEn = event[i].e() * 1_GeV;
- MomentumVector pyP(rootCS, {event[i].px() * 1_GeV, event[i].py() * 1_GeV,
- event[i].pz() * 1_GeV});
-
- std::cout << "particle: id=" << pyId
- << " momentum=" << pyP.GetComponents() / 1_GeV << " energy=" << pyEn
- << std::endl;
-
- vP.AddSecondary(
- std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
- corsika::Point, TimeType>{pyId, pyEn, pyP, decayPoint, t0});
+ auto const pyId = convert_from_PDG(static_cast<PDGCode>(event[i].id()));
+ HEPEnergyType const Erest = event[i].e() * 1_GeV;
+ MomentumVector const pRest(
+ rotatedCS,
+ {event[i].px() * 1_GeV, event[i].py() * 1_GeV, event[i].pz() * 1_GeV});
+ FourVector const fourMomRest{Erest, pRest};
+ auto const fourMomLab = boost.fromCoM(fourMomRest);
+
+ std::cout << "particle: id=" << pyId << " momentum="
+ << fourMomLab.getSpaceLikeComponents().getComponents(labCS) / 1_GeV
+ << " energy=" << fourMomLab.getTimeLikeComponent() << std::endl;
+
+ view.addSecondary(std::make_tuple(pyId, fourMomLab.getTimeLikeComponent(),
+ fourMomLab.getSpaceLikeComponents(), decayPoint,
+ t0));
}
// set particle stable
- Decay::SetStable(vP.GetPID());
+ Decay::setStable(particleId);
}
-
} // namespace corsika::pythia8
diff --git a/corsika/detail/modules/pythia8/Interaction.inl b/corsika/detail/modules/pythia8/Interaction.inl
index 687f40c37..e67d140a9 100644
--- a/corsika/detail/modules/pythia8/Interaction.inl
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -18,126 +18,116 @@
#include <tuple>
-using Projectile = corsika::setup::StackView::ParticleType;
-using Particle = corsika::setup::Stack::ParticleType;
+using Projectile = corsika::setup::StackView::particle_type;
+using Particle = corsika::setup::Stack::particle_type;
namespace corsika::pythia8 {
- typedef corsika::Vector<corsika::hepmomentum_d> MomentumVector;
-
Interaction::~Interaction() {
- std::cout << "Pythia::Interaction n=" << fCount << std::endl;
+ std::cout << "Pythia::Interaction n=" << count_ << std::endl;
}
- void Interaction::Init() {
+ Interaction::Interaction(const bool print_listing)
+ : print_listing_(print_listing) {
- using corsika::RNGManager;
+ std::cout << "Pythia::Interaction n=" << count_ << std::endl;
// initialize Pythia
- if (!fInitialized) {
+ if (!initialized_) {
- fPythia.readString("Print:quiet = on");
+ pythia_.readString("Print:quiet = on");
// TODO: proper process initialization for MinBias needed
- fPythia.readString("HardQCD:all = on");
- fPythia.readString("ProcessLevel:resonanceDecays = off");
+ pythia_.readString("HardQCD:all = on");
+ pythia_.readString("ProcessLevel:resonanceDecays = off");
- fPythia.init();
+ pythia_.init();
// any decays in pythia? if yes need to define which particles
- if (fInternalDecays) {
+ if (internalDecays_) {
// define which particles are passed to corsika, i.e. which particles make it into
// history even very shortlived particles like charm or pi0 are of interest here
- const std::vector<corsika::Code> HadronsWeWantTrackedByCorsika = {
- corsika::Code::PiPlus, corsika::Code::PiMinus,
- corsika::Code::Pi0, corsika::Code::KMinus,
- corsika::Code::KPlus, corsika::Code::K0Long,
- corsika::Code::K0Short, corsika::Code::SigmaPlus,
- corsika::Code::SigmaMinus, corsika::Code::Lambda0,
- corsika::Code::Xi0, corsika::Code::XiMinus,
- corsika::Code::OmegaMinus, corsika::Code::DPlus,
- corsika::Code::DMinus, corsika::Code::D0,
- corsika::Code::D0Bar};
-
- Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
+ const std::vector<Code> HadronsWeWantTrackedByCorsika = {
+ Code::PiPlus, Code::PiMinus, Code::Pi0, Code::KMinus,
+ Code::KPlus, Code::K0Long, Code::K0Short, Code::SigmaPlus,
+ Code::SigmaMinus, Code::Lambda0, Code::Xi0, Code::XiMinus,
+ Code::OmegaMinus, Code::DPlus, Code::DMinus, Code::D0,
+ Code::D0Bar};
+
+ Interaction::setStable(HadronsWeWantTrackedByCorsika);
}
// basic initialization of cross section routines
- fSigma.init(&fPythia.info, fPythia.settings, &fPythia.particleData, &fPythia.rndm);
+ sigma_.init(&pythia_.info, pythia_.settings, &pythia_.particleData, &pythia_.rndm);
- fInitialized = true;
+ initialized_ = true;
}
}
- void Interaction::SetParticleListStable(
- std::vector<corsika::Code> const& particleList) {
- for (auto p : particleList) Interaction::SetStable(p);
+ void Interaction::setStable(std::vector<Code> const& particleList) {
+ for (auto p : particleList) Interaction::setStable(p);
}
- void Interaction::SetUnstable(const corsika::Code pCode) {
+ void Interaction::setUnstable(Code const pCode) {
std::cout << "Pythia::Interaction: setting " << pCode << " unstable.." << std::endl;
- fPythia.particleData.mayDecay(static_cast<int>(corsika::get_PDG(pCode)), true);
+ pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), true);
}
- void Interaction::SetStable(const corsika::Code pCode) {
+ void Interaction::setStable(Code const pCode) {
std::cout << "Pythia::Interaction: setting " << pCode << " stable.." << std::endl;
- fPythia.particleData.mayDecay(static_cast<int>(corsika::get_PDG(pCode)), false);
+ pythia_.particleData.mayDecay(static_cast<int>(get_PDG(pCode)), false);
}
- void Interaction::ConfigureLabFrameCollision(const corsika::Code BeamId,
- const corsika::Code TargetId,
+ void Interaction::configureLabFrameCollision(const Code BeamId, const Code TargetId,
const HEPEnergyType BeamEnergy) {
// Pythia configuration of the current event
// very clumsy. I am sure this can be done better..
// set beam
// beam id for pythia
- auto const pdgBeam = static_cast<int>(corsika::get_PDG(BeamId));
+ auto const pdgBeam = static_cast<int>(get_PDG(BeamId));
std::stringstream stBeam;
stBeam << "Beams:idA = " << pdgBeam;
- fPythia.readString(stBeam.str());
+ pythia_.readString(stBeam.str());
// set target
- auto pdgTarget = static_cast<int>(corsika::get_PDG(TargetId));
+ auto pdgTarget = static_cast<int>(get_PDG(TargetId));
// replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
- if (TargetId == corsika::Code::Hydrogen)
- pdgTarget = static_cast<int>(corsika::get_PDG(corsika::Code::Proton));
+ if (TargetId == Code::Hydrogen) pdgTarget = static_cast<int>(get_PDG(Code::Proton));
std::stringstream stTarget;
stTarget << "Beams:idB = " << pdgTarget;
- fPythia.readString(stTarget.str());
+ pythia_.readString(stTarget.str());
// set frame to lab. frame
- fPythia.readString("Beams:frameType = 2");
+ pythia_.readString("Beams:frameType = 2");
// set beam energy
const double Elab = BeamEnergy / 1_GeV;
std::stringstream stEnergy;
stEnergy << "Beams:eA = " << Elab;
- fPythia.readString(stEnergy.str());
+ pythia_.readString(stEnergy.str());
// target at rest
- fPythia.readString("Beams:eB = 0.");
+ pythia_.readString("Beams:eB = 0.");
// initialize this config
- fPythia.init();
+ pythia_.init();
}
- bool Interaction::CanInteract(const corsika::Code pCode) {
- return pCode == corsika::Code::Proton || pCode == corsika::Code::Neutron ||
- pCode == corsika::Code::AntiProton || pCode == corsika::Code::AntiNeutron ||
- pCode == corsika::Code::PiMinus || pCode == corsika::Code::PiPlus;
+ bool Interaction::canInteract(Code const pCode) {
+ return pCode == Code::Proton || pCode == Code::Neutron || pCode == Code::AntiProton ||
+ pCode == Code::AntiNeutron || pCode == Code::PiMinus || pCode == Code::PiPlus;
}
- std::tuple<CrossSectionType, CrossSectionType> Interaction::GetCrossSection(
- const corsika::Code BeamId, const corsika::Code TargetId,
- const HEPEnergyType CoMenergy) {
+ std::tuple<CrossSectionType, CrossSectionType> Interaction::getCrossSection(
+ const Code BeamId, const Code TargetId, const HEPEnergyType CoMenergy) {
// interaction possible in pythia?
- if (TargetId == corsika::Code::Proton || TargetId == corsika::Code::Hydrogen) {
- if (CanInteract(BeamId) && ValidCoMEnergy(CoMenergy)) {
+ if (TargetId == Code::Proton || TargetId == Code::Hydrogen) {
+ if (canInteract(BeamId) && isValidCoMEnergy(CoMenergy)) {
// input particle PDG
- auto const pdgCodeBeam = static_cast<int>(corsika::get_PDG(BeamId));
- auto const pdgCodeTarget = static_cast<int>(corsika::get_PDG(TargetId));
+ auto const pdgCodeBeam = static_cast<int>(get_PDG(BeamId));
+ auto const pdgCodeTarget = static_cast<int>(get_PDG(TargetId));
const double ecm = CoMenergy / 1_GeV;
// calculate cross section
- fSigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
- if (fSigma.hasSigmaTot()) {
- const double sigEla = fSigma.sigmaEl();
- const double sigProd = fSigma.sigmaTot() - sigEla;
+ sigma_.calc(pdgCodeBeam, pdgCodeTarget, ecm);
+ if (sigma_.hasSigmaTot()) {
+ const double sigEla = sigma_.sigmaEl();
+ const double sigProd = sigma_.sigmaTot() - sigEla;
return std::make_tuple(sigProd * (1_fm * 1_fm), sigEla * (1_fm * 1_fm));
@@ -154,38 +144,38 @@ namespace corsika::pythia8 {
}
template <>
- GrammageType Interaction::GetInteractionLength(Particle& p) {
+ GrammageType Interaction::getInteractionLength(Particle& particle) {
// coordinate system, get global frame of reference
- CoordinateSystem& rootCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ MomentumVector const& pMomentum = particle.getMomentum();
+ CoordinateSystemPtr const& labCS = pMomentum.getCoordinateSystem();
- const corsika::Code corsikaBeamId = p.GetPID();
+ const Code corsikaBeamId = particle.getPID();
// beam particles for pythia : 1, 2, 3 for p, pi, k
// read from cross section code table
- const bool kInteraction = CanInteract(corsikaBeamId);
+ const bool kInteraction = canInteract(corsikaBeamId);
// FOR NOW: assume target is at rest
- corsika::MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
+ MomentumVector pTarget(labCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
- HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
- corsika::MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
- pTotLab += p.GetMomentum();
+ HEPEnergyType Elab = particle.getEnergy() + constants::nucleonMass;
+ MomentumVector pTotLab(labCS, {0_GeV, 0_GeV, 0_GeV});
+ pTotLab += pMomentum;
pTotLab += pTarget;
- auto const pTotLabNorm = pTotLab.norm();
+ auto const pTotLabNorm = pTotLab.getNorm();
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
std::cout << "Interaction: LambdaInt: \n"
- << " input energy: " << p.GetEnergy() / 1_GeV << std::endl
+ << " input energy: " << particle.getEnergy() / 1_GeV << std::endl
<< " beam can interact:" << kInteraction << std::endl
- << " beam pid:" << p.GetPID() << std::endl;
+ << " beam pid:" << particle.getPID() << std::endl;
// TODO: move limits into variables
- if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM)) {
+ if (kInteraction && Elab >= 8.5_GeV && isValidCoMEnergy(ECoM)) {
// get target from environment
/*
@@ -193,23 +183,23 @@ namespace corsika::pythia8 {
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
- const auto* currentNode = p.GetNode();
+ const auto* currentNode = particle.getNode();
const auto mediumComposition =
- currentNode->GetModelProperties().getNuclearComposition();
+ currentNode->getModelProperties().getNuclearComposition();
// determine average interaction length
auto const weightedProdCrossSection =
- mediumComposition.WeightedSum([=](auto vTargetID) {
- return std::get<0>(this->GetCrossSection(corsikaBeamId, vTargetID, ECoM));
+ mediumComposition.getWeightedSum([=](auto vTargetID) {
+ return std::get<0>(this->getCrossSection(corsikaBeamId, vTargetID, ECoM));
});
std::cout << "Interaction: IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mb << std::endl
<< "Interaction: IntLength: average mass number: "
- << mediumComposition.GetAverageMassNumber() << std::endl;
+ << mediumComposition.getAverageMassNumber() << std::endl;
// calculate interaction length in medium
- GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
+ GrammageType const int_length = mediumComposition.getAverageMassNumber() *
constants::u / weightedProdCrossSection;
std::cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
@@ -226,43 +216,43 @@ namespace corsika::pythia8 {
event is copied (and boosted) into the shower lab frame.
*/
- template <>
- void Interaction::doInteraction(Projectile& vP) {
+ template <class TView>
+ void Interaction::doInteraction(TView& view) {
+
+ auto projectile = view.getProjectile();
- const auto corsikaBeamId = vP.GetPID();
+ const auto corsikaBeamId = projectile.getPID();
std::cout << "Pythia::Interaction: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
- << corsika::pythia8::Interaction::CanInteract(corsikaBeamId) << std::endl;
+ << corsika::pythia8::Interaction::canInteract(corsikaBeamId) << std::endl;
- if (corsika::is_nucleus(corsikaBeamId)) {
+ if (is_nucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
}
- if (corsika::pythia8::Interaction::CanInteract(corsikaBeamId)) {
+ if (corsika::pythia8::Interaction::canInteract(corsikaBeamId)) {
- const CoordinateSystem& rootCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
+ // define projectile
+ HEPEnergyType const eProjectileLab = projectile.getEnergy();
+ auto const pProjectileLab = projectile.getMomentum();
+ CoordinateSystemPtr const& labCS = pProjectileLab.getCoordinateSystem();
// position and time of interaction, not used in Sibyll
- Point pOrig = vP.GetPosition();
- TimeType tOrig = vP.GetTime();
+ Point pOrig = projectile.getPosition();
+ TimeType tOrig = projectile.getTime();
// define target
// FOR NOW: target is always at rest
const auto eTargetLab = 0_GeV + constants::nucleonMass;
- const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
+ const auto pTargetLab = MomentumVector(labCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
- // define projectile
- HEPEnergyType const eProjectileLab = vP.GetEnergy();
- auto const pProjectileLab = vP.GetMomentum();
-
std::cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << std::endl
- << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
+ << "Interaction: pbeam lab: " << pProjectileLab.getComponents() / 1_GeV
<< std::endl;
std::cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << std::endl
- << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV
+ << "Interaction: ptarget lab: " << pTargetLab.getComponents() / 1_GeV
<< std::endl;
const FourVector PprojLab(eProjectileLab, pProjectileLab);
@@ -279,101 +269,98 @@ namespace corsika::pythia8 {
// boost target
auto const PtargCoM = boost.toCoM(PtargLab);
- std::cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
+ std::cout << "Interaction: ebeam CoM: " << PprojCoM.getTimeLikeComponent() / 1_GeV
<< std::endl
<< "Interaction: pbeam CoM: "
- << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << std::endl;
- std::cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
+ << PprojCoM.getSpaceLikeComponents().getComponents() / 1_GeV << std::endl;
+ std::cout << "Interaction: etarget CoM: " << PtargCoM.getTimeLikeComponent() / 1_GeV
<< std::endl
<< "Interaction: ptarget CoM: "
- << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << std::endl;
+ << PtargCoM.getSpaceLikeComponents().getComponents() / 1_GeV << std::endl;
- std::cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
+ std::cout << "Interaction: position of interaction: " << pOrig.getCoordinates()
<< std::endl;
std::cout << "Interaction: time: " << tOrig << std::endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = vP.GetMomentum();
+ MomentumVector Ptot = projectile.getMomentum();
// invariant mass, i.e. cm. energy
- HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
+ HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.getSquaredNorm());
// sample target mass number
- const auto* currentNode = vP.GetNode();
+ const auto* currentNode = projectile.getNode();
const auto& mediumComposition =
- currentNode->GetModelProperties().getNuclearComposition();
+ currentNode->getModelProperties().getNuclearComposition();
// get cross sections for target materials
/*
Here we read the cross section from the interaction model again,
- should be passed from GetInteractionLength if possible
+ should be passed from getInteractionLength if possible
*/
//#warning reading interaction cross section again, should not be necessary
- auto const& compVec = mediumComposition.GetComponents();
+ auto const& compVec = mediumComposition.getComponents();
std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ const auto [sigProd, sigEla] = getCrossSection(corsikaBeamId, targetId, Ecm);
[[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
}
const auto corsikaTargetId =
- mediumComposition.SampleTarget(cross_section_of_components, fRNG);
+ mediumComposition.sampleTarget(cross_section_of_components, RNG_);
std::cout << "Interaction: target selected: " << corsikaTargetId << std::endl;
- if (corsikaTargetId != corsika::Code::Hydrogen &&
- corsikaTargetId != corsika::Code::Neutron &&
- corsikaTargetId != corsika::Code::Proton)
+ if (corsikaTargetId != Code::Hydrogen && corsikaTargetId != Code::Neutron &&
+ corsikaTargetId != Code::Proton)
throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
std::cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
- if (eProjectileLab < 8.5_GeV || !ValidCoMEnergy(Ecm)) {
+ if (eProjectileLab < 8.5_GeV || !isValidCoMEnergy(Ecm)) {
std::cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << std::endl;
throw std::runtime_error("energy too low for PYTHIA");
} else {
- fCount++;
+ count_++;
- ConfigureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
+ configureLabFrameCollision(corsikaBeamId, corsikaTargetId, eProjectileLab);
// create event in pytia
- if (!fPythia.next()) throw std::runtime_error("Pythia::DoInteraction: failed!");
+ if (!pythia_.next()) throw std::runtime_error("Pythia::DoInteraction: failed!");
// link to pythia stack
- Pythia8::Event& event = fPythia.event;
+ Pythia8::Event& event = pythia_.event;
// print final state
event.list();
- MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ MomentumVector Plab_final(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
for (int i = 0; i < event.size(); ++i) {
Pythia8::Particle& p8p = event[i];
// skip particles that have decayed in pythia
if (!p8p.isFinal()) continue;
- auto const pyId =
- corsika::convert_from_PDG(static_cast<corsika::PDGCode>(p8p.id()));
+ auto const pyId = convert_from_PDG(static_cast<PDGCode>(p8p.id()));
const MomentumVector pyPlab(
- rootCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
+ labCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
HEPEnergyType const pyEn = p8p.e() * 1_GeV;
// add to corsika stack
- auto pnew = vP.AddSecondary(
- std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
- corsika::Point, TimeType>{pyId, pyEn, pyPlab, pOrig, tOrig});
+ auto pnew =
+ projectile.addSecondary(std::make_tuple(pyId, pyEn, pyPlab, pOrig, tOrig));
- Plab_final += pnew.GetMomentum();
- Elab_final += pnew.GetEnergy();
+ Plab_final += pnew.getMomentum();
+ Elab_final += pnew.getEnergy();
}
std::cout << "conservation (all GeV): "
<< "Elab_final=" << Elab_final / 1_GeV
- << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << std::endl;
+ << ", Plab_final=" << (Plab_final / 1_GeV).getComponents() << std::endl;
}
}
}
diff --git a/corsika/detail/modules/pythia8/Random.inl b/corsika/detail/modules/pythia8/Random.inl
index 980052613..c077e0bf7 100644
--- a/corsika/detail/modules/pythia8/Random.inl
+++ b/corsika/detail/modules/pythia8/Random.inl
@@ -12,6 +12,6 @@
namespace corsika::pythia8 {
- double Random::flat() { return fDist(fRNG); }
+ double Random::flat() { return Dist_(RNG_); }
} // namespace corsika::pythia8
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index f76772168..9909a6062 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -25,30 +25,25 @@ namespace corsika::urqmd {
UrQMD::UrQMD() { ::urqmd::iniurqmdc8_(); }
- using SetupStack = corsika::setup::Stack;
- using SetupParticle = corsika::setup::Stack::StackIterator;
- using SetupProjectile = corsika::setup::StackView::StackIterator;
-
- CrossSectionType UrQMD::GetCrossSection(corsika::Code vProjectileCode,
- corsika::Code vTargetCode,
+ CrossSectionType UrQMD::getCrossSection(Code vProjectileCode, Code vTargetCode,
HEPEnergyType vLabEnergy,
int vAProjectile = 1) {
// the following is a translation of ptsigtot() into C++
- if (vProjectileCode != corsika::Code::Nucleus &&
- !corsika::is_nucleus(vTargetCode)) { // both particles are "special"
- auto const mProj = corsika::get_mass(vProjectileCode);
- auto const mTar = corsika::get_mass(vTargetCode);
+ if (vProjectileCode != Code::Nucleus &&
+ !is_nucleus(vTargetCode)) { // both particles are "special"
+ auto const mProj = get_mass(vProjectileCode);
+ auto const mTar = get_mass(vTargetCode);
double sqrtS =
sqrt(static_pow<2>(mProj) + static_pow<2>(mTar) + 2 * vLabEnergy * mTar) *
(1 / 1_GeV);
// we must set some UrQMD globals first...
- auto const [ityp, iso3] = ConvertToUrQMD(vProjectileCode);
+ auto const [ityp, iso3] = convertToUrQMD(vProjectileCode);
::urqmd::inputs_.spityp[0] = ityp;
::urqmd::inputs_.spiso3[0] = iso3;
- auto const [itypTar, iso3Tar] = ConvertToUrQMD(vTargetCode);
+ auto const [itypTar, iso3Tar] = convertToUrQMD(vTargetCode);
::urqmd::inputs_.spityp[1] = itypTar;
::urqmd::inputs_.spiso3[1] = iso3Tar;
@@ -57,7 +52,7 @@ namespace corsika::urqmd {
return ::urqmd::sigtot_(one, two, sqrtS) * 1_mb;
} else {
int const Ap = vAProjectile;
- int const At = is_nucleus(vTargetCode) ? corsika::get_nucleus_A(vTargetCode) : 1;
+ int const At = is_nucleus(vTargetCode) ? get_nucleus_A(vTargetCode) : 1;
double const maxImpact = ::urqmd::nucrad_(Ap) + ::urqmd::nucrad_(At) +
2 * ::urqmd::options_.CTParam[30 - 1];
@@ -68,85 +63,86 @@ namespace corsika::urqmd {
template <typename TParticle> // need template here, as this is called both with
// SetupParticle as well as SetupProjectile
- CrossSectionType UrQMD::GetCrossSection(TParticle const& vProjectile,
- corsika::Code vTargetCode) const {
+ 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();
+ auto const projectileCode = vProjectile.getPID();
+ auto const projectileEnergyLab = vProjectile.getEnergy();
- if (projectileCode == corsika::Code::K0Long) {
- return 0.5 *
- (GetCrossSection(corsika::Code::K0, vTargetCode, projectileEnergyLab) +
- GetCrossSection(corsika::Code::K0Bar, vTargetCode, projectileEnergyLab));
+ if (projectileCode == Code::K0Long) {
+ return 0.5 * (getCrossSection(Code::K0, vTargetCode, projectileEnergyLab) +
+ getCrossSection(Code::K0Bar, vTargetCode, projectileEnergyLab));
}
- int const Ap =
- (projectileCode == corsika::Code::Nucleus) ? vProjectile.GetNuclearA() : 1;
- return GetCrossSection(projectileCode, vTargetCode, projectileEnergyLab, Ap);
+ int const Ap = (projectileCode == Code::Nucleus) ? vProjectile.getNuclearA() : 1;
+ return getCrossSection(projectileCode, vTargetCode, projectileEnergyLab, Ap);
}
- bool UrQMD::CanInteract(corsika::Code vCode) const {
+ bool UrQMD::canInteract(Code vCode) const {
// 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
// only the usual long-lived species as input.
// TODO: Charmed mesons should be added to the list, too
- static corsika::Code const validProjectileCodes[] = {
- corsika::Code::Nucleus, corsika::Code::Proton, corsika::Code::AntiProton,
- corsika::Code::Neutron, corsika::Code::AntiNeutron, corsika::Code::PiPlus,
- corsika::Code::PiMinus, corsika::Code::KPlus, corsika::Code::KMinus,
- corsika::Code::K0, corsika::Code::K0Bar, corsika::Code::K0Long};
+ 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};
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes);
}
- GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle) const {
+ template <typename TParticle>
+ GrammageType UrQMD::getInteractionLength(TParticle const& vParticle) const {
- if (!CanInteract(vParticle.GetPID())) {
- // we could do the canInteract check in GetCrossSection, too but if
+ if (!canInteract(vParticle.getPID())) {
+ // we could do the canInteract check in getCrossSection, too but if
// we do it here we have the advantage of avoiding the loop
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
auto const& mediumComposition =
- vParticle.GetNode()->GetModelProperties().getNuclearComposition();
+ vParticle.getNode()->getModelProperties().getNuclearComposition();
using namespace std::placeholders;
CrossSectionType const weightedProdCrossSection = mediumComposition.WeightedSum(
- std::bind(&UrQMD::GetCrossSection<decltype(vParticle)>, this, vParticle, _1));
+ std::bind(&UrQMD::getCrossSection<decltype(vParticle)>, this, vParticle, _1));
- return mediumComposition.GetAverageMassNumber() * constants::u /
+ return mediumComposition.getAverageMassNumber() * constants::u /
weightedProdCrossSection;
}
- void UrQMD::doInteraction(SetupProjectile& vProjectile) {
+ template <typename TView>
+ void UrQMD::doInteraction(TView& view) {
+
+ auto projectile = view.getProjectile();
- auto projectileCode = vProjectile.GetPID();
- auto const projectileEnergyLab = vProjectile.GetEnergy();
- auto const& projectileMomentumLab = vProjectile.GetMomentum();
- auto const& projectilePosition = vProjectile.GetPosition();
- auto const projectileTime = vProjectile.GetTime();
+ auto projectileCode = projectile.getPID();
+ auto const projectileEnergyLab = projectile.getEnergy();
+ auto const& projectileMomentumLab = projectile.getMomentum();
+ auto const& projectilePosition = projectile.getPosition();
+ auto const projectileTime = projectile.getTime();
// sample target particle
auto const& mediumComposition =
- vProjectile.GetNode()->GetModelProperties().getNuclearComposition();
+ projectile.getNode()->getModelProperties().getNuclearComposition();
auto const componentCrossSections = std::invoke([&]() {
- auto const& components = mediumComposition.GetComponents();
+ auto const& components = mediumComposition.getComponents();
std::vector<CrossSectionType> crossSections;
crossSections.reserve(components.size());
for (auto const c : components) {
- crossSections.push_back(GetCrossSection(vProjectile, c));
+ crossSections.push_back(getCrossSection(projectile, c));
}
return crossSections;
});
- auto const targetCode = mediumComposition.SampleTarget(componentCrossSections, fRNG);
- auto const targetA = corsika::get_nucleus_A(targetCode);
- auto const targetZ = corsika::get_nucleus_Z(targetCode);
+ auto const targetCode = mediumComposition.sampleTarget(componentCrossSections, RNG_);
+ auto const targetA = get_nucleus_A(targetCode);
+ auto const targetZ = get_nucleus_Z(targetCode);
::urqmd::inputs_.nevents = 1;
::urqmd::sys_.eos = 0; // could be configurable in principle
@@ -154,14 +150,14 @@ namespace corsika::urqmd {
::urqmd::sys_.nsteps = 1;
// initialization regarding projectile
- if (corsika::Code::Nucleus == projectileCode) {
+ if (Code::Nucleus == projectileCode) {
// is this everything?
::urqmd::inputs_.prspflg = 0;
- ::urqmd::sys_.Ap = vProjectile.GetNuclearA();
- ::urqmd::sys_.Zp = vProjectile.GetNuclearZ();
- ::urqmd::rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV) /
- vProjectile.GetNuclearA();
+ ::urqmd::sys_.Ap = projectile.getNuclearA();
+ ::urqmd::sys_.Zp = projectile.getNuclearZ();
+ ::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV) /
+ projectile.getNuclearA();
::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) +
::urqmd::nucrad_(::urqmd::sys_.Ap) +
@@ -175,22 +171,22 @@ namespace corsika::urqmd {
1; // even for non-baryons this has to be set, see vanilla UrQMD.f
::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) + ::urqmd::nucrad_(1) +
2 * ::urqmd::options_.CTParam[30 - 1];
- ::urqmd::rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV);
+ ::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV);
- if (projectileCode == corsika::Code::K0Long) {
- projectileCode = fBooleanDist(fRNG) ? corsika::Code::K0 : corsika::Code::K0Bar;
- } else if (projectileCode == corsika::Code::K0Short) {
+ if (projectileCode == Code::K0Long) {
+ projectileCode = booleanDist_(RNG_) ? Code::K0 : Code::K0Bar;
+ } else if (projectileCode == Code::K0Short) {
throw std::runtime_error("K0Short should not interact");
}
- auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
+ auto const [ityp, iso3] = convertToUrQMD(projectileCode);
// todo: conversion of K_long/short into strong eigenstates;
::urqmd::inputs_.spityp[0] = ityp;
::urqmd::inputs_.spiso3[0] = iso3;
}
// initilazation regarding target
- if (corsika::is_nucleus(targetCode)) {
+ if (is_nucleus(targetCode)) {
::urqmd::sys_.Zt = targetZ;
::urqmd::sys_.At = targetA;
::urqmd::inputs_.trspflg = 0; // nucleus as target
@@ -198,7 +194,7 @@ namespace corsika::urqmd {
::urqmd::cascinit_(::urqmd::sys_.Zt, ::urqmd::sys_.At, id);
} else {
::urqmd::inputs_.trspflg = 1; // special particle as target
- auto const [ityp, iso3] = ConvertToUrQMD(targetCode);
+ auto const [ityp, iso3] = convertToUrQMD(targetCode);
::urqmd::inputs_.spityp[1] = ityp;
::urqmd::inputs_.spiso3[1] = iso3;
}
@@ -207,49 +203,47 @@ namespace corsika::urqmd {
::urqmd::urqmd_(iflb);
// now retrieve secondaries from UrQMD
- auto const& originalCS = projectileMomentumLab.GetCoordinateSystem();
- corsika::CoordinateSystem const zAxisFrame =
- originalCS.RotateToZ(projectileMomentumLab);
+ auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
+ CoordinateSystemPtr const& zAxisFrame =
+ make_rotationToZ(originalCS, projectileMomentumLab);
for (int i = 0; i < ::urqmd::sys_.npart; ++i) {
- auto code = ConvertFromUrQMD(::urqmd::isys_.ityp[i], ::urqmd::isys_.iso3[i]);
- if (code == corsika::Code::K0 || code == corsika::Code::K0Bar) {
- code = fBooleanDist(fRNG) ? corsika::Code::K0Short : corsika::Code::K0Long;
+ auto code = convertFromUrQMD(::urqmd::isys_.ityp[i], ::urqmd::isys_.iso3[i]);
+ if (code == Code::K0 || code == Code::K0Bar) {
+ code = booleanDist_(RNG_) ? Code::K0Short : Code::K0Long;
}
// "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well
- auto momentum = corsika::Vector(
- zAxisFrame,
- corsika::QuantityVector<dimensionless_d>{
- ::urqmd::coor_.px[i], ::urqmd::coor_.py[i], ::urqmd::coor_.pz[i]} *
- 1_GeV);
+ auto momentum =
+ Vector(zAxisFrame,
+ QuantityVector<dimensionless_d>{
+ ::urqmd::coor_.px[i], ::urqmd::coor_.py[i], ::urqmd::coor_.pz[i]} *
+ 1_GeV);
- auto const energy = sqrt(momentum.squaredNorm() + square(corsika::get_mass(code)));
+ auto const energy = sqrt(momentum.getSquaredNorm() + square(get_mass(code)));
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(
- std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
- corsika::Point, TimeType>{code, energy, momentum, projectilePosition,
- projectileTime});
+ projectile.addSecondary(
+ std::make_tuple(code, energy, momentum, projectilePosition, projectileTime));
}
std::cout << "UrQMD generated " << ::urqmd::sys_.npart << " secondaries!"
<< std::endl;
}
- corsika::Code ConvertFromUrQMD(int vItyp, int vIso3) {
+ Code convertFromUrQMD(int vItyp, int vIso3) {
int const pdgInt =
::urqmd::pdgid_(vItyp, vIso3); // use the conversion function provided by UrQMD
if (pdgInt == 0) { // ::urqmd::pdgid_ returns 0 on error
throw std::runtime_error("UrQMD pdgid() returned 0");
}
- auto const pdg = static_cast<corsika::PDGCode>(pdgInt);
- return corsika::convert_from_PDG(pdg);
+ auto const pdg = static_cast<PDGCode>(pdgInt);
+ return convert_from_PDG(pdg);
}
- std::pair<int, int> ConvertToUrQMD(corsika::Code code) {
+ std::pair<int, int> convertToUrQMD(Code code) {
static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
// data mostly from github.com/afedynitch/ParticleDataTool
{22, {100, 0}}, // gamma
diff --git a/corsika/modules/pythia8/Decay.hpp b/corsika/modules/pythia8/Decay.hpp
index d80f0779b..4c3568fc7 100644
--- a/corsika/modules/pythia8/Decay.hpp
+++ b/corsika/modules/pythia8/Decay.hpp
@@ -17,29 +17,52 @@
namespace corsika::pythia8 {
- typedef corsika::Vector<hepmomentum_d> MomentumVector;
-
- class Decay : public corsika::DecayProcess<Decay> {
- const std::vector<corsika::Code> fTrackedParticles;
- int fCount = 0;
+ class Decay : public DecayProcess<Decay> {
public:
- Decay(std::vector<corsika::Code>);
+ Decay(bool const print_listing = false);
+ Decay(std::set<Code> const&);
~Decay();
- void Init();
- void SetParticleListStable(const std::vector<corsika::Code>);
- void SetUnstable(const corsika::Code);
- void SetStable(const corsika::Code);
+ // is Pythia::Decay set to handle the decay of this particle?
+ bool isDecayHandled(Code const);
+
+ //! is decay possible in principle?
+ bool canHandleDecay(Code const);
+
+ //! set Pythia::Decay to handle the decay of this particle!
+ void setHandleDecay(Code const);
+ //! set Pythia::Decay to handle the decay of this list of particles!
+ void setHandleDecay(std::vector<Code> const&);
+ //! set Pythia::Decay to handle all particle decays
+ void setHandleAllDecays();
+
+ //! print internal configuration for this particle
+ void printDecayConfig(Code const);
+ //! print configuration of decays in corsika
+ void printDecayConfig();
+
+ bool canDecay(Code const);
template <typename TParticle>
- TimeType GetLifetime(TParticle const&);
+ TimeType getLifetime(TParticle const&);
- template <typename TProjectile>
- void DoDecay(TProjectile&);
+ template <typename TView>
+ void doDecay(TView&);
private:
- Pythia8::Pythia fPythia;
+ bool isStable(Code const vCode);
+ void setStable(std::vector<Code> const&);
+ void setUnstable(Code const);
+ void setStable(Code const);
+
+ // data members
+ Pythia8::Pythia pythia_;
+
+ std::set<Code> handledDecays_;
+ int count_ = 0;
+ bool handleAllDecays_ = true;
+ bool print_listing_ = false;
};
} // namespace corsika::pythia8
diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp
index 71b58de38..eb06aefb4 100644
--- a/corsika/modules/pythia8/Interaction.hpp
+++ b/corsika/modules/pythia8/Interaction.hpp
@@ -20,31 +20,25 @@ namespace corsika::pythia8 {
class Interaction : public corsika::InteractionProcess<Interaction> {
- int fCount = 0;
- bool fInitialized = false;
-
public:
- Interaction() {}
+ Interaction(const bool print_listing = false);
~Interaction();
- void Init();
+ void setStable(std::vector<Code> const&);
+ void setUnstable(const Code);
+ void setStable(const Code);
- void SetParticleListStable(std::vector<corsika::Code> const&);
- void SetUnstable(const corsika::Code);
- void SetStable(const corsika::Code);
+ bool wasInitialized() { return initialized_; }
+ bool isValidCoMEnergy(HEPEnergyType ecm) { return (10_GeV < ecm) && (ecm < 1_PeV); }
- bool WasInitialized() { return fInitialized; }
- bool ValidCoMEnergy(HEPEnergyType ecm) { return (10_GeV < ecm) && (ecm < 1_PeV); }
+ bool canInteract(const Code);
+ void configureLabFrameCollision(const Code, const Code, const HEPEnergyType);
- bool CanInteract(const corsika::Code);
- void ConfigureLabFrameCollision(const corsika::Code, const corsika::Code,
- const HEPEnergyType);
- std::tuple<CrossSectionType, CrossSectionType> GetCrossSection(
- const corsika::Code BeamId, const corsika::Code TargetId,
- const HEPEnergyType CoMenergy);
+ std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
+ const Code BeamId, const Code TargetId, const HEPEnergyType CoMenergy);
template <typename TParticle>
- GrammageType GetInteractionLength(TParticle&);
+ GrammageType getInteractionLength(TParticle&);
/**
In this function PYTHIA is called to produce one event. The
@@ -55,11 +49,13 @@ namespace corsika::pythia8 {
void doInteraction(TProjectile&);
private:
- corsika::default_prng_type& fRNG =
- corsika::RNGManager::getInstance().getRandomStream("pythia");
- Pythia8::Pythia fPythia;
- Pythia8::SigmaTotal fSigma;
- const bool fInternalDecays = true;
+ default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("pythia");
+ Pythia8::Pythia pythia_;
+ Pythia8::SigmaTotal sigma_;
+ const bool internalDecays_ = true;
+ int count_ = 0;
+ bool initialized_ = false;
+ bool print_listing_ = false;
};
} // namespace corsika::pythia8
diff --git a/corsika/modules/pythia8/Random.hpp b/corsika/modules/pythia8/Random.hpp
index e23c4be73..0e04cdb6b 100644
--- a/corsika/modules/pythia8/Random.hpp
+++ b/corsika/modules/pythia8/Random.hpp
@@ -18,9 +18,8 @@ namespace corsika::pythia8 {
double flat();
private:
- std::uniform_real_distribution<double> fDist;
- corsika::default_prng_type& fRNG =
- corsika::RNGManager::getInstance().getRandomStream("pythia");
+ std::uniform_real_distribution<double> Dist_;
+ default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("pythia");
};
} // namespace corsika::pythia8
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
index 35bbe08f1..9fc261376 100644
--- a/corsika/modules/urqmd/UrQMD.hpp
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -22,26 +22,28 @@
namespace corsika::urqmd {
- class UrQMD : public corsika::InteractionProcess<UrQMD> {
+ class UrQMD : public InteractionProcess<UrQMD> {
public:
UrQMD();
- void Init() {}
- GrammageType GetInteractionLength(corsika::setup::Stack::StackIterator&) const;
template <typename TParticle>
- CrossSectionType GetCrossSection(TParticle const&, corsika::Code) const;
+ GrammageType getInteractionLength(TParticle const&) const;
- void doInteraction(corsika::setup::StackView::StackIterator&);
+ template <typename TParticle>
+ CrossSectionType getCrossSection(TParticle const&, Code) const;
+
+ template <typename TView>
+ void doInteraction(TView&);
- bool CanInteract(corsika::Code) const;
+ bool canInteract(Code) const;
private:
- static CrossSectionType GetCrossSection(corsika::Code, corsika::Code, HEPEnergyType,
- int);
- corsika::default_prng_type& fRNG =
- corsika::RNGManager::getInstance().getRandomStream("urqmd");
+ static CrossSectionType getCrossSection(Code, Code, HEPEnergyType, int);
+
+ // data members
+ default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("urqmd");
- std::uniform_int_distribution<int> fBooleanDist{0, 1};
+ std::uniform_int_distribution<int> booleanDist_{0, 1};
};
/**
@@ -49,8 +51,8 @@ namespace corsika::urqmd {
*
* In the current implementation a detour via the PDG code is made.
*/
- std::pair<int, int> ConvertToUrQMD(corsika::Code);
- corsika::Code ConvertFromUrQMD(int vItyp, int vIso3);
+ std::pair<int, int> convertToUrQMD(Code);
+ Code convertFromUrQMD(int vItyp, int vIso3);
} // namespace corsika::urqmd
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 335e68c7f..c96ac3a43 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -7,10 +7,10 @@ set (test_modules_sources
testObservationPlane.cpp
testQGSJetII.cpp
# testSwitchProcess.cpp -> gone
- # testPythia8.cpp
+ #testPythia8.cpp
#testUrQMD.cpp
#testCONEXSourceCut.cpp
- #testOnShellCheck.cpp
+ testOnShellCheck.cpp
testParticleCut.cpp
testSibyll.cpp
# testNullModel.cpp -> gone
diff --git a/tests/modules/testOnShellCheck.cpp b/tests/modules/testOnShellCheck.cpp
index 77c64e178..ab713f7b5 100644
--- a/tests/modules/testOnShellCheck.cpp
+++ b/tests/modules/testOnShellCheck.cpp
@@ -6,39 +6,35 @@
* the license.
*/
-#include <corsika/process/on_shell_check/OnShellCheck.h>
+#include <corsika/modules/OnShellCheck.hpp>
-#include <corsika/environment/Environment.h>
-#include <corsika/geometry/FourVector.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/media/Environment.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupStack.hpp>
#include <catch2/catch.hpp>
using namespace corsika;
-using namespace corsika::process::on_shell_check;
-using namespace corsika::units;
-using namespace corsika::units::si;
TEST_CASE("OnShellCheck", "[processes]") {
feenableexcept(FE_INVALID);
using EnvType = setup::Environment;
EnvType env;
- const geometry::CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
// setup empty particle stack
setup::Stack stack;
- stack.Clear();
+ stack.clear();
// two energies
const HEPEnergyType E = 10_GeV;
// list of arbitrary particles
- std::array const particleList{particles::Code::PiPlus, particles::Code::PiMinus,
- particles::Code::Helium, particles::Code::Gamma};
+ std::array const particleList{Code::PiPlus, Code::PiMinus, Code::Helium, Code::Gamma};
std::array const mass_shifts{1.1, 1.001, 1.0, 1.0};
@@ -46,43 +42,36 @@ TEST_CASE("OnShellCheck", "[processes]") {
OnShellCheck check(1.e-2, 0.01, false);
- check.Init();
-
// add primary particle to stack
- auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Proton, E,
- corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ auto particle = stack.addParticle(
+ std::make_tuple(Code::Proton, E, MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
// view on secondary particles
setup::StackView view{particle};
// ref. to primary particle through the secondary view.
// only this way the secondary view is populated
- auto projectile = view.GetProjectile();
+ auto projectile = view.getProjectile();
// add secondaries, all with energies above the threshold
// only cut is by species
int count = -1;
for (auto proType : particleList) {
count++;
- const auto pz = sqrt((E - particles::GetMass(proType) * mass_shifts[count]) *
- (E + particles::GetMass(proType) * mass_shifts[count]));
- auto const momentum = corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, pz});
- projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{
- proType, E, momentum, geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ const auto pz = sqrt((E - get_mass(proType) * mass_shifts[count]) *
+ (E + get_mass(proType) * mass_shifts[count]));
+ auto const momentum = MomentumVector(rootCS, {0_GeV, 0_GeV, pz});
+ projectile.addSecondary(
+ std::make_tuple(proType, E, momentum, Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
}
- check.DoSecondaries(view);
+ check.doSecondaries(view);
int i = -1;
for (auto& p : view) {
i++;
- auto const Plab = corsika::geometry::FourVector(p.GetEnergy(), p.GetMomentum());
- auto const m_kinetic = Plab.GetNorm();
+ auto const Plab = FourVector(p.getEnergy(), p.getMomentum());
+ auto const m_kinetic = Plab.getNorm();
if (i == 0)
- CHECK(m_kinetic / particles::PiPlus::GetMass() == Approx(1));
+ CHECK(m_kinetic / PiPlus::mass == Approx(1));
else if (i == 1)
- CHECK_FALSE(m_kinetic / particles::PiMinus::GetMass() == Approx(1));
+ CHECK_FALSE(m_kinetic / PiMinus::mass == Approx(1));
}
}
}
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
index e565edfdf..a3a0cffe1 100644
--- a/tests/modules/testPythia8.cpp
+++ b/tests/modules/testPythia8.cpp
@@ -13,8 +13,13 @@
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/random/RNGManager.hpp>
+#include <SetupTestEnvironment.hpp>
+#include <SetupTestStack.hpp>
+
#include <catch2/catch.hpp>
+using namespace corsika;
+
TEST_CASE("Pythia", "[processes]") {
SECTION("linking pythia") {
@@ -55,17 +60,11 @@ TEST_CASE("Pythia", "[processes]") {
}
SECTION("pythia interface") {
- using namespace corsika;
-
- const std::vector<corsika::Code> particleList = {
- corsika::Code::PiPlus, corsika::Code::PiMinus, corsika::Code::KPlus,
- corsika::Code::KMinus, corsika::Code::K0Long, corsika::Code::K0Short};
-
- corsika::RNGManager::getInstance().registerRandomStream("pythia");
+ std::set<Code> const particleList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
+ Code::KMinus, Code::K0Long, Code::K0Short};
+ RNGManager::getInstance().registerRandomStream("pythia");
corsika::pythia8::Decay model(particleList);
-
- model.Init();
}
}
@@ -84,89 +83,86 @@ TEST_CASE("Pythia", "[processes]") {
#include <corsika/media/NuclearComposition.hpp>
using namespace corsika;
-using namespace corsika::units::si;
template <typename TStackView>
-auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
- geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
-
- for (auto const& p : view) { sum += p.GetMomentum(); }
-
+auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
+ MomentumVector sum{vCS, 0_eV, 0_eV, 0_eV};
+ for (auto const& p : view) { sum += p.getMomentum(); }
return sum;
}
TEST_CASE("pythia process") {
- auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Proton);
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
[[maybe_unused]] auto const& node_dummy = nodePtr;
SECTION("pythia decay") {
- feenableexcept(FE_INVALID);
- auto [stackPtr, secViewPtr] =
- setup::testing::setupStack(particles::Code::PiPlus, 0, 0, P0, nodePtr, *csPtr);
-
const HEPEnergyType E0 = 10_GeV;
- HEPMomentumType P0 = sqrt(E0 * E0 - corsika::PiPlus::mass * corsika::PiPlus::mass);
- auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
- corsika::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle = stack.AddParticle(
- std::tuple<corsika::Code, units::si::HEPEnergyType, corsika::MomentumVector,
- corsika::Point, units::si::TimeType>{corsika::Code::PiPlus, E0, plab,
- pos, 0_ns});
+ HEPMomentumType P0 = sqrt(E0 * E0 - PiPlus::mass * PiPlus::mass);
+
+ // feenableexcept(FE_INVALID); \todo how does this work nowadays...???
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::PiPlus, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+ auto& stack = *stackPtr;
+ auto& view = *secViewPtr;
+
+ auto plab = MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ Point pos(cs, 0_m, 0_m, 0_m);
+ auto particle =
+ stackPtr->addParticle(std::make_tuple(Code::PiPlus, E0, plab, pos, 0_ns));
- const std::vector<corsika::Code> particleList = {
- corsika::Code::PiPlus, corsika::Code::PiMinus, corsika::Code::KPlus,
- corsika::Code::KMinus, corsika::Code::K0Long, corsika::Code::K0Short};
+ std::set<Code> const particleList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
+ Code::KMinus, Code::K0Long, Code::K0Short};
- corsika::RNGManager::getInstance().registerRandomStream("pythia");
+ RNGManager::getInstance().registerRandomStream("pythia");
corsika::pythia8::Decay model(particleList);
- [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
- model.DoDecay(view);
+ [[maybe_unused]] const TimeType time = model.getLifetime(particle);
+ model.doDecay(*secViewPtr);
CHECK(stack.getEntries() == 3);
auto const pSum = sumMomentum(view, cs);
- CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(1e-4));
- CHECK((pSum.norm() - plab.norm()) / 1_GeV == Approx(0).margin(1e-4));
+ CHECK((pSum - plab).getNorm() / 1_GeV == Approx(0).margin(1e-4));
+ CHECK((pSum.getNorm() - plab.getNorm()) / 1_GeV == Approx(0).margin(1e-4));
}
SECTION("pythia decay config") {
- process::pythia::Decay model({particles::Code::PiPlus, particles::Code::PiMinus});
- REQUIRE(model.IsDecayHandled(particles::Code::PiPlus));
- REQUIRE(model.IsDecayHandled(particles::Code::PiMinus));
- REQUIRE_FALSE(model.IsDecayHandled(particles::Code::KPlus));
+ corsika::pythia8::Decay model({Code::PiPlus, Code::PiMinus});
+ CHECK(model.isDecayHandled(Code::PiPlus));
+ CHECK(model.isDecayHandled(Code::PiMinus));
+ CHECK_FALSE(model.isDecayHandled(Code::KPlus));
- const std::vector<particles::Code> particleTestList = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::Lambda0Bar, particles::Code::D0Bar};
+ const std::vector<Code> particleTestList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
+ Code::Lambda0Bar, Code::D0Bar};
// setup decays
- model.SetHandleDecay(particleTestList);
- for (auto& pCode : particleTestList) REQUIRE(model.IsDecayHandled(pCode));
+ model.setHandleDecay(particleTestList);
+ for (auto& pCode : particleTestList) CHECK(model.isDecayHandled(pCode));
// individually
- model.SetHandleDecay(particles::Code::KMinus);
+ model.setHandleDecay(Code::KMinus);
// possible decays
- REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Proton));
- REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Electron));
- REQUIRE(model.CanHandleDecay(particles::Code::PiPlus));
- REQUIRE(model.CanHandleDecay(particles::Code::MuPlus));
+ CHECK_FALSE(model.canHandleDecay(Code::Proton));
+ CHECK_FALSE(model.canHandleDecay(Code::Electron));
+ CHECK(model.canHandleDecay(Code::PiPlus));
+ CHECK(model.canHandleDecay(Code::MuPlus));
}
SECTION("pythia interaction") {
- feenableexcept(FE_INVALID);
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::PiPlus, 0,
- 0, 100_GeV, nodePtr, *csPtr);
+ // feenableexcept(FE_INVALID); \todo how does this work nowadays
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::PiPlus, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+ *csPtr);
auto& view = *secViewPtr;
auto particle = stackPtr->first();
- process::pythia::Interaction model;
-
- [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
- [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+ corsika::pythia8::Interaction model;
+ model.doInteraction(view);
+ [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
+ CHECK(length == 50_g / square(1_cm));
}
}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index a51081fe4..13a3a572e 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -17,8 +17,8 @@
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
+#include <SetupTestStack.hpp>
+#include <SetupTestEnvironment.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -35,136 +35,126 @@ using namespace corsika::urqmd;
template <typename TStackView>
auto sumCharge(TStackView const& view) {
int totalCharge = 0;
-
- for (auto const& p : view) { totalCharge += corsika::charge_number(p.GetPID()); }
-
+ for (auto const& p : view) { totalCharge += get_charge_number(p.getPID()); }
return totalCharge;
}
template <typename TStackView>
-auto sumMomentum(TStackView const& view, corsika::CoordinateSystem const& vCS) {
- corsika::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
-
- for (auto const& p : view) { sum += p.GetMomentum(); }
-
+auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
+ MomentumVector sum{vCS, 0_eV, 0_eV, 0_eV};
+ for (auto const& p : view) { sum += p.getMomentum(); }
return sum;
}
TEST_CASE("UrQMD") {
SECTION("conversion") {
- REQUIRE_THROWS(corsika::urqmd::ConvertFromUrQMD(106, 0));
- REQUIRE(corsika::urqmd::ConvertFromUrQMD(101, 0) == corsika::Code::Pi0);
- REQUIRE(corsika::urqmd::ConvertToUrQMD(corsika::Code::PiPlus) ==
- std::make_pair<int, int>(101, 2));
+ CHECK_THROWS(corsika::urqmd::convertFromUrQMD(106, 0));
+ CHECK(corsika::urqmd::convertFromUrQMD(101, 0) == Code::Pi0);
+ CHECK(corsika::urqmd::convertToUrQMD(Code::PiPlus) ==
+ std::make_pair<int, int>(101, 2));
}
feenableexcept(FE_INVALID);
- corsika::RNGManager::getInstance().registerRandomStream("urqmd");
+ RNGManager::getInstance().registerRandomStream("urqmd");
UrQMD urqmd;
SECTION("interaction length") {
- auto [env, csPtr, nodePtr] =
- setup::testing::setupEnvironment(particles::Code::Nitrogen);
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Nitrogen);
auto const& cs = *csPtr;
{ [[maybe_unused]] auto const& env_dummy = env; }
- corsika::Code validProjectileCodes[] = {
- corsika::Code::PiPlus, corsika::Code::PiMinus, corsika::Code::Proton,
- corsika::Code::Neutron, corsika::Code::KPlus, corsika::Code::KMinus,
- corsika::Code::K0, corsika::Code::K0Bar, corsika::Code::K0Long};
+ Code validProjectileCodes[] = {Code::PiPlus, Code::PiMinus, Code::Proton,
+ Code::Neutron, Code::KPlus, Code::KMinus,
+ Code::K0, Code::K0Bar, Code::K0Long};
for (auto code : validProjectileCodes) {
- auto [stack, view] = setup::testing::setupStack(code, 0, 0, 100_GeV, nodePtr, cs);
- REQUIRE(stack->getEntries() == 1);
- REQUIRE(view->getEntries() == 0);
+ 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);
// simple check whether the cross-section is non-vanishing
- REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), corsika::Code::Proton) / 1_mb >
- 0);
- REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), corsika::Code::Nitrogen) /
- 1_mb >
- 0);
- REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), corsika::Code::Oxygen) / 1_mb >
- 0);
- REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), corsika::Code::Argon) / 1_mb >
- 0);
+ 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);
}
}
SECTION("nucleus projectile") {
- auto [env, csPtr, nodePtr] =
- setup::testing::setupEnvironment(particles::Code::Oxygen);
+ 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
unsigned short constexpr A = 14, Z = 7;
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
- Z, 400_GeV, nodePtr, *csPtr);
- REQUIRE(stackPtr->getEntries() == 1);
- REQUIRE(secViewPtr->getEntries() == 0);
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Nucleus, A, Z, 400_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();
- [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(*secViewPtr);
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
+ urqmd.doInteraction(*secViewPtr);
- REQUIRE(sumCharge(*secViewPtr) == Z + corsika::charge_number(corsika::Code::Oxygen));
+ CHECK(sumCharge(*secViewPtr) == Z + get_charge_number(Code::Oxygen));
auto const secMomSum =
- sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
- REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
- Approx(0).margin(1e-2));
+ sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
+ CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
+ Approx(0).margin(1e-2));
}
SECTION("\"special\" projectile") {
- auto [env, csPtr, nodePtr] =
- setup::testing::setupEnvironment(particles::Code::Oxygen);
+ 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
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::PiPlus, 0,
- 0, 400_GeV, nodePtr, *csPtr);
- REQUIRE(stackPtr->getEntries() == 1);
- REQUIRE(secViewPtr->getEntries() == 0);
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::PiPlus, 0, 0, 400_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();
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
- [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(*secViewPtr);
+ urqmd.doInteraction(*secViewPtr);
- REQUIRE(sumCharge(*secViewPtr) == corsika::charge_number(corsika::Code::PiPlus) +
- corsika::charge_number(corsika::Code::Oxygen));
+ CHECK(sumCharge(*secViewPtr) ==
+ get_charge_number(Code::PiPlus) + get_charge_number(Code::Oxygen));
auto const secMomSum =
- sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
- REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
- Approx(0).margin(1e-2));
+ sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
+ CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
+ Approx(0).margin(1e-2));
}
SECTION("K0Long projectile") {
- auto [env, csPtr, nodePtr] =
- setup::testing::setupEnvironment(particles::Code::Oxygen);
+ 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
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::K0Long, 0,
- 0, 400_GeV, nodePtr, *csPtr);
- REQUIRE(stackPtr->getEntries() == 1);
- REQUIRE(secViewPtr->getEntries() == 0);
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::K0Long, 0, 0, 400_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();
+ auto projectile = secViewPtr->getProjectile();
+ auto const projectileMomentum = projectile.getMomentum();
- [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(*secViewPtr);
+ urqmd.doInteraction(*secViewPtr);
- REQUIRE(sumCharge(*secViewPtr) == corsika::charge_number(corsika::Code::K0Long) +
- corsika::charge_number(corsika::Code::Oxygen));
+ CHECK(sumCharge(*secViewPtr) ==
+ get_charge_number(Code::K0Long) + get_charge_number(Code::Oxygen));
auto const secMomSum =
- sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
- REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
- Approx(0).margin(1e-2));
+ sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
+ CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
+ Approx(0).margin(1e-2));
}
}
--
GitLab