diff --git a/corsika/detail/modules/sophia/InteractionModel.inl b/corsika/detail/modules/sophia/InteractionModel.inl
new file mode 100644
index 0000000000000000000000000000000000000000..89560a3929451e4716bf5e9604a0df2510ad5764
--- /dev/null
+++ b/corsika/detail/modules/sophia/InteractionModel.inl
@@ -0,0 +1,142 @@
+/*
+ * (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
diff --git a/corsika/modules/Sophia.hpp b/corsika/modules/Sophia.hpp
index 85522208b8973d6571d885af7cb5a0d1bbd9e694..48b4b6706defad2bbd166cab2438ec2b84941343 100644
--- a/corsika/modules/Sophia.hpp
+++ b/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
diff --git a/corsika/modules/sophia/InteractionModel.hpp b/corsika/modules/sophia/InteractionModel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2f6d7636216242ac0c61a46115bd9a2f6a65773
--- /dev/null
+++ b/corsika/modules/sophia/InteractionModel.hpp
@@ -0,0 +1,105 @@
+/*
+ * (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
diff --git a/modules/sophia/sophia.hpp b/modules/sophia/sophia.hpp
index b2564f195cbffca8714bb6f99faf35d363bab952..fbf7e0e64424b8eaaa8f47c8a378e25499a62850 100644
--- a/modules/sophia/sophia.hpp
+++ b/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*);
diff --git a/tests/modules/testSophia.cpp b/tests/modules/testSophia.cpp
index 16170f1dabc86b08c9f8225433b15af9bd8277c0..fc843e653a122e6da4f4824349c2d8bf2859ad40 100644
--- a/tests/modules/testSophia.cpp
+++ b/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