From d48da02c85139e53d4dc2ac16fec49b05d843600 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 5 Apr 2022 18:39:10 +0100
Subject: [PATCH] added unit test for photo-hadronic in proposal

---
 .../detail/modules/proposal/Interaction.inl   | 93 ++++++++++---------
 .../modules/proposal/ContinuousProcess.hpp    |  1 +
 corsika/modules/proposal/Interaction.hpp      |  9 ++
 tests/modules/CMakeLists.txt                  |  1 +
 tests/modules/testProposal.cpp                | 61 ++++++++++++
 5 files changed, 123 insertions(+), 42 deletions(-)
 create mode 100644 tests/modules/testProposal.cpp

diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index 77833c1d8..ef75b79c9 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -107,48 +107,9 @@ namespace corsika::proposal {
             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!");
+        Code const targetId =
+            get_nucleus_code(target.GetAtomicNum(), target.GetNucCharge());
+        doHadronicInteraction(view, labCS, projectileP4, targetId);
       } else {
         auto sec =
             std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
@@ -165,6 +126,54 @@ namespace corsika::proposal {
     return ProcessReturn::Ok;
   }
 
+  template <typename TStackView>
+  inline ProcessReturn Interaction::doHadronicInteraction(
+      TStackView& view, CoordinateSystemPtr const& labCS,
+      FourMomentum const& projectileP4, Code const targetId) {
+    CORSIKA_LOG_INFO(
+        "HE photo-hadronic interaction! calling hadronic interaction model..");
+
+    //  copy from sibyll::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!");
+    return ProcessReturn::Ok;
+  }
+
   template <typename TParticle>
   inline CrossSectionType Interaction::getCrossSection(TParticle const& projectile,
                                                        Code const projectileId,
diff --git a/corsika/modules/proposal/ContinuousProcess.hpp b/corsika/modules/proposal/ContinuousProcess.hpp
index ec2532017..d7b0e49ec 100644
--- a/corsika/modules/proposal/ContinuousProcess.hpp
+++ b/corsika/modules/proposal/ContinuousProcess.hpp
@@ -17,6 +17,7 @@
 #include <corsika/framework/process/ProcessReturn.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/random/UniformRealDistribution.hpp>
+#include <corsika/modules/writers/WriterOff.hpp>
 
 #include <corsika/modules/proposal/ProposalProcessBase.hpp>
 
diff --git a/corsika/modules/proposal/Interaction.hpp b/corsika/modules/proposal/Interaction.hpp
index 569d5d9c9..d42970f66 100644
--- a/corsika/modules/proposal/Interaction.hpp
+++ b/corsika/modules/proposal/Interaction.hpp
@@ -59,6 +59,15 @@ namespace corsika::proposal {
     ProcessReturn doInteraction(TSecondaryView&, Code const projectileId,
                                 FourMomentum const& projectileP4);
 
+    //!
+    //! Calculate produce the hadronic secondaries in a hadronic photon interaction and
+    //! store them on the particle stack.
+    //!
+    template <typename TSecondaryView>
+    ProcessReturn doHadronicInteraction(TSecondaryView&, CoordinateSystemPtr const&,
+                                        FourMomentum const& projectileP4,
+                                        Code const targetId);
+
     //!
     //! Calculates and returns the cross section.
     //!
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 00f28c41c..2073a7d50 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -5,6 +5,7 @@ set (test_modules_sources
   ## testExecTime.cpp --> needs to be fixed, see #326
   testObservationPlane.cpp
   testQGSJetII.cpp
+  testProposal.cpp
   testPythia8.cpp
   testUrQMD.cpp
   testCONEX.cpp
diff --git a/tests/modules/testProposal.cpp b/tests/modules/testProposal.cpp
new file mode 100644
index 000000000..07d579699
--- /dev/null
+++ b/tests/modules/testProposal.cpp
@@ -0,0 +1,61 @@
+/*
+ * (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/modules/PROPOSAL.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+
+#include <SetupTestEnvironment.hpp>
+#include <SetupTestStack.hpp>
+#include <catch2/catch.hpp>
+#include <tuple>
+
+using namespace corsika;
+using namespace corsika::proposal;
+
+using DummyEnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using DummyEnvironment = Environment<DummyEnvironmentInterface>;
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+
+TEST_CASE("ProposalInterface", "modules") {
+
+  logging::set_level(logging::level::info);
+
+  // the test environment
+  auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
+  auto const& cs = *csPtr;
+
+  auto [stackPtr, viewPtr] = setup::testing::setup_stack(
+      Code::Electron, 10_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, cs);
+  test::StackView& view = *viewPtr;
+
+  RNGManager<>::getInstance().registerRandomStream("sibyll");
+  RNGManager<>::getInstance().registerRandomStream("proposal");
+
+  corsika::sibyll::Interaction hadModel;
+  corsika::proposal::Interaction emModel(*env, hadModel);
+
+  SECTION("InteractionInterface - cross section") {
+    auto& stack = *stackPtr;
+    auto particle = stack.first();
+    FourMomentum P4(
+        100_GeV,
+        {cs, {sqrt(static_pow<2>(100_MeV) - static_pow<2>(Proton::mass)), 0_eV, 0_eV}});
+    CHECK(emModel.getCrossSection(particle, Code::Proton, P4) == 0_mb);
+  }
+
+  SECTION("InteractionInterface - hadronic photon interaction") {
+    // auto& stack = *stackPtr;
+    // auto particle = stack.first();
+    FourMomentum P4(100_TeV, {cs, {100_TeV, 0_eV, 0_eV}});
+    CHECK(emModel.doHadronicInteraction(view, cs, P4, Code::Oxygen) == ProcessReturn::Ok);
+  }
+}
\ No newline at end of file
-- 
GitLab