From b02eef1fa2f3b956bfbc30b05c95202e2dc3ef52 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 15 Apr 2022 12:50:02 +0100
Subject: [PATCH] added threshold parameter

---
 .../modules/proposal/HadronicPhotonModel.inl  | 26 +++++++++++++++----
 .../modules/proposal/InteractionModel.inl     | 13 +++++-----
 corsika/modules/PROPOSAL.hpp                  |  4 +--
 .../modules/proposal/HadronicPhotonModel.hpp  |  6 ++---
 corsika/modules/proposal/InteractionModel.hpp |  2 +-
 examples/corsika.cpp                          |  8 ++++--
 examples/em_shower.cpp                        |  3 ++-
 examples/mars.cpp                             |  9 +++++--
 tests/modules/testProposal.cpp                |  7 +++--
 9 files changed, 53 insertions(+), 25 deletions(-)

diff --git a/corsika/detail/modules/proposal/HadronicPhotonModel.inl b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
index 5aa5aded3..c6e30a4d3 100644
--- a/corsika/detail/modules/proposal/HadronicPhotonModel.inl
+++ b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
@@ -16,15 +16,32 @@
 namespace corsika::proposal {
 
   template <typename THadronicModel>
-  inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(THadronicModel& _hadint)
-      : heHadronicInteraction_(_hadint){};
+  inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(
+      THadronicModel& _hadint, HEPEnergyType const& _heenthresholdNN)
+      : heHadronicInteraction_(_hadint)
+      , heHadronicModelThresholdLabNN_(_heenthresholdNN) {
+    // check validity of threshold assuming photon-nucleon
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtS =
+        calculate_com_energy(_heenthresholdNN, Rho0::mass, Proton::mass);
+    if (!heHadronicInteraction_.isValid(Code::Rho0, Code::Proton, sqrtS)) {
+      CORSIKA_LOG_ERROR(
+          "Invalid energy threshold for hadron interaction model. theshold_lab= {} GeV, "
+          "theshold_com={} GeV",
+          _heenthresholdNN / 1_GeV, sqrtS / 1_GeV);
+      throw std::runtime_error("Configuration error!");
+    }
+    CORSIKA_LOG_INFO(
+        "Threshold for HE hadronic interactions in proposal set to Elab={} GeV",
+        _heenthresholdNN / 1_GeV);
+  };
 
   template <typename THadronicModel>
   template <typename TStackView>
   inline ProcessReturn HadronicPhotonModel<THadronicModel>::doHadronicPhotonInteraction(
       TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
       Code const& targetId) {
-    if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLab_) {
+    if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLabNN_) {
       CORSIKA_LOG_INFO(
           "HE photo-hadronic interaction! calling hadronic interaction model..");
 
@@ -46,8 +63,7 @@ namespace corsika::proposal {
       TStackView photon_secondaries(hadronicPhoton);
 
       // call inner hadronic event generator
-      CORSIKA_LOG_TRACE("calling HadronicInteraction...");
-      CORSIKA_LOG_INFO("{} + {} interactions. Ekinlab = {} GeV", hadPhotonCode, targetId,
+      CORSIKA_LOG_INFO("{} + {} interaction. Ekinlab = {} GeV", hadPhotonCode, targetId,
                        photonP4.getTimeLikeComponent() / 1_GeV);
       heHadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
                                            photonP4, targetP4);
diff --git a/corsika/detail/modules/proposal/InteractionModel.inl b/corsika/detail/modules/proposal/InteractionModel.inl
index d1f63c685..867b03858 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -20,10 +20,11 @@ namespace corsika::proposal {
 
   template <typename THadronicModel>
   template <typename TEnvironment>
-  inline InteractionModel<THadronicModel>::InteractionModel(TEnvironment const& _env,
-                                                            THadronicModel& _hadint)
+  inline InteractionModel<THadronicModel>::InteractionModel(
+      TEnvironment const& _env, THadronicModel& _hadint,
+      HEPEnergyType const& _enthreshold)
       : ProposalProcessBase(_env)
-      , HadronicPhotonModel<THadronicModel>(_hadint) {}
+      , HadronicPhotonModel<THadronicModel>(_hadint, _enthreshold) {}
 
   template <typename THadronicModel>
   inline void InteractionModel<THadronicModel>::buildCalculator(
@@ -124,10 +125,8 @@ namespace corsika::proposal {
           auto const Z = int(target.GetNucCharge());
           Code const targetId = get_nucleus_code(A, Z);
           CORSIKA_LOG_INFO(
-              "photo-hadronic interaction! projectile={} target={} energy={} GeV",
-              projectileId,
-              (is_nucleus(targetId) ? get_nucleus_name(targetId) : get_name(targetId)),
-              E / 1_GeV);
+              "photo-hadronic interaction of projectile={} with target={}! Energy={} GeV",
+              projectileId, targetId, E / 1_GeV);
           this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
         } else {
           auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
diff --git a/corsika/modules/PROPOSAL.hpp b/corsika/modules/PROPOSAL.hpp
index a019ac63d..13e18c3de 100644
--- a/corsika/modules/PROPOSAL.hpp
+++ b/corsika/modules/PROPOSAL.hpp
@@ -18,7 +18,7 @@ namespace corsika::proposal {
                       public InteractionProcess<Interaction<THadronicModel>> {
   public:
     template <typename TEnvironment>
-    Interaction(TEnvironment const& env, THadronicModel& model)
-        : InteractionModel<THadronicModel>(env, model) {}
+    Interaction(TEnvironment const& env, THadronicModel& model, HEPEnergyType const& thr)
+        : InteractionModel<THadronicModel>(env, model, thr) {}
   };
 } // namespace corsika::proposal
diff --git a/corsika/modules/proposal/HadronicPhotonModel.hpp b/corsika/modules/proposal/HadronicPhotonModel.hpp
index a3d9b6064..befc2bd1b 100644
--- a/corsika/modules/proposal/HadronicPhotonModel.hpp
+++ b/corsika/modules/proposal/HadronicPhotonModel.hpp
@@ -24,7 +24,7 @@ namespace corsika::proposal {
   template <class THadronicModel>
   class HadronicPhotonModel {
   public:
-    HadronicPhotonModel(THadronicModel&);
+    HadronicPhotonModel(THadronicModel&, HEPEnergyType const&);
     //!
     //! Calculate produce the hadronic secondaries in a hadronic photon interaction and
     //! store them on the particle stack.
@@ -35,8 +35,8 @@ namespace corsika::proposal {
 
   private:
     THadronicModel& heHadronicInteraction_;
-    static HEPEnergyType constexpr heHadronicModelThresholdLab_ =
-        80. * 1e9 * electronvolt; //!< energy threshold between LE and HE model
+    //! threshold for high energy hadronic interaction model. Lab. energy per nucleon
+    HEPEnergyType heHadronicModelThresholdLabNN_;
   };
 } // namespace corsika::proposal
 
diff --git a/corsika/modules/proposal/InteractionModel.hpp b/corsika/modules/proposal/InteractionModel.hpp
index ec5739225..135f685d1 100644
--- a/corsika/modules/proposal/InteractionModel.hpp
+++ b/corsika/modules/proposal/InteractionModel.hpp
@@ -55,7 +55,7 @@ namespace corsika::proposal {
     //! compositions and stochastic description limited by the particle cut.
     //!
     template <typename TEnvironment>
-    InteractionModel(TEnvironment const& env, THadronicModel&);
+    InteractionModel(TEnvironment const& env, THadronicModel&, HEPEnergyType const&);
 
     //!
     //! Calculate the rates for the different targets and interactions. Sample a
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index 68d428ab0..9d4e05656 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -321,7 +321,10 @@ int main(int argc, char** argv) {
   HEPEnergyType const hadcut = 50_GeV;
   ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, dEdX);
 
-  corsika::proposal::Interaction emCascade(env, sibyll);
+  // energy threshold for high energy hadronic model. Affects LE/HE switch for hadron
+  // interactions and the hadronic photon model in proposal
+  HEPEnergyType heHadronModelThreshold = 63.1_GeV;
+  corsika::proposal::Interaction emCascade(env, sibyll, heHadronModelThreshold);
   // NOT available for PROPOSAL due to interface trouble:
   // InteractionCounter emCascadeCounted(emCascade);
   // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env);
@@ -342,7 +345,8 @@ int main(int argc, char** argv) {
         : cutE_(cutE) {}
     bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
   };
-  auto hadronSequence = make_select(EnergySwitch(63.1_GeV), urqmdCounted, heModelCounted);
+  auto hadronSequence =
+      make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, heModelCounted);
   auto decaySequence = make_sequence(decayPythia, decaySibyll);
 
   TrackWriter trackWriter{tracks};
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index f4e2c37a4..aec2e1a3e 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -164,7 +164,8 @@ int main(int argc, char** argv) {
   ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, 60_GeV, 100_PeV, 100_PeV, true,
                                              dEdX);
   corsika::sibyll::Interaction sibyll;
-  corsika::proposal::Interaction emCascade(env, sibyll);
+  HEPEnergyType heThresholdNN = 80_GeV;
+  corsika::proposal::Interaction emCascade(env, sibyll, heThresholdNN);
   corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
   //  BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
 
diff --git a/examples/mars.cpp b/examples/mars.cpp
index 23acc590e..df9e531fa 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -353,7 +353,11 @@ int main(int argc, char** argv) {
 
   // decaySibyll.printDecayConfig();
 
-  corsika::proposal::Interaction emCascade(env, sibyll);
+  // energy threshold for high energy hadronic model. Affects LE/HE switch for hadron
+  // interactions and the hadronic photon model in proposal
+  HEPEnergyType heHadronModelThreshold = 63.1_GeV;
+
+  corsika::proposal::Interaction emCascade(env, sibyll, heHadronModelThreshold);
   // NOT possible right now, due to interface difference for PROPOSAL:
   //  InteractionCounter emCascadeCounted(emCascade);
   // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
@@ -375,7 +379,8 @@ int main(int argc, char** argv) {
         : cutE_(cutE) {}
     bool operator()(Particle const& p) const { return (p.getKineticEnergy() < cutE_); }
   };
-  auto hadronSequence = make_select(EnergySwitch(63.1_GeV), urqmdCounted, heModelCounted);
+  auto hadronSequence =
+      make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, heModelCounted);
   auto decaySequence = make_sequence(decayPythia, decaySibyll);
 
   // track writer
diff --git a/tests/modules/testProposal.cpp b/tests/modules/testProposal.cpp
index d4b3c62dd..7b6f14812 100644
--- a/tests/modules/testProposal.cpp
+++ b/tests/modules/testProposal.cpp
@@ -41,6 +41,9 @@ public:
                           MomentumVector(csPrime, {0_GeV, 0_GeV, 0_GeV}).normalized()));
     }
   }
+  bool constexpr isValid(Code const, Code const, HEPEnergyType const) const {
+    return true;
+  };
 };
 
 TEST_CASE("ProposalInterface", "modules") {
@@ -58,8 +61,8 @@ TEST_CASE("ProposalInterface", "modules") {
   RNGManager<>::getInstance().registerRandomStream("proposal");
 
   DummyHadronicModel hadModel;
-
-  corsika::proposal::InteractionModel emModel(*env, hadModel);
+  HEPEnergyType heThresholdLab = 80_GeV;
+  corsika::proposal::InteractionModel emModel(*env, hadModel, heThresholdLab);
 
   SECTION("InteractionInterface - cross section") {
     auto& stack = *stackPtr;
-- 
GitLab