added sophia interaction model

corsika/detail/modules/sophia/InteractionModel.inl

+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+
+#include <corsika/modules/sophia/ParticleConversion.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/modules/sophia/SophiaStack.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <sophia.hpp>
+
+namespace corsika::sophia {
+
+  inline void InteractionModel::setVerbose(bool const flag) { sophia_listing_ = flag; }
+
+  inline InteractionModel::InteractionModel()
+      : sophia_listing_(false) {
+    // set all particles stable in SOPHIA
+    for (int i = 0; i < 49; ++i) so_csydec_.idb[i] = -abs(so_csydec_.idb[i]);
+  };
+
+  inline InteractionModel::~InteractionModel() {
+    CORSIKA_LOG_DEBUG("Sophia::Model n={}", count_);
+  }
+
+  inline bool constexpr InteractionModel::isValid(Code const projectileId,
+                                                  Code const targetId,
+                                                  HEPEnergyType const sqrtSnn) const {
+    if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
+
+    if (!(targetId == Code::Proton || targetId == Code::Neutron ||
+          targetId == Code::Hydrogen))
+      return false;
+
+    if (projectileId != Code::Photon) return false;
+
+    return true;
+  }
+
+  template <typename TSecondaryView>
+  inline void InteractionModel::doInteraction(TSecondaryView& secondaries,
+                                              Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "projectile: Id={}, E={} GeV, p3={} GeV", projectileId,
+                         projectileP4.getTimeLikeComponent() / 1_GeV,
+                         projectileP4.getSpaceLikeComponents().getComponents() / 1_GeV);
+    CORSIKA_LOGGER_DEBUG(logger_, "target: Id={}, E={} GeV, p3={} GeV", targetId,
+                         targetP4.getTimeLikeComponent() / 1_GeV,
+                         targetP4.getSpaceLikeComponents().getComponents() / 1_GeV);
+
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+
+    CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={}GeV", sqrtS / 1_GeV);
+
+    // accepts only photon-nucleon interactions
+    if (!isValid(projectileId, targetId, sqrtS)) {
+      CORSIKA_LOGGER_ERROR(logger_, "Invalid target/projectile/energy combination");
+      throw std::runtime_error("SOPHIA: Invalid target/projectile/energy combination");
+    }
+
+    auto const nucleonId = (projectileId == Code::Photon ? targetId : projectileId);
+    auto const nucleonP4 = (projectileId == Code::Photon ? targetP4 : projectileP4);
+    auto const photonP4 = (projectileId == Code::Photon ? projectileP4 : targetP4);
+    COMBoost const boost(projectileP4, targetP4);
+
+    int nucleonSophiaCode = convertToSophiaRaw(nucleonId); // either proton or neutron
+    // initialize resonance spectrum
+    initial_(nucleonSophiaCode);
+    double Enucleon = nucleonP4.getTimeLikeComponent() / 1_GeV;
+    double Ephoton = photonP4.getTimeLikeComponent() / 1_GeV;
+    double theta = 0.0; // set head on collision
+    int Imode;
+    CORSIKA_LOGGER_DEBUG(logger_,
+                        "calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}",
+                        nucleonSophiaCode, Enucleon, Ephoton, theta);
+    count_++;
+    // call sophia
+    eventgen_(nucleonSophiaCode, Enucleon, Ephoton, theta, Imode);
+
+    if (sophia_listing_) {
+      int arg = 3;
+      print_event_(arg);
+    }
+
+    auto const& originalCS = boost.getOriginalCS();
+    // SOPHIA has photon along -z  and nucleon along +z (GZK calc..)
+    COMBoost const boostInternal(nucleonP4, photonP4);
+    auto const& csPrime = boost.getRotatedCS();
+    CoordinateSystemPtr csPrimePrime =
+        make_rotation(csPrime, QuantityVector<length_d>{1_m, 0_m, 0_m}, M_PI);
+
+    SophiaStack ss;
+
+    MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType E_final = 0_GeV;
+    for (auto& psop : ss) {
+      // abort on particles that have decayed in SOPHIA. Should not happen!
+      if (psop.hasDecayed()) { // LCOV_EXCL_START
+        throw std::runtime_error("found particle that decayed in SOPHIA!");
+      } // LCOV_EXCL_STOP
+
+      auto momentumSophia = psop.getMomentum(csPrimePrime);
+      momentumSophia.rebase(csPrime);
+      auto const energySophia = psop.getEnergy();
+      auto const P4com = boostInternal.toCoM(FourVector{energySophia, momentumSophia});
+      auto const P4lab = boost.fromCoM(P4com);
+      SophiaCode const pidSophia = psop.getPID();
+      Code const pid = convertFromSophia(pidSophia);
+      auto momentum = P4lab.getSpaceLikeComponents();
+      momentum.rebase(originalCS);
+      HEPEnergyType const Ekin =
+          calculate_kinetic_energy(momentum.getNorm(), get_mass(pid));
+
+      CORSIKA_LOGGER_TRACE(logger_, "SOPHIA: pid={}, p={} GeV", pidSophia,
+                          momentumSophia.getComponents() / 1_GeV);
+
+      CORSIKA_LOGGER_TRACE(logger_, "CORSIKA: pid={}, p={} GeV", pid,
+                          momentum.getComponents() / 1_GeV);
+
+      auto pnew =
+          secondaries.addSecondary(std::make_tuple(pid, Ekin, momentum.normalized()));
+
+      P_final += pnew.getMomentum();
+      E_final += pnew.getEnergy();
+    }
+    CORSIKA_LOGGER_TRACE(logger_, "Efinal={} GeV,Pfinal={} GeV", E_final / 1_GeV,
+                        P_final.getComponents() / 1_GeV);
+  }
+
+} // namespace corsika::sophia
\ No newline at end of file

corsika/modules/Sophia.hpp

@@ -28,10 +28,10 @@ namespace corsika::sophia {
    * The sophia::InteractionModel is wrapped as an InteractionProcess here in order
    * to provide all the functions for ProcessSequence.
    */
-  struct Interaction : public InteractionModel, public InteractionProcess<Interaction> {
-    template <typename TEnvironment>
-    Interaction(TEnvironment const& env)
-        : InteractionModel{env} {}
-  };
+  // struct Interaction : public InteractionModel, public InteractionProcess<Interaction> {
+  //   template <typename TEnvironment>
+  //   Interaction(TEnvironment const& env)
+  //       : InteractionModel{env} {}
+  // };
 
 } // namespace corsika::sophia

corsika/modules/sophia/InteractionModel.hpp

+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+
+#include <tuple>
+
+namespace corsika::sophia {
+
+  /**
+   * @brief Provides the SOPHIA photon-nucleon interaction model.
+   *
+   * This is a TModel argument for InteractionProcess<TModel>.
+   */
+
+  class InteractionModel {
+
+  public:
+    InteractionModel();
+    ~InteractionModel();
+
+    /**
+     * @brief Set the Verbose flag.
+     *
+     * If flag is true, SOPHIA will printout additional secondary particle information
+     * lists, etc.
+     *
+     * @param flag to switch.
+     */
+    void setVerbose(bool const flag);
+
+    /**
+     * @brief evaluated validity of collision system.
+     *
+     * SOPHIA only accepts nucleons as targets, or protons aka Hydrogen or
+     * neutrons (p,n == nucleon).
+     */
+    bool constexpr isValid(Code const projectileId, Code const targetId,
+                           HEPEnergyType const sqrtSnn) const;
+
+    /**
+     * Returns inelastic (production) cross section.
+     *
+     * This cross section must correspond to the process described in doInteraction.
+     * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4: four-momentum of projectile
+     * @param targetP4: four-momentum of target
+     *
+     * @return inelastic cross section
+     * elastic cross section
+     */
+    CrossSectionType getCrossSection(Code const projectile, Code const target,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const {
+      return CrossSectionType::zero();
+    }
+    /**
+     * In this function SOPHIA is called to produce one event. The
+     * event is copied (and boosted) into the frame of the incoming particles.
+     *
+     * @param view is the stack object for the secondaries
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4: four-momentum of projectile
+     * @param targetP4: four-momentum of target
+     */
+
+    template <typename TSecondaries>
+    void doInteraction(TSecondaries& view, Code const projectile, Code const target,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+  private:
+    HEPEnergyType constexpr getMinEnergyCoM() const { return minEnergyCoM_; }
+    HEPEnergyType constexpr getMaxEnergyCoM() const { return maxEnergyCoM_; }
+
+    // hard model limits
+    static HEPEnergyType constexpr minEnergyCoM_ = 1.079166345 * 1e9 * electronvolt;
+    static HEPEnergyType constexpr maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
+
+    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sophia");
+
+    // data members
+    int count_ = 0;
+    bool sophia_listing_;
+
+    std::shared_ptr<spdlog::logger> logger_ =
+        get_logger("corsika_sophia_InteractionModel");
+  };
+
+} // namespace corsika::sophia
+
+#include <corsika/detail/modules/sophia/InteractionModel.inl>
\ No newline at end of file

modules/sophia/sophia.hpp

@@ -76,11 +76,14 @@ extern struct {
   double am2[49];
 } so_mass1_;
 
+// sophia initialization
+void initial_(const int&);
+
 // sophia main subroutine
-void eventgen_(const int&, const double&, const double&, const double&, const int&);
+void eventgen_(const int&, const double&, const double&, const double&, int&);
 
 // print event
-void print_event_(int&);
+void print_event_(const int&);
 
 // decay routine (LA,P0,ND,LL,P)
 // void decpar_(const int&, const double*, int&, int*, double*);

tests/modules/testSophia.cpp

@@ -6,7 +6,7 @@
  * the license.
  */
 
-//#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/Sophia.hpp>
 #include <corsika/modules/sophia/ParticleConversion.hpp>
 
 #include <corsika/framework/core/ParticleProperties.hpp>
@@ -42,12 +42,13 @@ TEST_CASE("Sophia", "modules") {
   SECTION("Sophia -> Corsika") {
     CHECK(Code::Electron ==
           corsika::sophia::convertFromSophia(corsika::sophia::SophiaCode::Electron));
+    CHECK_THROWS(convertFromSophia(corsika::sophia::SophiaCode::Unknown));
   }
 
   SECTION("Corsika -> Sophia") {
     CHECK(corsika::sophia::convertToSophia(Electron::code) ==
           corsika::sophia::SophiaCode::Electron);
-    CHECK(corsika::sophia::convertToSophiaRaw(Proton::code) == 13);    
+    CHECK(corsika::sophia::convertToSophiaRaw(Proton::code) == 13);
   }
 
   SECTION("canInteractInSophia") {
@@ -56,7 +57,7 @@ TEST_CASE("Sophia", "modules") {
     CHECK_FALSE(corsika::sophia::canInteract(Code::XiCPlus));
 
     CHECK_FALSE(corsika::sophia::canInteract(Code::Electron));
-    
+
     CHECK_FALSE(corsika::sophia::canInteract(Code::Iron));
     CHECK_FALSE(corsika::sophia::canInteract(Code::Helium));
   }
@@ -74,3 +75,71 @@ TEST_CASE("Sophia", "modules") {
     CHECK_THROWS(corsika::sophia::getSophiaMass(Code::H0));
   }
 }
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+
+#include <SetupTestEnvironment.hpp>
+#include <SetupTestStack.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+
+template <typename TStackView>
+auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
+  Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
+  for (auto const& p : view) { sum += p.getMomentum(); }
+  return sum;
+}
+
+TEST_CASE("SophiaInterface", "modules") {
+
+  logging::set_level(logging::level::debug);
+
+  // the environment and stack should eventually disappear from here
+  auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
+  auto const& cs = *csPtr;
+  { [[maybe_unused]] auto const& env_dummy = env; }
+
+  auto [stack, viewPtr] = setup::testing::setup_stack(
+      Code::Photon, 10_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, cs);
+  test::StackView& view = *viewPtr;
+
+  RNGManager<>::getInstance().registerRandomStream("sophia");
+
+  SECTION("InteractionInterface - valid targets/projectiles") {
+
+    corsika::sophia::InteractionModel model;
+
+    CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 100_GeV));
+    CHECK(model.isValid(Code::Photon, Code::Hydrogen, 3_GeV));
+
+    FourMomentum const aP4(100_GeV, {cs, 99_GeV, 0_GeV, 0_GeV});
+    FourMomentum const bP4(1_GeV, {cs, 0_GeV, 0_GeV, 0_GeV});
+    CHECK(0_mb == model.getCrossSection(Code::Photon, Code::Proton, aP4, bP4));
+  }
+
+  SECTION("InteractionInterface - interaction") {
+    const HEPEnergyType P0 = 1.2_GeV;
+    MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
+    // also print particles after SOPHIA was called
+    corsika::sophia::InteractionModel model;
+    model.setVerbose(true);
+    HEPEnergyType const Elab = P0;
+    FourMomentum const projectileP4(Elab, plab);
+    FourMomentum const nucleonP4(Proton::mass, MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
+    view.clear();
+    model.doInteraction(view, Code::Photon, Code::Proton, projectileP4, nucleonP4);
+
+    auto const pSum = sumMomentum(view, cs);
+    CHECK(pSum.getComponents(cs).getX() / P0 == Approx(1).margin(0.05));
+    CHECK(pSum.getComponents(cs).getY() / 1_GeV == Approx(0).margin(1e-3));
+    CHECK(pSum.getComponents(cs).getZ() / 1_GeV == Approx(0).margin(1e-3));
+  }
+}
\ No newline at end of file