From c777957e664ff02ae59bf049aa3e4f85bebda914 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 1 Apr 2022 20:48:23 +0100
Subject: [PATCH] first try photohadronic

---
 .../detail/modules/proposal/Interaction.inl   | 76 ++++++++++++++++---
 corsika/modules/proposal/Interaction.hpp      | 11 ++-
 examples/em_shower.cpp                        |  8 +-
 modules/data                                  |  2 +-
 src/modules/sibyll/sibyll_codes.dat           |  2 +-
 tests/modules/testSibyll.cpp                  |  1 +
 6 files changed, 81 insertions(+), 19 deletions(-)

diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index 062aebe19..77833c1d8 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -8,7 +8,6 @@
 
 #include <corsika/media/IMediumModel.hpp>
 #include <corsika/media/NuclearComposition.hpp>
-#include <corsika/modules/proposal/Interaction.hpp>
 #include <corsika/framework/utility/COMBoost.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
@@ -20,8 +19,10 @@
 namespace corsika::proposal {
 
   template <typename TEnvironment>
-  inline Interaction::Interaction(TEnvironment const& _env)
-      : ProposalProcessBase(_env) {}
+  inline Interaction::Interaction(TEnvironment const& _env,
+                                  corsika::sibyll::Interaction& hadint)
+      : ProposalProcessBase(_env)
+      , hadronicInteraction_(hadint) {}
 
   inline void Interaction::buildCalculator(Code code, NuclearComposition const& comp) {
     // search crosssection builder for given particle
@@ -99,15 +100,66 @@ namespace corsika::proposal {
       PROPOSAL::Component target;
       if (type != PROPOSAL::InteractionType::Ioniz)
         target = PROPOSAL::Component::GetComponentForHash(target_hash);
-      auto sec =
-          std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
-      for (auto& s : sec) {
-        auto E = s.energy * 1_MeV;
-        auto vecProposal = s.direction;
-        auto dir = DirectionVector(
-            labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
-        auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
-        view.addSecondary(std::make_tuple(sec_code, E - get_mass(sec_code), dir));
+      if (type == PROPOSAL::InteractionType::Photonuclear)
+        CORSIKA_LOG_INFO(
+            "HE photo-hadronic interaction! calling hadronic interaction model. Energy = "
+            "{} GeV, v = {}, v * E = {}",
+            projectile.getEnergy() / 1_GeV, v, v * projectile.getEnergy());
+      if (type == PROPOSAL::InteractionType::Photonuclear &&
+          projectile.getEnergy() > heHadronicModelThresholdLab_) {
+        CORSIKA_LOG_INFO(
+            "HE photo-hadronic interaction! calling hadronic interaction model..");
+
+        //  copy from NuclearInteractionModel
+        //  temporarily add to stack, will be removed after interaction in DoInteraction
+        typename TStackView::inner_stack_value_type photonStack;
+        Point const pDummy(labCS, {0_m, 0_m, 0_m});
+        TimeType const tDummy = 0_ns;
+        Code const hadPhotonCode = Code::Rho0; // stand in for hadronic-photon
+        auto const A = target.GetAtomicNum();
+        auto const Z = target.GetNucCharge();
+        Code const targetId = get_nucleus_code(A, Z);
+        // target at rest
+        FourMomentum const targetP4(get_mass(targetId),
+                                    MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
+        auto const p3PhotonLab = projectileP4.getSpaceLikeComponents();
+        HEPEnergyType const Ekin = projectileP4.getTimeLikeComponent();
+        auto hadronicPhoton = photonStack.addParticle(std::make_tuple(
+            hadPhotonCode, Ekin, p3PhotonLab.normalized(), pDummy, tDummy));
+        hadronicPhoton.setNode(view.getProjectile().getNode());
+        CORSIKA_LOG_INFO("gam - had interaction: kin. energy of gamma: {} GeV",
+                         Ekin / 1_GeV);
+
+        // create inelastic interaction of the hadronic photon
+        // create new StackView for the photon
+        TStackView photon_secondaries(hadronicPhoton);
+
+        // call inner hadronic event generator
+        CORSIKA_LOG_TRACE("calling HadronicInteraction...");
+        CORSIKA_LOG_INFO("{} + {} interactions. Ekinlab = {}", hadPhotonCode, targetId,
+                         Ekin / 1_GeV);
+        hadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
+                                           projectileP4, targetP4);
+        for (const auto& pSec : photon_secondaries) {
+
+          auto const p3lab = pSec.getMomentum();
+          Code const pid = pSec.getPID();
+          HEPEnergyType const mass = get_mass(pid);
+          HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+          view.addSecondary(std::make_tuple(pid, Ekin, p3lab.normalized()));
+        }
+        // throw std::runtime_error("photo-hadronic interaction!");
+      } else {
+        auto sec =
+            std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
+        for (auto& s : sec) {
+          auto E = s.energy * 1_MeV;
+          auto vecProposal = s.direction;
+          auto dir = DirectionVector(
+              labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
+          auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
+          view.addSecondary(std::make_tuple(sec_code, E - get_mass(sec_code), dir));
+        }
       }
     }
     return ProcessReturn::Ok;
diff --git a/corsika/modules/proposal/Interaction.hpp b/corsika/modules/proposal/Interaction.hpp
index 3b769bc34..569d5d9c9 100644
--- a/corsika/modules/proposal/Interaction.hpp
+++ b/corsika/modules/proposal/Interaction.hpp
@@ -16,7 +16,7 @@
 #include <corsika/framework/geometry/FourVector.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/random/UniformRealDistribution.hpp>
-
+#include <corsika/modules/Sibyll.hpp>
 #include <corsika/modules/proposal/ProposalProcessBase.hpp>
 
 namespace corsika::proposal {
@@ -25,7 +25,9 @@ namespace corsika::proposal {
   //! Electro-magnetic and photon stochastic losses produced by proposal. It makes
   //! use of interpolation tables which are runtime intensive calculation, but can be
   //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
+  //! @tparam THadronModel
   //!
+
   class Interaction : public InteractionProcess<Interaction>, ProposalProcessBase {
 
     enum { eSECONDARIES, eINTERACTION };
@@ -46,7 +48,7 @@ namespace corsika::proposal {
     //! compositions and stochastic description limited by the particle cut.
     //!
     template <typename TEnvironment>
-    Interaction(TEnvironment const& env);
+    Interaction(TEnvironment const& env, corsika::sibyll::Interaction&);
 
     //!
     //! Calculate the rates for the different targets and interactions. Sample a
@@ -63,6 +65,11 @@ namespace corsika::proposal {
     template <typename TParticle>
     CrossSectionType getCrossSection(TParticle const& p, Code const projectileId,
                                      FourMomentum const& projectileP4);
+
+  private:
+    corsika::sibyll::Interaction& hadronicInteraction_;
+    static HEPEnergyType constexpr heHadronicModelThresholdLab_ =
+        80. * 1e9 * electronvolt;
   };
 } // namespace corsika::proposal
 
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 018666a46..f4e2c37a4 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -38,6 +38,7 @@
 #include <corsika/modules/ObservationPlane.hpp>
 #include <corsika/modules/ParticleCut.hpp>
 #include <corsika/modules/TrackWriter.hpp>
+#include <corsika/modules/Sibyll.hpp>
 #include <corsika/modules/PROPOSAL.hpp>
 
 #include <corsika/setup/SetupStack.hpp>
@@ -57,8 +58,7 @@
   executable. If you include the header below multiple times and
   link this togehter, it will fail.
  */
-#include <corsika/modules/sibyll/Random.hpp>
-#include <corsika/modules/urqmd/Random.hpp>
+#include <corsika/modules/Random.hpp>
 
 using namespace corsika;
 using namespace std;
@@ -66,6 +66,7 @@ using namespace std;
 void registerRandomStreams(int seed) {
   RNGManager<>::getInstance().registerRandomStream("cascade");
   RNGManager<>::getInstance().registerRandomStream("proposal");
+  RNGManager<>::getInstance().registerRandomStream("sibyll");
   if (seed == 0) {
     std::random_device rd;
     seed = rd();
@@ -162,7 +163,8 @@ int main(int argc, char** argv) {
 
   ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, 60_GeV, 100_PeV, 100_PeV, true,
                                              dEdX);
-  corsika::proposal::Interaction emCascade(env);
+  corsika::sibyll::Interaction sibyll;
+  corsika::proposal::Interaction emCascade(env, sibyll);
   corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
   //  BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
 
diff --git a/modules/data b/modules/data
index eded83360..c42a119ea 160000
--- a/modules/data
+++ b/modules/data
@@ -1 +1 @@
-Subproject commit eded833608bbc6845f9d7e2af4e151683f27818f
+Subproject commit c42a119ea1dc5053cb748994a9d8a4561f2206d5
diff --git a/src/modules/sibyll/sibyll_codes.dat b/src/modules/sibyll/sibyll_codes.dat
index a445af00d..9b44025ad 100644
--- a/src/modules/sibyll/sibyll_codes.dat
+++ b/src/modules/sibyll/sibyll_codes.dat
@@ -22,7 +22,7 @@ NuTauBar                 93      0       CannotInteract
 Photon                   1       0       CannotInteract
 Pi0                      6       1       Pion
 #                        rho0    could interact but sibyll has no cross section/interaction length. was used for gamma had int
-Rho0                     27      0       CannotInteract
+Rho0                     27      1       Pion
 K0Long                   11      1       Kaon
 K0                       21      0       Kaon
 K0Bar                    22      0       Kaon
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index a4dbdd4a1..dffba30ca 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -130,6 +130,7 @@ TEST_CASE("SibyllInterface", "modules") {
     CHECK_FALSE(model.isValid(Code::Proton, Code::Helium3, 100_GeV));
     CHECK_FALSE(model.isValid(Code::Proton, Code::Iron, 100_GeV));
     CHECK(model.isValid(Code::Proton, Code::Oxygen, 100_GeV));
+    CHECK(model.isValid(Code::Rho0, Code::Oxygen, 100_GeV));
     // beam particles
     CHECK_FALSE(model.isValid(Code::Electron, Code::Oxygen, 100_GeV));
     CHECK_FALSE(model.isValid(Code::Iron, Code::Oxygen, 100_GeV));
-- 
GitLab