Interaction model selection at run-time

corsika/framework/process/DynamicInteractionProcess.hpp

+/*
+ * (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

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;
 }

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

tests/framework/testDynamicInteractionProcess.cpp

+/*
+ * (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);
+  }
+}