fix interface

corsika/detail/modules/proposal/HadronicPhotonModel.inl

@@ -47,7 +47,7 @@ namespace corsika::proposal {
 
       // call inner hadronic event generator
       CORSIKA_LOG_TRACE("calling HadronicInteraction...");
-      CORSIKA_LOG_INFO("{} + {} interactions. Ekinlab = {}", hadPhotonCode, targetId,
+      CORSIKA_LOG_INFO("{} + {} interactions. Ekinlab = {} GeV", hadPhotonCode, targetId,
                        photonP4.getTimeLikeComponent() / 1_GeV);
       heHadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
                                            photonP4, targetP4);

corsika/detail/modules/proposal/InteractionModel.inl

@@ -77,7 +77,8 @@ namespace corsika::proposal {
         CORSIKA_LOG_WARN(
             "PROPOSAL: No particle interaction possible. "
             "Put initial particle back on stack.");
-        view.addSecondary(std::make_tuple(projectileId, projectile.getEnergy(),
+        view.addSecondary(std::make_tuple(projectileId,
+                                          projectile.getEnergy() - get_mass(projectileId),
                                           projectile.getDirection()));
         return ProcessReturn::Ok;
       }
@@ -104,34 +105,31 @@ namespace corsika::proposal {
       if (type != PROPOSAL::InteractionType::Ioniz)
         target = PROPOSAL::Component::GetComponentForHash(target_hash);
 
-      auto const photonEnergy = projectile.getEnergy() * v;
-
-      if (type == PROPOSAL::InteractionType::Photonuclear) {
-        auto const photonDirection =
-            projectileP4.getSpaceLikeComponents()
-                .normalized(); // photon collinear with projectile
-        FourMomentum const photonP4(photonEnergy, photonEnergy * photonDirection);
-        Code const targetId =
-            get_nucleus_code(target.GetAtomicNum(), target.GetNucCharge());
-        CORSIKA_LOG_INFO(
-            "photo-hadronic interaction ({} + {})! Energy = "
-            "{} GeV, v = {}, Photon energy (v*E) = {} GeV",
-            projectileId, targetId, projectile.getEnergy() / 1_GeV, v,
-            photonEnergy / 1_GeV);
-        this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
-        if (projectileId != Code::Photon)
-          // add lepton, apply energy loss to kinetic energy
-          view.addSecondary(std::make_tuple(
-              projectileId, (1 - v) * (projectile.getEnergy() - get_mass(projectileId)),
-              photonDirection));
-      } 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 =
+          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()});
+
+        if (static_cast<PROPOSAL::ParticleType>(s.type) ==
+            PROPOSAL::ParticleType::Hadron) {
+          FourMomentum const photonP4(E, E * dir);
+          // for PROPOSAL media target.GetAtomicNum() returns the atomic number in units
+          // of atomic mass, so Nitrogen is 14.0067. When using media in CORSIKA the same
+          // field is filled with the number of nucleons (ie. 14 for Nitrogen) To be sure
+          // we explicitly cast to int
+          auto const A = int(target.GetAtomicNum());
+          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);
+          this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
+        } else {
           auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
           view.addSecondary(std::make_tuple(sec_code, E - get_mass(sec_code), dir));
         }

tests/modules/testProposal.cpp

@@ -65,9 +65,14 @@ TEST_CASE("ProposalInterface", "modules") {
     auto& stack = *stackPtr;
     auto particle = stack.first();
     FourMomentum P4(
-        100_GeV,
+        100_MeV,
         {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);
+
+    FourMomentum eleP4(
+        100_MeV,
+        {cs, {sqrt(static_pow<2>(100_MeV) - static_pow<2>(Electron::mass)), 0_eV, 0_eV}});
+    CHECK(emModel.getCrossSection(particle, Code::Electron, eleP4) > 0_mb);
   }
 
   SECTION("InteractionInterface - LE hadronic photon interaction") {