IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 17924bb8 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

Merge branch '517-sophia-for-low-energy-photo-hadronic-interaction' into 'master' 

Resolve "SOPHIA for low energy photo-hadronic interaction"

Closes #517

See merge request !465
parents 3da34a87 07762640
No related branches found
No related tags found
1 merge request!465Resolve "SOPHIA for low energy photo-hadronic interaction"
Pipeline #9653 passed with warnings
Stage: quality
Stage: config
Stage: build_test_example
Stage: python
Hide whitespace changes
Inline Side-by-side
Showing
with 724 additions and 162 deletions
CMakeLists.txt
+ 1
− 0
View file @ 17924bb8
......@@ -253,6 +253,7 @@ set (CORSIKA_DATA_WITH_TEST ON) # we want to run the corsika-data unit test
add_subdirectory (modules/data) # this is corsika-data (submodule)
add_subdirectory (modules/pythia8)
add_subdirectory (modules/sibyll)
add_subdirectory (modules/sophia)
add_subdirectory (modules/qgsjetII)
add_subdirectory (modules/urqmd)
add_subdirectory (modules/conex)
......
corsika/detail/modules/proposal/HadronicPhotonModel.inl
+ 74
− 46
View file @ 17924bb8
......@@ -12,13 +12,16 @@
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <tuple>
#include <random>
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 +30,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 +40,90 @@ 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");
// sample nucleon from nucleus A,Z
double const fProtons = get_nucleus_Z(targetId) / double(get_nucleus_A(targetId));
double const fNeutrons = 1. - fProtons;
std::discrete_distribution<int> nucleonChannelDist{fProtons, fNeutrons};
corsika::default_prng_type& rng =
corsika::RNGManager<>::getInstance().getRandomStream("proposal");
Code const nucleonId = (nucleonChannelDist(rng) ? Code::Neutron : Code::Proton);
// target passed to SOPHIA needs to be exactly on-shell!
FourMomentum const nucleonP4(get_mass(nucleonId),
MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
CORSIKA_LOGGER_DEBUG(logger_,
"selected {} as target nucleon (f_proton, f_neutron)={},{}",
nucleonId, fProtons, fNeutrons);
if (!leHadronicInteraction_.isValid(Code::Photon, nucleonId, 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, nucleonId,
photonP4, nucleonP4);
}
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;
}
......
corsika/detail/modules/proposal/InteractionModel.inl
+ 13
− 10
View file @ 17924bb8
......@@ -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) {
......
corsika/detail/modules/sophia/InteractionModel.inl 0 → 100644
+ 139
− 0
View file @ 17924bb8
/*
* (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/geometry/Point.hpp>
#include <corsika/modules/sophia/ParticleConversion.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <corsika/modules/sophia/SophiaStack.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <sophia.hpp>
namespace corsika::sophia {
inline void InteractionModel::setVerbose(bool const flag) { sophia_listing_ = flag; }
inline InteractionModel::InteractionModel()
: sophia_listing_(false) {
// set all particles stable in SOPHIA
for (int i = 0; i < 49; ++i) so_csydec_.idb[i] = -abs(so_csydec_.idb[i]);
}
inline InteractionModel::~InteractionModel() {
CORSIKA_LOG_DEBUG("Sophia::Model n={}", count_);
}
inline bool constexpr InteractionModel::isValid(Code const projectileId,
Code const targetId,
HEPEnergyType const sqrtSnn) const {
if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
if (!(targetId == Code::Proton || targetId == Code::Neutron ||
targetId == Code::Hydrogen))
return false;
if (projectileId != Code::Photon) return false;
return true;
}
template <typename TSecondaryView>
inline void InteractionModel::doInteraction(TSecondaryView& secondaries,
Code const projectileId,
Code const targetId,
FourMomentum const& projectileP4,
FourMomentum const& targetP4) {
CORSIKA_LOGGER_DEBUG(logger_, "projectile: Id={}, E={} GeV, p3={} GeV", projectileId,
projectileP4.getTimeLikeComponent() / 1_GeV,
projectileP4.getSpaceLikeComponents().getComponents() / 1_GeV);
CORSIKA_LOGGER_DEBUG(logger_, "target: Id={}, E={} GeV, p3={} GeV", targetId,
targetP4.getTimeLikeComponent() / 1_GeV,
targetP4.getSpaceLikeComponents().getComponents() / 1_GeV);
// sqrtS per target nucleon
HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={}GeV", sqrtS / 1_GeV);
// accepts only photon-nucleon interactions
if (!isValid(projectileId, targetId, sqrtS)) {
CORSIKA_LOGGER_ERROR(logger_, "Invalid target/projectile/energy combination");
throw std::runtime_error("SOPHIA: Invalid target/projectile/energy combination");
}
COMBoost const boost(projectileP4, targetP4);
int nucleonSophiaCode = convertToSophiaRaw(targetId); // either proton or neutron
// initialize resonance spectrum
initial_(nucleonSophiaCode);
double Enucleon = targetP4.getTimeLikeComponent() / 1_GeV;
double Ephoton = projectileP4.getTimeLikeComponent() / 1_GeV;
double theta = 0.0; // set nucleon at rest in collision
int Imode = -1; // overwritten inside SOPHIA
CORSIKA_LOGGER_DEBUG(logger_,
"calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}",
nucleonSophiaCode, Enucleon, Ephoton, theta);
count_++;
// call sophia
eventgen_(nucleonSophiaCode, Enucleon, Ephoton, theta, Imode);
if (sophia_listing_) {
int arg = 3;
print_event_(arg);
}
auto const& originalCS = boost.getOriginalCS();
// SOPHIA has photon along -z and nucleon along +z (GZK calc..)
COMBoost const boostInternal(targetP4, projectileP4);
auto const& csPrime = boost.getRotatedCS();
CoordinateSystemPtr csPrimePrime =
make_rotation(csPrime, QuantityVector<length_d>{1_m, 0_m, 0_m}, M_PI);
SophiaStack ss;
MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType E_final = 0_GeV;
for (auto& psop : ss) {
// abort on particles that have decayed in SOPHIA. Should not happen!
if (psop.hasDecayed()) { // LCOV_EXCL_START
throw std::runtime_error("found particle that decayed in SOPHIA!");
} // LCOV_EXCL_STOP
auto momentumSophia = psop.getMomentum(csPrimePrime);
momentumSophia.rebase(csPrime);
auto const energySophia = psop.getEnergy();
auto const P4com = boostInternal.toCoM(FourVector{energySophia, momentumSophia});
auto const P4lab = boost.fromCoM(P4com);
SophiaCode const pidSophia = psop.getPID();
Code const pid = convertFromSophia(pidSophia);
auto momentum = P4lab.getSpaceLikeComponents();
momentum.rebase(originalCS);
HEPEnergyType const Ekin =
calculate_kinetic_energy(momentum.getNorm(), get_mass(pid));
CORSIKA_LOGGER_TRACE(logger_, "SOPHIA: pid={}, p={} GeV", pidSophia,
momentumSophia.getComponents() / 1_GeV);
CORSIKA_LOGGER_TRACE(logger_, "CORSIKA: pid={}, p={} GeV", pid,
momentum.getComponents() / 1_GeV);
auto pnew =
secondaries.addSecondary(std::make_tuple(pid, Ekin, momentum.normalized()));
P_final += pnew.getMomentum();
E_final += pnew.getEnergy();
}
CORSIKA_LOGGER_TRACE(logger_, "Efinal={} GeV,Pfinal={} GeV", E_final / 1_GeV,
P_final.getComponents() / 1_GeV);
}
} // namespace corsika::sophia
\ No newline at end of file
corsika/detail/modules/sophia/ParticleConversion.inl 0 → 100644
+ 25
− 0
View file @ 17924bb8
/*
* (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 <sophia.hpp>
namespace corsika::sophia {
inline HEPMassType getSophiaMass(Code const pCode) {
if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
auto sCode = convertToSophiaRaw(pCode);
if (sCode == 0)
throw std::runtime_error("getSophiaMass: unknown particle!");
else
return sqrt(get_sophia_mass2(sCode)) * 1_GeV;
}
} // namespace corsika::sophia
\ No newline at end of file
corsika/modules/PROPOSAL.hpp
+ 8
− 5
View file @ 17924bb8
......@@ -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
corsika/modules/Random.hpp
+ 1
− 0
View file @ 17924bb8
......@@ -17,6 +17,7 @@
link this togehter, it will fail.
*/
#include <corsika/modules/sibyll/Random.hpp>
#include <corsika/modules/sophia/Random.hpp>
#include <corsika/modules/epos/Random.hpp>
#include <corsika/modules/urqmd/Random.hpp>
#include <corsika/modules/qgsjetII/Random.hpp>
......
corsika/modules/Sophia.hpp 0 → 100644
+ 15
− 0
View file @ 17924bb8
/*
* (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/modules/sophia/ParticleConversion.hpp>
#include <corsika/modules/sophia/InteractionModel.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
corsika/modules/proposal/HadronicPhotonModel.hpp
+ 9
− 7
View file @ 17924bb8
......@@ -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_;
};
......
corsika/modules/proposal/InteractionModel.hpp
+ 6
− 4
View file @ 17924bb8
......@@ -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
......
corsika/modules/sibyll/HadronInteractionModel.hpp
+ 6
− 0
View file @ 17924bb8
......@@ -91,6 +91,12 @@ namespace corsika::sibyll {
/**
* In this function SIBYLL is called to produce one event. The
* event is copied (and boosted) into the shower lab frame.
*
* @param view is the stack object for the secondaries
* @param projectile is the Code of the projectile
* @param target is the Code of the target
* @param projectileP4: four-momentum of projectile
* @param targetP4: four-momentum of target
*/
template <typename TSecondaries>
......
corsika/modules/sophia/InteractionModel.hpp 0 → 100644
+ 107
− 0
View file @ 17924bb8
/*
* (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/core/PhysicalUnits.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <tuple>
namespace corsika::sophia {
/**
* @brief Provides the SOPHIA photon-nucleon interaction model.
*
* This is a TModel argument for InteractionProcess<TModel>.
*/
class InteractionModel {
public:
InteractionModel();
~InteractionModel();
/**
* @brief Set the Verbose flag.
*
* If flag is true, SOPHIA will printout additional secondary particle information
* lists, etc.
*
* @param flag to switch.
*/
void setVerbose(bool const flag);
/**
* @brief evaluated validity of collision system.
*
* SOPHIA only accepts nucleons as targets, that is protons (Hydrogen) or
* neutrons.
*/
bool constexpr isValid(Code const projectileId, Code const targetId,
HEPEnergyType const sqrtSnn) const;
/**
* Returns inelastic (production) cross section.
*
* This cross section must correspond to the process described in doInteraction.
* Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
*
* @param projectile is the Code of the projectile
* @param target is the Code of the target
* @param projectileP4: four-momentum of projectile
* @param targetP4: four-momentum of target
*
* @return inelastic cross section
* elastic cross section
*/
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();
}
/**
* In this function SOPHIA is called to produce one event. The
* event is copied (and boosted) into the frame of the incoming particles.
*
* @param view is the stack object for the secondaries
* @param projectile is the Code of the projectile
* @param target is the Code of the target
* @param projectileP4: four-momentum of projectile
* @param targetP4: four-momentum of target
*/
template <typename TSecondaries>
void doInteraction(TSecondaries& view, Code const projectile, Code const target,
FourMomentum const& projectileP4, FourMomentum const& targetP4);
private:
HEPEnergyType constexpr getMinEnergyCoM() const { return minEnergyCoM_; }
HEPEnergyType constexpr getMaxEnergyCoM() const { return maxEnergyCoM_; }
// hard model limits
static HEPEnergyType constexpr minEnergyCoM_ = 1.079166345 * 1e9 * electronvolt;
static HEPEnergyType constexpr maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sophia");
// data members
int count_ = 0;
bool sophia_listing_;
std::shared_ptr<spdlog::logger> logger_ =
get_logger("corsika_sophia_InteractionModel");
};
} // namespace corsika::sophia
#include <corsika/detail/modules/sophia/InteractionModel.inl>
\ No newline at end of file
corsika/modules/sophia/ParticleConversion.hpp 0 → 100644
+ 52
− 0
View file @ 17924bb8
/*
* (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/core/PhysicalUnits.hpp>
#include <sophia.hpp>
#include <string>
namespace corsika::sophia {
enum class SophiaCode : int8_t;
using SophiaCodeIntType = std::underlying_type<SophiaCode>::type;
#include <corsika/modules/sophia/Generated.inc>
SophiaCode constexpr convertToSophia(Code const pCode) {
return corsika2sophia[static_cast<CodeIntType>(pCode)];
}
Code constexpr convertFromSophia(SophiaCode const pCode) {
auto const s = static_cast<SophiaCodeIntType>(pCode);
auto const corsikaCode = sophia2corsika[s - minSophia];
if (corsikaCode == Code::Unknown) {
throw std::runtime_error(std::string("SOPHIA/CORSIKA conversion of ")
.append(std::to_string(s))
.append(" impossible"));
}
return corsikaCode;
}
int constexpr convertToSophiaRaw(Code const code) {
return static_cast<int>(convertToSophia(code));
}
bool constexpr canInteract(Code const pCode) {
return (pCode == Code::Photon ? true : false);
}
HEPMassType getSophiaMass(Code const);
} // namespace corsika::sophia
#include <corsika/detail/modules/sophia/ParticleConversion.inl>
corsika/modules/sophia/Random.hpp 0 → 100644
+ 31
− 0
View file @ 17924bb8
/*
* (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/random/RNGManager.hpp>
#include <random>
/**
* \file sophia/Random.hpp
*
* This file is an integral part of the sophia interface. It must be
* linked to the executable linked to sophia exactly once
*
*/
namespace sophia {
double rndm_interface() {
static corsika::default_prng_type& rng =
corsika::RNGManager<>::getInstance().getRandomStream("sophia");
std::uniform_real_distribution<double> dist;
return dist(rng);
}
} // namespace sophia
corsika/modules/sophia/SophiaStack.hpp 0 → 100644
+ 126
− 0
View file @ 17924bb8
/*
* (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/PhysicalUnits.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/stack/Stack.hpp>
#include <corsika/modules/sophia/ParticleConversion.hpp>
#include <sophia.hpp>
namespace corsika::sophia {
typedef corsika::Vector<hepmomentum_d> MomentumVector;
class SophiaStackData {
public:
void dump() const {}
void clear() { so_plist_.np = 0; }
unsigned int getSize() const { return so_plist_.np; }
unsigned int getCapacity() const { return 2000; }
void setId(const unsigned int i, const int v) { so_plist_.llist[i] = v; }
void setEnergy(const unsigned int i, const HEPEnergyType v) {
so_plist_.p[3][i] = v / 1_GeV;
}
void setMass(const unsigned int i, const HEPMassType v) {
so_plist_.p[4][i] = v / 1_GeV;
}
void setMomentum(const unsigned int i, const MomentumVector& v) {
auto tmp = v.getComponents();
for (int idx = 0; idx < 3; ++idx) so_plist_.p[idx][i] = tmp[idx] / 1_GeV;
}
int getId(const unsigned int i) const { return so_plist_.llist[i]; }
HEPEnergyType getEnergy(const int i) const { return so_plist_.p[3][i] * 1_GeV; }
HEPEnergyType getMass(const unsigned int i) const {
return so_plist_.p[4][i] * 1_GeV;
}
MomentumVector getMomentum(const unsigned int i,
const CoordinateSystemPtr& CS) const {
QuantityVector<hepmomentum_d> components = {so_plist_.p[0][i] * 1_GeV,
so_plist_.p[1][i] * 1_GeV,
so_plist_.p[2][i] * 1_GeV};
return MomentumVector(CS, components);
}
void copy(const unsigned int i1, const unsigned int i2) {
so_plist_.llist[i2] = so_plist_.llist[i1];
for (unsigned int i = 0; i < 5; ++i) so_plist_.p[i][i2] = so_plist_.p[i][i1];
}
void swap(const unsigned int i1, const unsigned int i2) {
std::swap(so_plist_.llist[i1], so_plist_.llist[i2]);
for (unsigned int i = 0; i < 5; ++i)
std::swap(so_plist_.p[i][i1], so_plist_.p[i][i2]);
}
void incrementSize() { so_plist_.np++; }
void decrementSize() {
if (so_plist_.np > 0) { so_plist_.np--; }
}
};
template <typename TStackIterator>
class ParticleInterface : public corsika::ParticleBase<TStackIterator> {
using corsika::ParticleBase<TStackIterator>::getStackData;
using corsika::ParticleBase<TStackIterator>::getIndex;
public:
void setParticleData(const int vID, const HEPEnergyType vE, const MomentumVector& vP,
const HEPMassType vM) {
setPID(vID);
setEnergy(vE);
setMomentum(vP);
setMass(vM);
}
void setParticleData(ParticleInterface<TStackIterator>& /*parent*/, const int vID,
const HEPEnergyType vE, const MomentumVector& vP,
const HEPMassType vM) {
setPID(vID);
setEnergy(vE);
setMomentum(vP);
setMass(vM);
}
void setEnergy(const HEPEnergyType v) { getStackData().setEnergy(getIndex(), v); }
HEPEnergyType getEnergy() const { return getStackData().getEnergy(getIndex()); }
bool hasDecayed() const { return abs(getStackData().getId(getIndex())) > 100; }
void setMass(const HEPMassType v) { getStackData().setMass(getIndex(), v); }
HEPEnergyType getMass() const { return getStackData().getMass(getIndex()); }
void setPID(const int v) { getStackData().setId(getIndex(), v); }
corsika::sophia::SophiaCode getPID() const {
return static_cast<corsika::sophia::SophiaCode>(getStackData().getId(getIndex()));
}
MomentumVector getMomentum(const CoordinateSystemPtr& CS) const {
return getStackData().getMomentum(getIndex(), CS);
}
void setMomentum(const MomentumVector& v) {
getStackData().setMomentum(getIndex(), v);
}
};
typedef corsika::Stack<SophiaStackData, ParticleInterface> SophiaStack;
} // namespace corsika::sophia
examples/corsika.cpp
+ 36
− 28
View file @ 17924bb8
......@@ -12,57 +12,58 @@
// to include it first...
#include <corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/GeomagneticModel.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/Epos.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/OnShellCheck.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/QGSJetII.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/Epos.hpp>
#include <corsika/modules/Sophia.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/UrQMD.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/QGSJetII.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <CLI/App.hpp>
#include <CLI/Formatter.hpp>
#include <CLI/Config.hpp>
#include <CLI/Formatter.hpp>
#include <iomanip>
#include <limits>
......@@ -90,6 +91,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("qgsjet");
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("epos");
RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("urqmd");
......@@ -313,18 +315,22 @@ int main(int argc, char** argv) {
// decaySibyll.printDecayConfig();
// hadronic photon interactions in resonance region
corsika::sophia::InteractionModel sophia;
HEPEnergyType const emcut = 50_GeV;
HEPEnergyType const hadcut = 50_GeV;
ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, dEdX);
// energy threshold for high energy hadronic model. Affects LE/HE switch for hadron
// interactions and the hadronic photon model in proposal
// energy threshold for high energy hadronic model. Affects LE/HE switch for
// hadron interactions and the hadronic photon model in proposal
HEPEnergyType heHadronModelThreshold = 63.1_GeV;
corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
heHadronModelThreshold);
corsika::proposal::Interaction emCascade(
env, sophia, sibyll.getHadronInteractionModel(), heHadronModelThreshold);
// NOT available for PROPOSAL due to interface trouble:
// InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
// emContinuous(env);
BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
LongitudinalWriter profile{showerAxis, 200, 10_g / square(1_cm)};
......@@ -361,7 +367,8 @@ int main(int argc, char** argv) {
observationLevel, longprof);
/* === END: SETUP PROCESS LIST === */
// create the cascade object using the default stack and tracking implementation
// create the cascade object using the default stack and tracking
// implementation
setup::Tracking tracking;
setup::Stack<EnvType> stack;
Cascade EAS(env, tracking, sequence, output, stack);
......@@ -416,7 +423,8 @@ int main(int argc, char** argv) {
Efinal / 1_GeV, dEdX.getEnergyLost() / 1_GeV,
observationLevel.getEnergyGround() / 1_GeV, (Efinal / E0 - 1) * 100);
// auto const hists = heModelCounted.getHistogram() + urqmdCounted.getHistogram();
// auto const hists = heModelCounted.getHistogram() +
// urqmdCounted.getHistogram();
auto const hists = sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
save_hist(hists.labHist(), labHist_file, true);
......
examples/em_shower.cpp
+ 18
− 16
View file @ 17924bb8
......@@ -6,41 +6,41 @@
* the license.
*/
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/Sophia.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
......@@ -68,6 +68,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("sophia");
if (seed == 0) {
std::random_device rd;
seed = rd();
......@@ -168,9 +169,10 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(2_MeV, 2_MeV, 100_GeV, 100_GeV, true, dEdX);
corsika::sibyll::Interaction sibyll{env};
corsika::sophia::InteractionModel sophia;
HEPEnergyType heThresholdNN = 60_GeV;
corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
heThresholdNN);
corsika::proposal::Interaction emCascade(
env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
......
examples/mars.cpp
+ 24
− 21
View file @ 17924bb8
......@@ -12,56 +12,57 @@
// to include it first...
#include <corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/SlidingPlanarExponential.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/QGSJetII.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/Sophia.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/UrQMD.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/QGSJetII.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <CLI/App.hpp>
#include <CLI/Formatter.hpp>
#include <CLI/Config.hpp>
#include <CLI/Formatter.hpp>
#include <iomanip>
#include <iostream>
......@@ -118,6 +119,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("qgsjet");
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("urqmd");
RNGManager<>::getInstance().registerRandomStream("proposal");
......@@ -350,12 +352,12 @@ int main(int argc, char** argv) {
// decaySibyll.printDecayConfig();
// energy threshold for high energy hadronic model. Affects LE/HE switch for hadron
// interactions and the hadronic photon model in proposal
// energy threshold for high energy hadronic model. Affects LE/HE switch for
// hadron interactions and the hadronic photon model in proposal
HEPEnergyType heHadronModelThreshold = 63.1_GeV;
corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
heHadronModelThreshold);
corsika::sophia::InteractionModel sophia;
corsika::proposal::Interaction emCascade(
env, sophia, sibyll.getHadronInteractionModel(), heHadronModelThreshold);
// NOT possible right now, due to interface difference for PROPOSAL:
// InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
......@@ -398,7 +400,8 @@ int main(int argc, char** argv) {
cut, trackWriter, observationLevel, profile);
/* === END: SETUP PROCESS LIST === */
// create the cascade object using the default stack and tracking implementation
// create the cascade object using the default stack and tracking
// implementation
setup::Tracking tracking;
setup::Stack<EnvType> stack;
Cascade EAS(env, tracking, sequence, output, stack);
......
examples/radio_em_shower.cpp
+ 6
− 2
View file @ 17924bb8
......@@ -40,6 +40,7 @@
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/Sophia.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/radio/RadioProcess.hpp>
......@@ -78,6 +79,7 @@ using namespace std;
void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("sibyll");
if (seed == 0) {
std::random_device rd;
......@@ -239,10 +241,12 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
corsika::sophia::InteractionModel sophia;
corsika::sibyll::Interaction sibyll{env};
HEPEnergyType heThresholdNN = 80_GeV;
corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
heThresholdNN);
corsika::proposal::Interaction emCascade(
env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
......
examples/vertical_EAS.cpp
+ 27
− 23
View file @ 17924bb8
......@@ -14,47 +14,48 @@
// to include it first...
#include <corsika/framework/process/InteractionCounter.hpp>
/* clang-format on */
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/GeomagneticModel.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/Epos.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/Sophia.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/UrQMD.hpp>
#include <corsika/modules/Epos.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
......@@ -85,6 +86,7 @@ using Particle = setup::Stack<EnvType>::particle_type;
void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("sibyll");
// RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("urqmd");
RNGManager<>::getInstance().registerRandomStream("proposal");
......@@ -217,8 +219,8 @@ int main(int argc, char** argv) {
// construct the continuous energy loss model
BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
// construct a particle cut - cuts are set to values close to reality, put higher
// values for faster runs
// construct a particle cut - cuts are set to values close to reality, put
// higher values for faster runs
ParticleCut<SubWriter<decltype(dEdX)>> cut{2_MeV, 2_MeV, 2_GeV, 300_MeV, true, dEdX};
// setup longitudinal profile
......@@ -238,8 +240,10 @@ int main(int argc, char** argv) {
HEPEnergyType heThresholdNN = 60_GeV;
// PROPOSAL is disabled for this example
// corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
// heThresholdNN); corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
// corsika::sophia::InteractionModel sophia;
// corsika::proposal::Interaction emCascade(env, sophia,
// sibyll.getHadronInteractionModel(), heThresholdNN);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
// emContinuous(env, dEdX);
corsika::pythia8::Decay decayPythia;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum