added hadronic photon model

corsika/detail/modules/proposal/InteractionModel.inl

@@ -22,9 +22,9 @@ namespace corsika::proposal {
   template <typename THadronicModel>
   template <typename TEnvironment>
   inline InteractionModel<THadronicModel>::InteractionModel(TEnvironment const& _env,
-                                                            THadronicModel& hadint)
+                                                            THadronicModel& _hadint)
       : ProposalProcessBase(_env)
-      , hadronicInteraction_(hadint) {}
+      , HadronicPhotonModel<THadronicModel>(_hadint) {}
 
   template <typename THadronicModel>
   inline void InteractionModel<THadronicModel>::buildCalculator(
@@ -107,14 +107,11 @@ namespace corsika::proposal {
 
       auto const photonEnergy = projectile.getEnergy() * v;
 
-      if (type == PROPOSAL::InteractionType::Photonuclear)
+      if (type == PROPOSAL::InteractionType::Photonuclear) {
         CORSIKA_LOG_INFO(
             "HE photo-hadronic interaction! Energy = "
             "{} GeV, v = {}, v * E = {}",
             projectile.getEnergy() / 1_GeV, v, photonEnergy);
-
-      if (type == PROPOSAL::InteractionType::Photonuclear &&
-          photonEnergy > heHadronicModelThresholdLab_) {
         auto const photonDirection =
             projectileP4.getSpaceLikeComponents()
                 .normalized(); // photon collinear with projectile
@@ -143,45 +140,54 @@ namespace corsika::proposal {
     return ProcessReturn::Ok;
   }
 
+  template <typename THadronicModel>
+  inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(THadronicModel& _hadint)
+      : heHadronicInteraction_(_hadint){};
+
   template <typename THadronicModel>
   template <typename TStackView>
-  inline ProcessReturn InteractionModel<THadronicModel>::doHadronicPhotonInteraction(
+  inline ProcessReturn HadronicPhotonModel<THadronicModel>::doHadronicPhotonInteraction(
       TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
       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
-    // target at rest
-    FourMomentum const targetP4(get_mass(targetId),
-                                MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
-    auto hadronicPhoton = photonStack.addParticle(
-        std::make_tuple(hadPhotonCode, photonP4.getTimeLikeComponent(),
-                        photonP4.getSpaceLikeComponents().normalized(), pDummy, tDummy));
-    hadronicPhoton.setNode(view.getProjectile().getNode());
-    // 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,
-                     photonP4.getTimeLikeComponent() / 1_GeV);
-    hadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
-                                       photonP4, targetP4);
-    for (const auto& pSec : photon_secondaries) {
-      auto const p3lab = pSec.getMomentum();
-      Code const pid = pSec.getPID();
-      HEPEnergyType const secEkin =
-          calculate_kinetic_energy(p3lab.getNorm(), get_mass(pid));
-      view.addSecondary(std::make_tuple(pid, secEkin, p3lab.normalized()));
+    if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLab_) {
+      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
+      // target at rest
+      FourMomentum const targetP4(get_mass(targetId),
+                                  MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
+      auto hadronicPhoton = photonStack.addParticle(std::make_tuple(
+          hadPhotonCode, photonP4.getTimeLikeComponent(),
+          photonP4.getSpaceLikeComponents().normalized(), pDummy, tDummy));
+      hadronicPhoton.setNode(view.getProjectile().getNode());
+      // 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,
+                       photonP4.getTimeLikeComponent() / 1_GeV);
+      heHadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
+                                           photonP4, targetP4);
+      for (const auto& pSec : photon_secondaries) {
+        auto const p3lab = pSec.getMomentum();
+        Code const pid = pSec.getPID();
+        HEPEnergyType const secEkin =
+            calculate_kinetic_energy(p3lab.getNorm(), get_mass(pid));
+        view.addSecondary(std::make_tuple(pid, secEkin, p3lab.normalized()));
+      }
+      CORSIKA_LOG_INFO("number of particles produced: {}", view.getEntries());
+    } else {
+      CORSIKA_LOG_INFO(
+          "LE photo-hadronic interaction! production of secondaries not implemented..");
     }
-    CORSIKA_LOG_INFO("number of particles produced: {}", view.getEntries());
     return ProcessReturn::Ok;
   }
 

corsika/modules/proposal/InteractionModel.hpp

@@ -33,7 +33,26 @@ namespace corsika::proposal {
   //!
 
   template <class THadronicModel>
-  class InteractionModel : public ProposalProcessBase {
+  class HadronicPhotonModel {
+  public:
+    HadronicPhotonModel(THadronicModel&);
+    //!
+    //! Calculate produce the hadronic secondaries in a hadronic photon interaction and
+    //! store them on the particle stack.
+    //!
+    template <typename TSecondaryView>
+    ProcessReturn doHadronicPhotonInteraction(TSecondaryView&, CoordinateSystemPtr const&,
+                                              FourMomentum const&, Code const&);
+
+  private:
+    THadronicModel& heHadronicInteraction_;
+    static HEPEnergyType constexpr heHadronicModelThresholdLab_ =
+        80. * 1e9 * electronvolt;
+  };
+
+  template <class THadronicModel>
+  class InteractionModel : public ProposalProcessBase,
+                           public HadronicPhotonModel<THadronicModel> {
 
     enum { eSECONDARIES, eINTERACTION };
     using calculator_t = std::tuple<std::unique_ptr<PROPOSAL::SecondariesCalculator>,
@@ -64,25 +83,12 @@ 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 doHadronicPhotonInteraction(TSecondaryView&, CoordinateSystemPtr const&,
-                                              FourMomentum const&, Code const&);
-
     //!
     //! Calculates and returns the cross section.
     //!
     template <typename TParticle>
     CrossSectionType getCrossSection(TParticle const& p, Code const projectileId,
                                      FourMomentum const& projectileP4);
-
-  private:
-    THadronicModel& hadronicInteraction_;
-    static HEPEnergyType constexpr heHadronicModelThresholdLab_ =
-        80. * 1e9 * electronvolt;
   };
 
 } // namespace corsika::proposal

tests/modules/testProposal.cpp

@@ -70,7 +70,20 @@ TEST_CASE("ProposalInterface", "modules") {
     CHECK(emModel.getCrossSection(particle, Code::Proton, P4) == 0_mb);
   }
 
-  SECTION("InteractionInterface - hadronic photon interaction") {
+  SECTION("InteractionInterface - LE hadronic photon interaction") {
+    auto& stack = *stackPtr;
+    // auto particle = stack.first();
+    FourMomentum P4(10_GeV, {cs, {10_GeV, 0_eV, 0_eV}});
+    // finish successfully
+    CHECK(emModel.doHadronicPhotonInteraction(view, cs, P4, Code::Oxygen) ==
+          ProcessReturn::Ok);
+    // no LE interactions
+    CHECK(stack.getEntries() == 1);
+    CORSIKA_LOG_INFO("Number of particles produced in hadronic photon interaction: {}",
+                     stack.getEntries()-1);
+  }
+
+  SECTION("InteractionInterface - HE hadronic photon interaction") {
     auto& stack = *stackPtr;
     // auto particle = stack.first();
     FourMomentum P4(100_TeV, {cs, {100_TeV, 0_eV, 0_eV}});
@@ -79,6 +92,6 @@ TEST_CASE("ProposalInterface", "modules") {
           ProcessReturn::Ok);
     CHECK(stack.getEntries() > 1);
     CORSIKA_LOG_INFO("Number of particles produced in hadronic photon interaction: {}",
-                     stack.getEntries());
+                     stack.getEntries()-1);
   }
 }
\ No newline at end of file