IAP GITLAB

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

restructured test

parent 9e9d3cfc
No related branches found
No related tags found
1 merge request!116Some improvements here and there
...@@ -25,22 +25,22 @@ template <typename TParticle> // need template here, as this is called both with ...@@ -25,22 +25,22 @@ template <typename TParticle> // need template here, as this is called both with
// SetupParticle as well as SetupProjectile // SetupParticle as well as SetupProjectile
CrossSectionType UrQMD::GetCrossSection( CrossSectionType UrQMD::GetCrossSection(
TParticle const& vProjectile, TParticle const& vProjectile,
corsika::particles::Code vTargetCode, corsika::particles::Code vTargetCode)
corsika::units::si::HEPEnergyType vProjectileEnergyLab)
const { const {
using namespace units::si; using namespace units::si;
// TODO: energy cuts, return 0 for non-hadrons // TODO: energy cuts, return 0 for non-hadrons
auto const projectileCode = vProjectile.GetPID(); auto const projectileCode = vProjectile.GetPID();
auto const projectileEnergyLab = vProjectile.GetEnergy();
// the following is a translation of ptsigtot() into C++ // the following is a translation of ptsigtot() into C++
if (projectileCode != particles::Code::Nucleus && if (projectileCode != particles::Code::Nucleus &&
!IsNucleus(vTargetCode)) { // both particles are "special" !IsNucleus(vTargetCode)) { // both particles are "special"
auto const mProj = particles::GetMass(projectileCode); auto const mProj = particles::GetMass(projectileCode);
auto const mTar = particles::GetMass(vTargetCode); auto const mTar = particles::GetMass(vTargetCode);
double sqrtS = sqrt(detail::static_pow<2>(mProj) + detail::static_pow<2>(mTar) + double sqrtS = sqrt(units::si::detail::static_pow<2>(mProj) + units::si::detail::static_pow<2>(mTar) +
2 * vProjectileEnergyLab * mTar) * 2 * projectileEnergyLab * mTar) *
(1 / 1_GeV); (1 / 1_GeV);
// we must set some UrQMD globals first... // we must set some UrQMD globals first...
...@@ -61,7 +61,7 @@ template <typename TParticle> // need template here, as this is called both with ...@@ -61,7 +61,7 @@ template <typename TParticle> // need template here, as this is called both with
int const At = IsNucleus(vTargetCode) ? particles::GetNucleusA(vTargetCode) : 1; int const At = IsNucleus(vTargetCode) ? particles::GetNucleusA(vTargetCode) : 1;
double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1]; double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1];
return 10_mb * M_PI * detail::static_pow<2>(maxImpact); return 10_mb * M_PI * units::si::detail::static_pow<2>(maxImpact);
// is a constant cross-section really reasonable? // is a constant cross-section really reasonable?
} }
} }
...@@ -72,8 +72,7 @@ GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&) ...@@ -72,8 +72,7 @@ GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&)
using namespace std::placeholders; using namespace std::placeholders;
CrossSectionType const weightedProdCrossSection = mediumComposition.WeightedSum( CrossSectionType const weightedProdCrossSection = mediumComposition.WeightedSum(
std::bind(&UrQMD::GetCrossSection<decltype(vParticle)>, this, vParticle, _1, std::bind(&UrQMD::GetCrossSection<decltype(vParticle)>, this, vParticle, _1));
vParticle.GetEnergy()));
return mediumComposition.GetAverageMassNumber() * units::constants::u / return mediumComposition.GetAverageMassNumber() * units::constants::u /
weightedProdCrossSection; weightedProdCrossSection;
...@@ -97,7 +96,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti ...@@ -97,7 +96,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti
crossSections.reserve(components.size()); crossSections.reserve(components.size());
for (auto const c : components) { for (auto const c : components) {
crossSections.push_back(GetCrossSection(vProjectile, c, projectileEnergyLab)); crossSections.push_back(GetCrossSection(vProjectile, c));
} }
return crossSections; return crossSections;
......
...@@ -20,8 +20,7 @@ namespace corsika::process::UrQMD { ...@@ -20,8 +20,7 @@ namespace corsika::process::UrQMD {
template <typename TParticle> template <typename TParticle>
corsika::units::si::CrossSectionType GetCrossSection( corsika::units::si::CrossSectionType GetCrossSection(
TParticle const&, corsika::particles::Code, TParticle const&, corsika::particles::Code) const;
corsika::units::si::HEPEnergyType) const;
corsika::process::EProcessReturn DoInteraction( corsika::process::EProcessReturn DoInteraction(
corsika::setup::StackView::StackIterator&); corsika::setup::StackView::StackIterator&);
......
...@@ -28,6 +28,8 @@ ...@@ -28,6 +28,8 @@
#include <corsika/environment/HomogeneousMedium.h> #include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h> #include <corsika/environment/NuclearComposition.h>
#include <tuple>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one #define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file // cpp file
#include <catch2/catch.hpp> #include <catch2/catch.hpp>
...@@ -52,15 +54,12 @@ auto sumCharge(TStack& stack) { ...@@ -52,15 +54,12 @@ auto sumCharge(TStack& stack) {
return totalCharge; return totalCharge;
} }
TEST_CASE("UrQMD") {
feenableexcept(FE_INVALID);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
UrQMD urqmd;
// setup environment, geometry auto setupEnvironment(particles::Code vTargetCode) {
environment::Environment<environment::IMediumModel> env; // setup environment, geometry
auto& universe = *(env.GetUniverse()); auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem(); auto& universe = *(env->GetUniverse());
const geometry::CoordinateSystem& cs = env->GetCoordinateSystem();
auto theMedium = auto theMedium =
environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>( environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
...@@ -71,59 +70,90 @@ TEST_CASE("UrQMD") { ...@@ -71,59 +70,90 @@ TEST_CASE("UrQMD") {
theMedium->SetModelProperties<MyHomogeneousModel>( theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m), 1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition( environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.})); std::vector<particles::Code>{vTargetCode}, std::vector<float>{1.}));
auto const* nodePtr = theMedium.get(); auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium)); universe.AddChild(std::move(theMedium));
return std::make_tuple(std::move(env), &cs, nodePtr);
}
geometry::Point const origin(cs, {0_m, 0_m, 0_m}); template <typename TNodeType>
geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second, auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr, geometry::CoordinateSystem const& cs) {
1_m / second); auto stack = std::make_unique<setup::Stack>();
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
const HEPEnergyType P0 = 1000_GeV;
auto pLab = corsika::stack::MomentumVector(cs, {P0, 0_GeV, 0_GeV});
SECTION("nucleon projectile") {
setup::Stack stack;
unsigned short constexpr A = 16, Z = 8;
auto constexpr mN = corsika::units::constants::nucleonMass; auto constexpr mN = corsika::units::constants::nucleonMass;
HEPMomentumType E0 = sqrt(A * A * mN * mN + P0 * P0);
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 = sqrt(units::si::detail::static_pow<2>(mN*vA) + pLab.squaredNorm());
auto particle = auto particle =
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, stack->AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{ units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, pLab, origin, 0_ns, A, Z}); particles::Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ});
particle.SetNode(nodePtr); particle.SetNode(vNodePtr);
corsika::stack::SecondaryView view(particle); return std::make_tuple(std::move(stack), std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
auto projectile = view.GetProjectile(); }
template <typename TNodeType>
auto setupStack(particles::Code vProjectileType, HEPEnergyType vMomentum, TNodeType* vNodePtr, geometry::CoordinateSystem const& cs) {
auto stack = std::make_unique<setup::Stack>();
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 = sqrt(units::si::detail::static_pow<2>(particles::GetMass(vProjectileType)) + pLab.squaredNorm());
auto particle =
stack->AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType>{
vProjectileType, E0, pLab, origin, 0_ns});
particle.SetNode(vNodePtr);
return std::make_tuple(std::move(stack), std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
}
//~ int main() {
TEST_CASE("UrQMD") {
feenableexcept(FE_INVALID);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
UrQMD urqmd;
SECTION("cross sections") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Invalid);
auto const& cs = *csPtr;
particles::Code validProjectileCodes[] = {particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::Proton, particles::Code::Neutron};
for (auto code: validProjectileCodes) {
auto [stack, view] = setupStack(code, 100_GeV, nodePtr, cs);
REQUIRE( stack->GetSize() == 1);
// simple check whether the cross-section is non-vanishing
REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Proton) / 1_mb > 0);
REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Nitrogen) / 1_mb > 0);
REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Oxygen) / 1_mb > 0);
REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Nitrogen) / 1_mb > 0);
}
}
SECTION("nucleon projectile") {
unsigned short constexpr A = 16, Z = 8;
[[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile); [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(stack) == Z + particles::GetChargeNumber(particles::Code::Oxygen)); REQUIRE(sumCharge(stack) == Z + particles::GetChargeNumber(particles::Code::Oxygen));
} }
SECTION("\"special\" projectile") { //~ SECTION("\"special\" projectile") {
setup::Stack stack;
auto constexpr code = particles::Code::PiPlus; //~ [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
auto constexpr mass = particles::GetMass(code);
HEPMomentumType E0 = sqrt(mass * mass + pLab.squaredNorm());
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
code, E0, pLab, origin, 0_ns});
particle.SetNode(nodePtr);
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
[[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile); //~ REQUIRE(sumCharge(stack) == particles::GetChargeNumber(particles::Code::PiPlus) +
//~ particles::GetChargeNumber(particles::Code::Oxygen));
REQUIRE(sumCharge(stack) == particles::GetChargeNumber(particles::Code::PiPlus) + //~ }
particles::GetChargeNumber(particles::Code::Oxygen));
}
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment