diff --git a/Data b/Data index b65a9961760359531f4ed8668e61d4e1350c7437..7c276813de7b4f7635eb024b3806936a4bc9c770 160000 --- a/Data +++ b/Data @@ -1 +1 @@ -Subproject commit b65a9961760359531f4ed8668e61d4e1350c7437 +Subproject commit 7c276813de7b4f7635eb024b3806936a4bc9c770 diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc index 8258b3c5eae4c52ab94633adae74f48a673a3d4b..ef22b9dd350cda32a41d9f350da99d0ef3ed4d36 100644 --- a/Processes/QGSJetII/testQGSJetII.cc +++ b/Processes/QGSJetII/testQGSJetII.cc @@ -22,6 +22,26 @@ using namespace corsika; using namespace corsika::process::qgsjetII; +using namespace corsika::units::si; + +template <typename TStackView> +auto sumCharge(TStackView const& view) { + int totalCharge = 0; + + for (auto const& p : view) { totalCharge += particles::GetChargeNumber(p.GetPID()); } + + return totalCharge; +} + +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(); } + + return sum; +} + TEST_CASE("QgsjetII", "[processes]") { @@ -113,18 +133,27 @@ TEST_CASE("QgsjetIIInterface", "[processes]") { auto particle = stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, geometry::Point, - units::si::TimeType, unsigned int, unsigned int>{ - particles::Code::Nucleus, E0, plab, pos, 0_ns, 16, 8}); - // corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ - // particles::Code::PiPlus, E0, plab, pos, 0_ns}); + units::si::TimeType>{ + particles::Code::Proton, E0, plab, pos, 0_ns}); particle.SetNode(nodePtr); corsika::stack::SecondaryView view(particle); auto projectile = view.GetProjectile(); + auto const projectileMomentum = projectile.GetMomentum(); Interaction model; model.Init(); [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile); [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle); + + REQUIRE( length/(1_g / square(1_cm)) == Approx(93.47).margin(0.1) ); + REQUIRE( view.GetSize() == 13 ); + /*REQUIRE( sumCharge(view) == + 1 + particles::GetChargeNumber(particles::Code::Oxygen) );*/ + auto const secMomSum = + sumMomentum(view, projectileMomentum.GetCoordinateSystem()); + REQUIRE( (secMomSum - projectileMomentum).norm() / projectileMomentum.norm() == + Approx(0).margin(1e-2)); + } } diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc index e48fb43814caaa50fbd341eaa615573f1df4bc20..ecc75e185ec36d08bc530d7bceb1481965514380 100644 --- a/Processes/UrQMD/UrQMD.cc +++ b/Processes/UrQMD/UrQMD.cc @@ -27,15 +27,13 @@ using namespace corsika::process::UrQMD; using namespace corsika::units::si; -UrQMD::UrQMD() { iniurqmdc8_(); } - using SetupStack = corsika::setup::Stack; using SetupParticle = corsika::setup::Stack::StackIterator; using SetupProjectile = corsika::setup::StackView::StackIterator; UrQMD::UrQMD(std::string const& xs_file) { readXSFile(xs_file); - iniurqmd_(); + iniurqmdc8_(); } CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode, diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h index ac24beb6f0a6198da4664702905745ff14711165..ef1f0909dc648bda8e3e418dae3ed5ef3672158a 100644 --- a/Processes/UrQMD/UrQMD.h +++ b/Processes/UrQMD/UrQMD.h @@ -77,74 +77,74 @@ namespace corsika::process::UrQMD { using nmaxDoubleArray = nmaxArray<double>; extern "C" { - // FORTRAN functions defined in UrQMD - void iniurqmdc8_(); - double ranf_(int&); - void cascinit_(int const&, int const&, int const&); - double nucrad_(int const&); - void urqmd_(int&); - int pdgid_(int const&, int const&); - double sigtot_(int&, int&, double&); - int collclass_(int&, int&, int&, int&); - double crossx_(int const&, double const&, int const&, int const&, double const&, - int const&, int const&, double const&, double&); - int readsigmaln_(int const&, int const&, int const&); - - // defined in coms.f - extern struct { - int npart, nbar, nmes, ctag, nsteps, uid_cnt, ranseed, event; - int Ap; // projectile mass number (in case of nucleus) - int At; // target mass number (in case of nucleus) - int Zp; // projectile charge number (in case of nucleus) - int Zt; // target charge number (in case of nucleus) - int eos, dectag, NHardRes, NSoftRes, NDecRes, NElColl, NBlColl; - } sys_; - - extern struct { - double time, acttime, bdist, bimp, bmin; - double ebeam; // lab-frame energy of projectile - double ecm; - } rsys_; - - // defined in coms.f - extern struct { - nmaxIntArray spin, ncoll, charge, ityp, lstcoll, iso3, origin, strid, uid; - } isys_; - - // defined in coor.f - extern struct { - nmaxDoubleArray r0, rx, ry, rz, p0, px, py, pz, fmass, rww, dectime; - } coor_; - - // defined in inputs.f - extern struct { - int nevents; - std::array<int, 2> spityp; // particle codes of: [0]: projectile, [1]: target - int prspflg; // projectile special flag - int trspflg; // target special flag, set to 1 unless target is nucleus > H - std::array<int, 2> spiso3; // particle codes of: [0]: projectile, [1]: target - int outsteps, bflag, srtflag, efuncflag, nsrt, npb, firstev; - } inputs_; - - // defined in inputs.f - extern struct { - double srtmin, srtmax, pbeam, betann, betatar, betapro, pbmin, pbmax; - } input2_; - - // defined in options.f - extern struct { - std::array<double, details::constants::numcto> CTOption; - std::array<double, details::constants::numctp> CTParam; - } options_; - - extern struct { - int fixedseed, bf13, bf14, bf15, bf16, bf17, bf18, bf19, bf20; - } loptions_; - - // defined in urqmdInterface.F - extern struct { std::array<double, 3> xs, bim; } cxs_u2_; + // FORTRAN functions defined in UrQMD + void iniurqmdc8_(); + double ranf_(int&); + void cascinit_(int const&, int const&, int const&); + double nucrad_(int const&); + void urqmd_(int&); + int pdgid_(int const&, int const&); + double sigtot_(int&, int&, double&); + int collclass_(int&, int&, int&, int&); + double crossx_(int const&, double const&, int const&, int const&, double const&, + int const&, int const&, double const&, double&); + int readsigmaln_(int const&, int const&, int const&); + + // defined in coms.f + extern struct { + int npart, nbar, nmes, ctag, nsteps, uid_cnt, ranseed, event; + int Ap; // projectile mass number (in case of nucleus) + int At; // target mass number (in case of nucleus) + int Zp; // projectile charge number (in case of nucleus) + int Zt; // target charge number (in case of nucleus) + int eos, dectag, NHardRes, NSoftRes, NDecRes, NElColl, NBlColl; + } sys_; + + extern struct { + double time, acttime, bdist, bimp, bmin; + double ebeam; // lab-frame energy of projectile + double ecm; + } rsys_; + + // defined in coms.f + extern struct { + nmaxIntArray spin, ncoll, charge, ityp, lstcoll, iso3, origin, strid, uid; + } isys_; + + // defined in coor.f + extern struct { + nmaxDoubleArray r0, rx, ry, rz, p0, px, py, pz, fmass, rww, dectime; + } coor_; + + // defined in inputs.f + extern struct { + int nevents; + std::array<int, 2> spityp; // particle codes of: [0]: projectile, [1]: target + int prspflg; // projectile special flag + int trspflg; // target special flag, set to 1 unless target is nucleus > H + std::array<int, 2> spiso3; // particle codes of: [0]: projectile, [1]: target + int outsteps, bflag, srtflag, efuncflag, nsrt, npb, firstev; + } inputs_; + + // defined in inputs.f + extern struct { + double srtmin, srtmax, pbeam, betann, betatar, betapro, pbmin, pbmax; + } input2_; + + // defined in options.f + extern struct { + std::array<double, details::constants::numcto> CTOption; + std::array<double, details::constants::numctp> CTParam; + } options_; + + extern struct { + int fixedseed, bf13, bf14, bf15, bf16, bf17, bf18, bf19, bf20; + } loptions_; + + // defined in urqmdInterface.F + extern struct { std::array<double, 3> xs, bim; } cxs_u2_; } - + /** * convert CORSIKA code to UrQMD code tuple * diff --git a/ThirdParty/CMakeLists.txt b/ThirdParty/CMakeLists.txt index b8a1f16815167c93e5c18eecf8ef99f94929e797..58d65940da36242540218ec3a64ae8037c03af66 100644 --- a/ThirdParty/CMakeLists.txt +++ b/ThirdParty/CMakeLists.txt @@ -227,10 +227,13 @@ endif (NOT (${USE_CONEX_C8} IN_LIST ThirdPartyChoiceValues)) message (STATUS "USE_CONEX_C8='${USE_CONEX_C8}'") add_library (C8::ext::conex STATIC IMPORTED GLOBAL) +# this is from the corsika-data/readLib submodule: get_directory_property (Boost_iostreams_FOUND DIRECTORY ${CORSIKA_DATA}/readLib DEFINITION Boost_iostreams_FOUND) set (conex_boost "") +set (boost_lib "boost_iostreams") if (NOT Boost_iostreams_FOUND) set (conex_boost "NO_BOOST=1") + set (boost_lib "") endif (NOT Boost_iostreams_FOUND) if ("x_${USE_CONEX_C8}" STREQUAL "x_SYSTEM") @@ -238,9 +241,12 @@ if ("x_${USE_CONEX_C8}" STREQUAL "x_SYSTEM") message (STATUS "Using system-level CONEX version at ${CONEX_INCLUDE_DIR}") set_target_properties ( C8::ext::conex PROPERTIES - IMPORTED_LOCATION ${CONEX_LIBRARY} - IMPORTED_LINK_INTERFACE_LIBRARIES dl - INTERFACE_INCLUDE_DIRECTORIES ${CONEX_INCLUDE_DIR} + IMPORTED_LOCATION ${CONEX_PREFIX}/lib/${CMAKE_SYSTEM_NAME}/libCONEXsibyll.a + IMPORTED_NO_SONAME 1 + SKIP_BUILD_RPATH FALSE + IMPORTED_LINK_INTERFACE_LIBRARIES ${boost_lib} + INTERFACE_INCLUDE_DIRECTORIES + $<BUILD_INTERFACE:${CONEX_INCLUDE_DIR}> ) set (CONEX_FOUND 1 PARENT_SCOPE) @@ -283,6 +289,7 @@ else () IMPORTED_LOCATION ${CONEX_PREFIX}/lib/${CMAKE_SYSTEM_NAME}/libCONEXsibyll.a IMPORTED_NO_SONAME 1 SKIP_BUILD_RPATH FALSE + IMPORTED_LINK_INTERFACE_LIBRARIES ${boost_lib} INTERFACE_INCLUDE_DIRECTORIES $<BUILD_INTERFACE:${CONEX_INCLUDE_DIR}> )