diff --git a/corsika/framework/process/DynamicInteractionProcess.hpp b/corsika/framework/process/DynamicInteractionProcess.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6cf76684d3c30f8c5d953c0e42b10959c29c07be
--- /dev/null
+++ b/corsika/framework/process/DynamicInteractionProcess.hpp
@@ -0,0 +1,111 @@
+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/ProcessTraits.hpp>
+
+#include <boost/type_index.hpp>
+
+#include <memory>
+
+namespace corsika {
+
+ /**
+ * This class allows selecting/using different InteractionProcesses at runtime without
+ * recompiling the process sequence. The implementation is based on the "type-erasure"
+ * technique.
+ *
+ * @tparam TStack the stack type; has to match the one used by Cascade together with
+ * the ProcessSequence containing the DynamicInteractionProcess.
+ */
+ template <typename TStack>
+ class DynamicInteractionProcess
+ : InteractionProcess<DynamicInteractionProcess<TStack>> {
+ public:
+ using stack_view_type = typename TStack::stack_view_type;
+
+ DynamicInteractionProcess() = default;
+
+ /**
+ * Create new DynamicInteractionProcess. Calls to this instance will be forwared
+ * to the process referred to by obj. The newly created DynamicInteractionProcess
+ * shares ownership of the underlying process, so that you do not need to worry
+ * about it going out of scope.
+ */
+ template <typename TInteractionProcess>
+ DynamicInteractionProcess(std::shared_ptr<TInteractionProcess> obj)
+ : concept_{std::make_unique<ConcreteModel<TInteractionProcess>>(std::move(obj))} {
+ CORSIKA_LOG_DEBUG("creating DynamicInteractionProcess from {}",
+ boost::typeindex::type_id<TInteractionProcess>().pretty_name());
+
+ // poor man's check while we don't have C++20 concepts
+ static_assert(is_interaction_process_v<TInteractionProcess>);
+
+ // since we only forward interaction process API, all other types of corsika
+ // processes are prohibited
+ static_assert(!is_continuous_process_v<TInteractionProcess>);
+ static_assert(!is_decay_process_v<TInteractionProcess>);
+ static_assert(!is_stack_process_v<TInteractionProcess>);
+ static_assert(!is_cascade_equations_process_v<TInteractionProcess>);
+ static_assert(!is_secondaries_process_v<TInteractionProcess>);
+ static_assert(!is_boundary_process_v<TInteractionProcess>);
+ static_assert(!is_cascade_equations_process_v<TInteractionProcess>);
+ }
+
+ /// forwards arguments to doInteraction() of wrapped instance
+ void doInteraction(stack_view_type& view, Code projectileId, Code targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4) {
+ return concept_->doInteraction(view, projectileId, targetId, projectileP4,
+ targetP4);
+ }
+
+ /// forwards arguments to getCrossSection() of wrapped instance
+ CrossSectionType getCrossSection(Code projectileId, Code targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
+ return concept_->getCrossSection(projectileId, targetId, projectileP4, targetP4);
+ }
+
+ private:
+ struct IInteractionModel {
+ virtual ~IInteractionModel() = default;
+ virtual void doInteraction(stack_view_type&, Code, Code, FourMomentum const&,
+ FourMomentum const&) = 0;
+ virtual CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
+ FourMomentum const&) const = 0;
+ };
+
+ template <typename TModel>
+ struct ConcreteModel final : IInteractionModel {
+ ConcreteModel(std::shared_ptr<TModel> obj) noexcept
+ : model_{std::move(obj)} {}
+
+ void doInteraction(stack_view_type& view, Code projectileId, Code targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) override {
+ return model_->doInteraction(view, projectileId, targetId, projectileP4,
+ targetP4);
+ }
+
+ CrossSectionType getCrossSection(Code projectileId, Code targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const override {
+ return model_->getCrossSection(projectileId, targetId, projectileP4, targetP4);
+ }
+
+ private:
+ std::shared_ptr<TModel> model_;
+ };
+
+ std::unique_ptr<IInteractionModel> concept_;
+ };
+} // namespace corsika
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index a7e9dca6b2666b4dcda9176eaadba0575af135f1..48bedadae945365a92cbab82f69af5055de500bb 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -19,6 +19,7 @@
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/process/DynamicInteractionProcess.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
@@ -65,8 +66,10 @@
#include <CLI/Config.hpp>
#include <CLI/Formatter.hpp>
+#include <cstdlib>
#include <iomanip>
#include <limits>
+#include <string>
#include <string_view>
/*
@@ -165,6 +168,10 @@ int main(int argc, char** argv) {
->default_val("info")
->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
->group("Misc.");
+ app.add_option("-M,--hadronModel", "High-energy hadronic interaction model")
+ ->default_val("SIBYLL-2.3d")
+ ->check(CLI::IsMember({"SIBYLL-2.3d", "QGSJet-II.04", "EPOS-LHC"}))
+ ->group("Misc.");
// parse the command line options into the variables
CLI11_PARSE(app, argc, argv);
@@ -287,8 +294,26 @@ int main(int argc, char** argv) {
TrackWriter tracks;
output.add("tracks", tracks);
- corsika::sibyll::Interaction sibyll{env};
- InteractionCounter sibyllCounted{sibyll};
+ DynamicInteractionProcess<setup::Stack<EnvType>> heModel;
+
+ // have SIBYLL always for PROPOSAL photo-hadronic interactions
+ auto sibyll = std::make_shared<corsika::sibyll::Interaction>(env);
+
+ if (auto const modelStr = app["--hadronModel"]->as<std::string_view>();
+ modelStr == "SIBYLL-2.3d") {
+ heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{sibyll};
+ } else if (modelStr == "QGSJet-II.04") {
+ heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{
+ std::make_shared<corsika::qgsjetII::Interaction>()};
+ } else if (modelStr == "EPOS-LHC") {
+ heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{
+ std::make_shared<corsika::epos::Interaction>()};
+ } else {
+ CORSIKA_LOG_CRITICAL("invalid choice \"{}\"; also check argument parser", modelStr);
+ return EXIT_FAILURE;
+ }
+
+ InteractionCounter heCounted{heModel};
corsika::pythia8::Decay decayPythia;
@@ -327,7 +352,7 @@ int main(int argc, char** argv) {
HEPEnergyType heHadronModelThreshold = 63.1_GeV;
corsika::proposal::Interaction emCascade(
- env, sophia, sibyll.getHadronInteractionModel(), heHadronModelThreshold);
+ env, sophia, sibyll->getHadronInteractionModel(), heHadronModelThreshold);
// use BetheBlochPDG for hadronic continuous losses, and proposal otherwise
corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuousProposal(
@@ -356,11 +381,9 @@ int main(int argc, char** argv) {
bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
};
auto hadronSequence =
- make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, sibyllCounted);
+ make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, heCounted);
auto decaySequence = make_sequence(decayPythia, decaySibyll);
- TrackWriter trackWriter{tracks};
-
// observation plane
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
@@ -430,9 +453,7 @@ int main(int argc, char** argv) {
Efinal / 1_GeV, dEdX.getEnergyLost() / 1_GeV,
observationLevel.getEnergyGround() / 1_GeV, (Efinal / E0 - 1) * 100);
- // auto const hists = heModelCounted.getHistogram() +
- // urqmdCounted.getHistogram();
- auto const hists = sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
+ auto const hists = heCounted.getHistogram() + urqmdCounted.getHistogram();
save_hist(hists.labHist(), labHist_file, true);
save_hist(hists.CMSHist(), cMSHist_file, true);
@@ -440,4 +461,6 @@ int main(int argc, char** argv) {
// and finalize the output on disk
output.endOfLibrary();
+
+ return EXIT_SUCCESS;
}
diff --git a/tests/framework/CMakeLists.txt b/tests/framework/CMakeLists.txt
index caaec9e4cc04763a30ceec719107e57f9b91a578..14670134dba78d4c20cfae3c4d4930bc3eeeb431 100644
--- a/tests/framework/CMakeLists.txt
+++ b/tests/framework/CMakeLists.txt
@@ -4,6 +4,7 @@ set (test_framework_sources
testFourVector.cpp
testSaveBoostHistogram.cpp
testClassTimer.cpp
+ testDynamicInteractionProcess.cpp;
testLogging.cpp
testParticles.cpp
testStackInterface.cpp
diff --git a/tests/framework/testDynamicInteractionProcess.cpp b/tests/framework/testDynamicInteractionProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7c68e341babab575ce04bcf8975c09e455d5ef8f
--- /dev/null
+++ b/tests/framework/testDynamicInteractionProcess.cpp
@@ -0,0 +1,143 @@
+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/framework/process/DynamicInteractionProcess.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+
+struct DummyStack {
+ using stack_view_type = int&;
+};
+
+struct DummyProcessA : public InteractionProcess<DummyProcessA> {
+ int id_;
+
+ DummyProcessA(int id)
+ : id_{id} {}
+
+ // mockup to "identify" the process via its ID
+ CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
+ FourMomentum const&) const {
+ return id_ * 10_mb;
+ }
+
+ // mockup: set argument to our ID
+ template <typename TArgument>
+ void doInteraction(TArgument& argument, Code, Code, FourMomentum const&,
+ FourMomentum const&) {
+ argument = id_ * 10;
+ }
+
+ // to make sure we can test proper handling of dangling references and
+ // stuff like that
+ ~DummyProcessA() { id_ = 0; }
+
+ // prevent copying
+ DummyProcessA(DummyProcessA const&) = delete;
+ DummyProcessA& operator=(DummyProcessA const&) = delete;
+
+ // allow moving
+ DummyProcessA(DummyProcessA&&) = default;
+ DummyProcessA& operator=(DummyProcessA&&) = default;
+};
+
+// same as above to have another class
+struct DummyProcessB : public InteractionProcess<DummyProcessB> {
+ int id_;
+
+ DummyProcessB(int id)
+ : id_{id} {}
+
+ CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
+ FourMomentum const&) const {
+ return id_ * 1_mb;
+ }
+
+ template <typename TArgument>
+ void doInteraction(TArgument& argument, Code, Code, FourMomentum const&,
+ FourMomentum const&) {
+ argument = id_;
+ }
+
+ // to make sure we can test proper handling of dangling references and
+ // stuff like that
+ ~DummyProcessB() { id_ = 0; }
+
+ // prevent copying
+ DummyProcessB(DummyProcessB const&) = delete;
+ DummyProcessB& operator=(DummyProcessB const&) = delete;
+
+ // allow moving
+ DummyProcessB(DummyProcessB&&) = default;
+ DummyProcessB& operator=(DummyProcessB&&) = default;
+};
+
+TEST_CASE("DynamicInteractionProcess", "[process]") {
+ auto const rootCS = get_root_CoordinateSystem();
+ FourMomentum const fourVec{1_GeV, {rootCS, 1_GeV, 1_GeV, 1_GeV}};
+
+ std::shared_ptr<DummyProcessA> a = std::make_shared<DummyProcessA>(1);
+ std::shared_ptr<DummyProcessB> b = std::make_shared<DummyProcessB>(2);
+
+ SECTION("getCrossSection") {
+ DynamicInteractionProcess<DummyStack> dynamic;
+
+ dynamic = a;
+ REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+ 10_mb ==
+ Approx(1).epsilon(1e-6));
+
+ dynamic = b;
+ REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+ 2_mb ==
+ Approx(1).epsilon(1e-6));
+
+ dynamic = DynamicInteractionProcess<DummyStack>{std::make_shared<DummyProcessA>(3)};
+ REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+ 30_mb ==
+ Approx(1).epsilon(1e-6));
+
+ dynamic = std::make_shared<DummyProcessB>(4);
+ REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+ 4_mb ==
+ Approx(1).epsilon(1e-6));
+
+ // test process going out of scope
+ { dynamic = std::make_shared<DummyProcessB>(5); }
+ REQUIRE(dynamic.getCrossSection(Code::Electron, Code::Electron, fourVec, fourVec) /
+ 5_mb ==
+ Approx(1).epsilon(1e-6));
+ }
+
+ SECTION("doInteraction") {
+ int value{};
+
+ DynamicInteractionProcess<DummyStack> dynamic;
+
+ dynamic = a;
+ dynamic.doInteraction(value, Code::Electron, Code::Electron, fourVec, fourVec);
+ REQUIRE(value == 10);
+
+ dynamic = b;
+ dynamic.doInteraction(value, Code::Electron, Code::Electron, fourVec, fourVec);
+ REQUIRE(value == 2);
+
+ dynamic = DynamicInteractionProcess<DummyStack>{std::make_shared<DummyProcessA>(3)};
+ dynamic.doInteraction(value, Code::Electron, Code::Electron, fourVec, fourVec);
+ REQUIRE(value == 30);
+
+ dynamic = std::make_shared<DummyProcessB>(4);
+ dynamic.doInteraction(value, Code::Electron, Code::Electron, fourVec, fourVec);
+ REQUIRE(value == 4);
+ }
+}