IAP GITLAB

Commit 3c70ce61 authored by Ralf Ulrich's avatar Ralf Ulrich

onshellcheck

parent 945e2aee
......@@ -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/
......
......@@ -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; }
}
......@@ -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;
}
......
This diff is collapsed.
......@@ -12,6 +12,6 @@
namespace corsika::pythia8 {
double Random::flat() { return fDist(fRNG); }
double Random::flat() { return Dist_(RNG_); }
} // namespace corsika::pythia8
This diff is collapsed.
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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
......
......@@ -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));
}
}
}
......@@ -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));
}
}
......@@ -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,