From d05e531b68c9e6ffe048459dbb763149340196da Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sun, 9 Oct 2022 22:55:04 +0200
Subject: [PATCH] added Sophia as LE model to hadronicPhoton in proposal
---
.../modules/proposal/HadronicPhotonModel.inl | 106 ++++++++++--------
.../modules/proposal/InteractionModel.inl | 23 ++--
.../modules/sophia/InteractionModel.inl | 2 +-
corsika/modules/PROPOSAL.hpp | 13 ++-
.../modules/proposal/HadronicPhotonModel.hpp | 16 +--
corsika/modules/proposal/InteractionModel.hpp | 10 +-
corsika/modules/sophia/InteractionModel.hpp | 8 +-
tests/modules/testProposal.cpp | 25 +++--
8 files changed, 119 insertions(+), 84 deletions(-)
diff --git a/corsika/detail/modules/proposal/HadronicPhotonModel.inl b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
index 0eeb4afc6..1275fe35a 100644
--- a/corsika/detail/modules/proposal/HadronicPhotonModel.inl
+++ b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
@@ -15,10 +15,12 @@
namespace corsika::proposal {
- template <typename THadronicModel>
- inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(
- THadronicModel& _hadint, HEPEnergyType const& _heenthresholdNN)
- : heHadronicInteraction_(_hadint)
+ template <typename THadronicLEModel, typename THadronicHEModel>
+ inline HadronicPhotonModel<THadronicLEModel, THadronicHEModel>::HadronicPhotonModel(
+ THadronicLEModel& _hadintLE, THadronicHEModel& _hadintHE,
+ HEPEnergyType const& _heenthresholdNN)
+ : leHadronicInteraction_(_hadintLE)
+ , heHadronicInteraction_(_hadintHE)
, heHadronicModelThresholdLabNN_(_heenthresholdNN) {
// check validity of threshold assuming photon-nucleon
// sqrtS per target nucleon
@@ -27,8 +29,8 @@ namespace corsika::proposal {
if (!heHadronicInteraction_.isValid(Code::Rho0, Code::Proton, sqrtS)) {
CORSIKA_LOGGER_CRITICAL(
logger_,
- "Invalid energy threshold for hadron interaction model. theshold_lab= {} GeV, "
- "theshold_com={} GeV",
+ "Invalid energy threshold for hadron interaction model. threshold_lab= {} GeV, "
+ "threshold_com={} GeV",
_heenthresholdNN / 1_GeV, sqrtS / 1_GeV);
throw std::runtime_error("Configuration error!");
}
@@ -37,65 +39,77 @@ namespace corsika::proposal {
_heenthresholdNN / 1_GeV);
}
- template <typename THadronicModel>
+ template <typename THadronicLEModel, typename THadronicHEModel>
template <typename TStackView>
- inline ProcessReturn HadronicPhotonModel<THadronicModel>::doHadronicPhotonInteraction(
+ inline ProcessReturn
+ HadronicPhotonModel<THadronicLEModel, THadronicHEModel>::doHadronicPhotonInteraction(
TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
Code const& targetId) {
- if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLabNN_) {
- CORSIKA_LOGGER_TRACE(
- logger_, "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);
+ // 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;
+ // target at rest
+ FourMomentum const targetP4(get_mass(targetId),
+ MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
+ auto hadronicPhoton = photonStack.addParticle(
+ std::make_tuple(Code::Photon, 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_LOGGER_TRACE(logger_, "{} + {} interaction. Ekinlab = {} GeV",
- hadPhotonCode, targetId,
- photonP4.getTimeLikeComponent() / 1_GeV);
- // check if had. model can handle configuration
- auto const sqrtSNN = (photonP4 + targetP4 / get_nucleus_A(targetId)).getNorm();
+ // call inner hadronic event generator
+ CORSIKA_LOGGER_TRACE(logger_, "{} + {} interaction. Ekinlab = {} GeV", Code::Photon,
+ targetId, photonP4.getTimeLikeComponent() / 1_GeV);
+ // check if had. model can handle configuration
+ auto const sqrtSNN = (photonP4 + targetP4 / get_nucleus_A(targetId)).getNorm();
+ CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={} GeV", sqrtSNN / 1_GeV);
+ if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLabNN_) {
+ CORSIKA_LOGGER_TRACE(logger_, "HE photo-hadronic interaction!");
// when Sibyll is used for hadronic interactions Argon cannot be used as target
// nucleus. Since PROPOSAL has a non-zero cross section for Argon
// targets we have to check here if the model can handle Argon (see Issue #498)
- if (!heHadronicInteraction_.isValid(hadPhotonCode, targetId, sqrtSNN)) {
+ if (!heHadronicInteraction_.isValid(Code::Rho0, targetId, sqrtSNN)) {
CORSIKA_LOGGER_WARN(
logger_,
"HE interaction model cannot handle configuration in photo-hadronic "
"interaction! projectile={}, target={} (A={}, Z={}), sqrt(S) per "
"nuc.={:8.2f} "
"GeV. Skipping secondary production!",
- hadPhotonCode, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
+ Code::Rho0, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
sqrtSNN / 1_GeV);
return ProcessReturn::Ok;
}
- heHadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
+ heHadronicInteraction_.doInteraction(photon_secondaries, Code::Rho0, 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()));
- }
} else {
- CORSIKA_LOGGER_TRACE(
- logger_,
- "LE photo-hadronic interaction! Production of secondaries not implemented..");
+ CORSIKA_LOGGER_TRACE(logger_,
+ "LE photo-hadronic interaction! implemented via SOPHIA "
+ "assuming a single nucleon as target");
+ // as LE model we only have SOPHIA which only handles photon-nucleon interactions
+ if (!leHadronicInteraction_.isValid(Code::Photon, Code::Proton, sqrtSNN)) {
+ CORSIKA_LOGGER_WARN(
+ logger_,
+ "LE interaction model cannot handle configuration in photo-hadronic "
+ "interaction! projectile={}, target={} (A={}, Z={}), sqrt(S) per "
+ "nuc.={:8.2f} "
+ "GeV. Skipping secondary production!",
+ Code::Photon, targetId, get_nucleus_A(targetId), get_nucleus_Z(targetId),
+ sqrtSNN / 1_GeV);
+ return ProcessReturn::Ok;
+ }
+ leHadronicInteraction_.doInteraction(photon_secondaries, Code::Photon, Code::Proton,
+ photonP4, targetP4 / get_nucleus_A(targetId));
+ }
+ 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()));
}
return ProcessReturn::Ok;
}
diff --git a/corsika/detail/modules/proposal/InteractionModel.inl b/corsika/detail/modules/proposal/InteractionModel.inl
index 8a0390e90..d8198bab5 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -18,16 +18,17 @@
namespace corsika::proposal {
- template <typename THadronicModel>
+ template <typename THadronicLEModel, typename THadronicHEModel>
template <typename TEnvironment>
- inline InteractionModel<THadronicModel>::InteractionModel(
- TEnvironment const& _env, THadronicModel& _hadint,
+ inline InteractionModel<THadronicLEModel, THadronicHEModel>::InteractionModel(
+ TEnvironment const& _env, THadronicLEModel& _hadintLE, THadronicHEModel& _hadintHE,
HEPEnergyType const& _enthreshold)
: ProposalProcessBase(_env)
- , HadronicPhotonModel<THadronicModel>(_hadint, _enthreshold) {}
+ , HadronicPhotonModel<THadronicLEModel, THadronicHEModel>(_hadintLE, _hadintHE,
+ _enthreshold) {}
- template <typename THadronicModel>
- inline void InteractionModel<THadronicModel>::buildCalculator(
+ template <typename THadronicLEModel, typename THadronicHEModel>
+ inline void InteractionModel<THadronicLEModel, THadronicHEModel>::buildCalculator(
Code code, NuclearComposition const& comp) {
// search crosssection builder for given particle
auto p_cross = cross.find(code);
@@ -52,9 +53,10 @@ namespace corsika::proposal {
PROPOSAL::make_interaction(c, true, true));
}
- template <typename THadronicModel>
+ template <typename THadronicLEModel, typename THadronicHEModel>
template <typename TStackView>
- inline ProcessReturn InteractionModel<THadronicModel>::doInteraction(
+ inline ProcessReturn
+ InteractionModel<THadronicLEModel, THadronicHEModel>::doInteraction(
TStackView& view, Code const projectileId, FourMomentum const& projectileP4) {
auto const projectile = view.getProjectile();
@@ -138,9 +140,10 @@ namespace corsika::proposal {
return ProcessReturn::Ok;
}
- template <typename THadronicModel>
+ template <typename THadronicLEModel, typename THadronicHEModel>
template <typename TParticle>
- inline CrossSectionType InteractionModel<THadronicModel>::getCrossSection(
+ inline CrossSectionType
+ InteractionModel<THadronicLEModel, THadronicHEModel>::getCrossSection(
TParticle const& projectile, Code const projectileId,
FourMomentum const& projectileP4) {
diff --git a/corsika/detail/modules/sophia/InteractionModel.inl b/corsika/detail/modules/sophia/InteractionModel.inl
index 89560a392..fee1146cb 100644
--- a/corsika/detail/modules/sophia/InteractionModel.inl
+++ b/corsika/detail/modules/sophia/InteractionModel.inl
@@ -80,7 +80,7 @@ namespace corsika::sophia {
initial_(nucleonSophiaCode);
double Enucleon = nucleonP4.getTimeLikeComponent() / 1_GeV;
double Ephoton = photonP4.getTimeLikeComponent() / 1_GeV;
- double theta = 0.0; // set head on collision
+ double theta = 0.0; // set nucleon at rest in collision
int Imode;
CORSIKA_LOGGER_DEBUG(logger_,
"calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}",
diff --git a/corsika/modules/PROPOSAL.hpp b/corsika/modules/PROPOSAL.hpp
index 13e18c3de..81066911f 100644
--- a/corsika/modules/PROPOSAL.hpp
+++ b/corsika/modules/PROPOSAL.hpp
@@ -13,12 +13,15 @@
namespace corsika::proposal {
- template <typename THadronicModel>
- class Interaction : public InteractionModel<THadronicModel>,
- public InteractionProcess<Interaction<THadronicModel>> {
+ template <typename THadronicLEModel, typename THadronicHEModel>
+ class Interaction
+ : public InteractionModel<THadronicLEModel, THadronicHEModel>,
+ public InteractionProcess<Interaction<THadronicLEModel, THadronicHEModel>> {
public:
template <typename TEnvironment>
- Interaction(TEnvironment const& env, THadronicModel& model, HEPEnergyType const& thr)
- : InteractionModel<THadronicModel>(env, model, thr) {}
+ Interaction(TEnvironment const& env, THadronicLEModel& modelLE,
+ THadronicHEModel& modelHE, HEPEnergyType const& thr)
+ : InteractionModel<THadronicLEModel, THadronicHEModel>(env, modelLE, modelHE,
+ thr) {}
};
} // namespace corsika::proposal
diff --git a/corsika/modules/proposal/HadronicPhotonModel.hpp b/corsika/modules/proposal/HadronicPhotonModel.hpp
index 6fe6b88bb..602d3caf5 100644
--- a/corsika/modules/proposal/HadronicPhotonModel.hpp
+++ b/corsika/modules/proposal/HadronicPhotonModel.hpp
@@ -15,17 +15,18 @@
namespace corsika::proposal {
//! Implements the production of secondary hadrons for the hadronic interaction of real
- //! and virtual photons. At high energies an external model
- //! is needed that implements the doInteraction(TSecondaries& view, Code const
- //! projectile, Code const target,FourMomentum const& projectileP4, FourMomentum const&
- //! targetP4) routine. Low energy interactions are currently not implemented. The
+ //! and virtual photons for PROPOSAL. The model distinguishes between resonance
+ //! production (LE) and continuum (HE). HE production replaces the photon with a rho0.
+ //! External models are needed that implement the hadronic particle production via the
+ //! doInteraction(TSecondaries& view, Code const projectile, Code const target,
+ //! FourMomentum const& projectileP4, FourMomentum const& targetP4) routine. The
//! threshold between LE and HE interactions is defined in lab energy.
//! @tparam THadronicModel
- template <class THadronicModel>
+ template <class THadronicLEModel, class THadronicHEModel>
class HadronicPhotonModel {
public:
- HadronicPhotonModel(THadronicModel&, HEPEnergyType const&);
+ HadronicPhotonModel(THadronicLEModel&, THadronicHEModel&, HEPEnergyType const&);
//!
//! Calculate produce the hadronic secondaries in a hadronic photon interaction and
//! store them on the particle stack.
@@ -36,7 +37,8 @@ namespace corsika::proposal {
private:
inline static auto logger_{get_logger("corsika_proposal_HadronicPhotonModel")};
- THadronicModel& heHadronicInteraction_;
+ THadronicLEModel& leHadronicInteraction_;
+ THadronicHEModel& heHadronicInteraction_;
//! threshold for high energy hadronic interaction model. Lab. energy per nucleon
HEPEnergyType const heHadronicModelThresholdLabNN_;
};
diff --git a/corsika/modules/proposal/InteractionModel.hpp b/corsika/modules/proposal/InteractionModel.hpp
index f25815322..ae9d3b7e7 100644
--- a/corsika/modules/proposal/InteractionModel.hpp
+++ b/corsika/modules/proposal/InteractionModel.hpp
@@ -33,9 +33,10 @@ namespace corsika::proposal {
//! @tparam THadronicModel
//!
- template <class THadronicModel>
- class InteractionModel : public ProposalProcessBase,
- public HadronicPhotonModel<THadronicModel> {
+ template <class THadronicLEModel, class THadronicHEModel>
+ class InteractionModel
+ : public ProposalProcessBase,
+ public HadronicPhotonModel<THadronicLEModel, THadronicHEModel> {
enum { eSECONDARIES, eINTERACTION };
using calculator_t = std::tuple<std::unique_ptr<PROPOSAL::SecondariesCalculator>,
@@ -57,7 +58,8 @@ namespace corsika::proposal {
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
- InteractionModel(TEnvironment const& env, THadronicModel&, HEPEnergyType const&);
+ InteractionModel(TEnvironment const& env, THadronicLEModel&, THadronicHEModel&,
+ HEPEnergyType const&);
//!
//! Calculate the rates for the different targets and interactions. Sample a
diff --git a/corsika/modules/sophia/InteractionModel.hpp b/corsika/modules/sophia/InteractionModel.hpp
index e2f6d7636..f6e6d8fd7 100644
--- a/corsika/modules/sophia/InteractionModel.hpp
+++ b/corsika/modules/sophia/InteractionModel.hpp
@@ -62,9 +62,11 @@ namespace corsika::sophia {
* @return inelastic cross section
* elastic cross section
*/
- CrossSectionType getCrossSection(Code const projectile, Code const target,
- FourMomentum const& projectileP4,
- FourMomentum const& targetP4) const {
+ CrossSectionType getCrossSection(
+ [[maybe_unused]] Code const projectile, [[maybe_unused]] Code const target,
+ [[maybe_unused]] FourMomentum const& projectileP4,
+ [[maybe_unused]] FourMomentum const& targetP4) const {
+ CORSIKA_LOGGER_ERROR(logger_, "cross section not implemented in SOPHIA!");
return CrossSectionType::zero();
}
/**
diff --git a/tests/modules/testProposal.cpp b/tests/modules/testProposal.cpp
index eb186c358..299c913fe 100644
--- a/tests/modules/testProposal.cpp
+++ b/tests/modules/testProposal.cpp
@@ -12,6 +12,7 @@
#include <SetupTestStack.hpp>
#include <catch2/catch.hpp>
#include <tuple>
+#include "corsika/framework/core/PhysicalUnits.hpp"
using namespace corsika;
using namespace corsika::proposal;
@@ -26,7 +27,8 @@ using DummyEnvironment = Environment<DummyEnvironmentInterface>;
class DummyHadronicModel {
public:
- DummyHadronicModel(){};
+ DummyHadronicModel(HEPEnergyType thr)
+ : threshold_(thr){};
template <typename TSecondaryView>
void doInteraction(TSecondaryView& view, Code const, Code const,
@@ -42,8 +44,11 @@ public:
}
}
bool constexpr isValid(Code const, Code const, HEPEnergyType const sqrsNN) const {
- return (sqrsNN >= 10_GeV);
+ return (sqrsNN >= threshold_);
};
+
+private:
+ HEPEnergyType threshold_;
};
TEST_CASE("ProposalInterface", "modules") {
@@ -61,14 +66,18 @@ TEST_CASE("ProposalInterface", "modules") {
RNGManager<>::getInstance().registerRandomStream("proposal");
SECTION("InteractionInterface - hadronic photon model threshold") {
- DummyHadronicModel hadModel;
- HEPEnergyType heThresholdLab1 = 10_GeV;
- CHECK_THROWS(corsika::proposal::InteractionModel(*env, hadModel, heThresholdLab1));
+ DummyHadronicModel hadModelLE(100_MeV);
+ DummyHadronicModel hadModelHE(10_GeV);
+ HEPEnergyType heThresholdLab1 = 12_GeV;
+ CHECK_THROWS(corsika::proposal::InteractionModel(*env, hadModelLE, hadModelHE,
+ heThresholdLab1));
}
- DummyHadronicModel hadModel;
+ DummyHadronicModel hadModelLE(100_MeV);
+ DummyHadronicModel hadModelHE(10_GeV);
HEPEnergyType heThresholdLab = 80_GeV;
- corsika::proposal::InteractionModel emModel(*env, hadModel, heThresholdLab);
+ corsika::proposal::InteractionModel emModel(*env, hadModelLE, hadModelHE,
+ heThresholdLab);
SECTION("InteractionInterface - cross section") {
auto& stack = *stackPtr;
@@ -92,7 +101,7 @@ TEST_CASE("ProposalInterface", "modules") {
CHECK(emModel.doHadronicPhotonInteraction(view, cs, P4, Code::Oxygen) ==
ProcessReturn::Ok);
// no LE interactions
- CHECK(stack.getEntries() == 1);
+ CHECK(stack.getEntries() == 6);
CORSIKA_LOG_INFO("Number of particles produced in hadronic photon interaction: {}",
stack.getEntries() - 1);
}
--
GitLab