IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 1db03312 authored by Felix Riehn's avatar Felix Riehn Committed by Maximilian Reininghaus
Browse files

split off hadronic photon model

parent ffe2508e
No related branches found
No related tags found
No related merge requests found
#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
......@@ -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(
......
/*
* (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
......@@ -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> {
......
......@@ -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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment