From 1db03312c4705534c8759e323d441b3ae0c781cb Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 12 Apr 2022 12:34:26 +0100
Subject: [PATCH] split off hadronic photon model
---
.../modules/proposal/HadronicPhotonModel.inl | 65 +++++++++++++++++++
.../modules/proposal/InteractionModel.inl | 63 ++----------------
.../modules/proposal/HadronicPhotonModel.hpp | 36 ++++++++++
corsika/modules/proposal/InteractionModel.hpp | 19 +-----
tests/modules/testProposal.cpp | 4 +-
5 files changed, 110 insertions(+), 77 deletions(-)
create mode 100644 corsika/detail/modules/proposal/HadronicPhotonModel.inl
create mode 100644 corsika/modules/proposal/HadronicPhotonModel.hpp
diff --git a/corsika/detail/modules/proposal/HadronicPhotonModel.inl b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
new file mode 100644
index 000000000..2b63b7a52
--- /dev/null
+++ b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
@@ -0,0 +1,65 @@
+#pragma once
+
+//#include <corsika/media/IMediumModel.hpp>
+//#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+//#include <limits>
+//#include <memory>
+//#include <random>
+#include <tuple>
+
+namespace corsika::proposal {
+
+ template <typename THadronicModel>
+ inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(THadronicModel& _hadint)
+ : heHadronicInteraction_(_hadint){};
+
+ 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_) {
+ 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..");
+ }
+ return ProcessReturn::Ok;
+ }
+} // namespace corsika::proposal
\ No newline at end of file
diff --git a/corsika/detail/modules/proposal/InteractionModel.inl b/corsika/detail/modules/proposal/InteractionModel.inl
index 2838c1676..44db596d8 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -10,7 +10,6 @@
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <limits>
#include <memory>
@@ -108,17 +107,18 @@ namespace corsika::proposal {
auto const photonEnergy = projectile.getEnergy() * v;
if (type == PROPOSAL::InteractionType::Photonuclear) {
- CORSIKA_LOG_INFO(
- "HE photo-hadronic interaction! Energy = "
- "{} GeV, v = {}, v * E = {}",
- projectile.getEnergy() / 1_GeV, v, photonEnergy);
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());
- doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
+ 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(
@@ -140,57 +140,6 @@ 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 HadronicPhotonModel<THadronicModel>::doHadronicPhotonInteraction(
- TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
- Code const& targetId) {
- 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..");
- }
- return ProcessReturn::Ok;
- }
-
template <typename THadronicModel>
template <typename TParticle>
inline CrossSectionType InteractionModel<THadronicModel>::getCrossSection(
diff --git a/corsika/modules/proposal/HadronicPhotonModel.hpp b/corsika/modules/proposal/HadronicPhotonModel.hpp
new file mode 100644
index 000000000..08c94d510
--- /dev/null
+++ b/corsika/modules/proposal/HadronicPhotonModel.hpp
@@ -0,0 +1,36 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/process/ProcessReturn.hpp>
+
+namespace corsika::proposal {
+
+ template <class THadronicModel>
+ 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;
+ };
+} // namespace corsika::proposal
+
+#include <corsika/detail/modules/proposal/HadronicPhotonModel.inl>
\ No newline at end of file
diff --git a/corsika/modules/proposal/InteractionModel.hpp b/corsika/modules/proposal/InteractionModel.hpp
index 9c3381db2..ec5739225 100644
--- a/corsika/modules/proposal/InteractionModel.hpp
+++ b/corsika/modules/proposal/InteractionModel.hpp
@@ -17,6 +17,7 @@
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/random/UniformRealDistribution.hpp>
#include <corsika/modules/proposal/ProposalProcessBase.hpp>
+#include <corsika/modules/proposal/HadronicPhotonModel.hpp>
namespace corsika::proposal {
@@ -32,24 +33,6 @@ namespace corsika::proposal {
//! @tparam THadronicModel
//!
- template <class THadronicModel>
- 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> {
diff --git a/tests/modules/testProposal.cpp b/tests/modules/testProposal.cpp
index 8718f1c52..650cfaab7 100644
--- a/tests/modules/testProposal.cpp
+++ b/tests/modules/testProposal.cpp
@@ -80,7 +80,7 @@ TEST_CASE("ProposalInterface", "modules") {
// no LE interactions
CHECK(stack.getEntries() == 1);
CORSIKA_LOG_INFO("Number of particles produced in hadronic photon interaction: {}",
- stack.getEntries()-1);
+ stack.getEntries() - 1);
}
SECTION("InteractionInterface - HE hadronic photon interaction") {
@@ -92,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()-1);
+ stack.getEntries() - 1);
}
}
\ No newline at end of file
--
GitLab