diff --git a/Processes/InteractionCounter/CMakeLists.txt b/Processes/InteractionCounter/CMakeLists.txt index 4ce877826cc206270220754cad4bcb206bd3fe97..30b8815312afe511b787c84a5dbee1efff274cf6 100644 --- a/Processes/InteractionCounter/CMakeLists.txt +++ b/Processes/InteractionCounter/CMakeLists.txt @@ -33,11 +33,15 @@ install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE}) # -------------------- # code unit testing -#~ CORSIKA_ADD_TEST(testLuminosityCounter) -#~ target_link_libraries ( - #~ testLuminosityCounter - #~ ProcessLuminosityCounter - #~ CORSIKAstackinterface - #~ CORSIKAtesting -#~ ) +CORSIKA_ADD_TEST(testInteractionCounter SOURCES + testInteractionCounter.cc + ${MODEL_HEADERS} +) +target_link_libraries ( + testInteractionCounter + ProcessInteractionCounter + CORSIKAunits + CORSIKAenvironment + CORSIKAtesting + ) diff --git a/Processes/InteractionCounter/InteractionCounter.h b/Processes/InteractionCounter/InteractionCounter.h index c65fd1f61d6d55dee535428161ef453f1d80eecf..a6fb6b43b8f8bb245697b935d0d54810ad2a1506 100644 --- a/Processes/InteractionCounter/InteractionCounter.h +++ b/Processes/InteractionCounter/InteractionCounter.h @@ -75,8 +75,6 @@ namespace corsika::process::interaction_counter { using namespace units::si; if (projectile_id == particles::Code::Nucleus) { - std::cerr << "NUCLEUS " << lab_energy << " " << A << " " << Z << std::endl; - auto const sqrtS = sqrt(A * A * (units::constants::nucleonMass * units::constants::nucleonMass) + mass_target * mass_target + 2 * lab_energy * mass_target); diff --git a/Processes/InteractionCounter/testInteractionCounter.cc b/Processes/InteractionCounter/testInteractionCounter.cc new file mode 100644 index 0000000000000000000000000000000000000000..797d0f87d1205fb0db28ae66b1ca21c50e709d3b --- /dev/null +++ b/Processes/InteractionCounter/testInteractionCounter.cc @@ -0,0 +1,156 @@ + +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * See file AUTHORS for a list of contributors. + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <corsika/process/interaction_counter/InteractionCounter.h> + +#include <corsika/environment/Environment.h> +#include <corsika/environment/HomogeneousMedium.h> +#include <corsika/environment/NuclearComposition.h> +#include <corsika/geometry/Point.h> +#include <corsika/geometry/RootCoordinateSystem.h> +#include <corsika/geometry/Vector.h> +#include <corsika/units/PhysicalUnits.h> + +#include <corsika/setup/SetupStack.h> + +#include <catch2/catch.hpp> + +#include <numeric> + +using namespace corsika; +using namespace corsika::process::interaction_counter; +using namespace corsika::units; +using namespace corsika::units::si; + +auto setupEnvironment(particles::Code target_code) { + // setup environment, geometry + auto env = std::make_unique<environment::Environment<environment::IMediumModel>>(); + auto& universe = *(env->GetUniverse()); + const geometry::CoordinateSystem& cs = env->GetCoordinateSystem(); + + auto theMedium = + environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>( + geometry::Point{cs, 0_m, 0_m, 0_m}, + 1_km * std::numeric_limits<double>::infinity()); + + using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>; + theMedium->SetModelProperties<MyHomogeneousModel>( + 1_kg / (1_m * 1_m * 1_m), + environment::NuclearComposition(std::vector<particles::Code>{target_code}, + std::vector<float>{1.})); + + auto const* nodePtr = theMedium.get(); + universe.AddChild(std::move(theMedium)); + + return std::make_tuple(std::move(env), &cs, nodePtr); +} + +template <typename TNodeType> +auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr, + geometry::CoordinateSystem const& cs) { + auto stack = std::make_unique<setup::Stack>(); + auto constexpr mN = corsika::units::constants::nucleonMass; + + 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 = + stack->AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, + corsika::stack::MomentumVector, geometry::Point, + units::si::TimeType, unsigned short, unsigned short>{ + particles::Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ}); + + particle.SetNode(vNodePtr); + return std::make_tuple( + std::move(stack), + std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle)); +} + +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)); +} + +struct DummyProcess { + template <typename TParticle> + GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) { + return 100_g / 1_cm / 1_cm; + } + + template <typename TParticle> + auto DoInteraction([[maybe_unused]] TParticle& projectile) { + return nullptr; + } +}; + +TEST_CASE("InteractionCounter") { + DummyProcess d; + InteractionCounter countedProcess(d); + + SECTION("GetInteractionLength") { + REQUIRE(countedProcess.GetInteractionLength(nullptr) == 100_g / 1_cm / 1_cm); + } + + auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen); + SECTION("DoInteraction nucleus") { + unsigned short constexpr A = 14, Z = 7; + auto [stackPtr, secViewPtr] = setupStack(A, Z, 105_TeV, nodePtr, *csPtr); + + auto projectile = secViewPtr->GetProjectile(); + auto const ret = countedProcess.DoInteraction(projectile); + REQUIRE(ret == nullptr); + + auto const& h = countedProcess.GetHistogram().labHists().second.at(1'000'070'140); + REQUIRE(h.at(50) == 1); + REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1); + + auto const& h2 = countedProcess.GetHistogram().CMSHists().second.at(1'000'070'140); + REQUIRE(h2.at(32) == 1); // bin 1.584 .. 1.995 TeV √s + REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1); + } + + SECTION("DoInteraction Lambda") { + auto constexpr code = particles::Code::Lambda0; + auto constexpr codeInt = static_cast<particles::CodeIntType>(code); + auto [stackPtr, secViewPtr] = setupStack(code, 105_TeV, nodePtr, *csPtr); + + auto projectile = secViewPtr->GetProjectile(); + auto const ret = countedProcess.DoInteraction(projectile); + REQUIRE(ret == nullptr); + + auto const& h = countedProcess.GetHistogram().labHists().first; + REQUIRE(h.at(codeInt, 50) == 1); + REQUIRE(std::accumulate(h.cbegin(), h.cend(), 0) == 1); + + auto const& h2 = countedProcess.GetHistogram().CMSHists().first; + REQUIRE(h2.at(codeInt, 32) == 1); // bin 1.584 .. 1.995 TeV √s + REQUIRE(std::accumulate(h2.cbegin(), h2.cend(), 0) == 1); + } +}