diff --git a/Processes/Proposal/CMakeLists.txt b/Processes/Proposal/CMakeLists.txt
index 970abe29dc0ed296479e2ba20df2f134f0c0d007..8d2a19788215b616ca2d6dd28199212ce8fc2918 100644
--- a/Processes/Proposal/CMakeLists.txt
+++ b/Processes/Proposal/CMakeLists.txt
@@ -15,6 +15,11 @@ SET_TARGET_PROPERTIES ( ProcessProposal PROPERTIES VERSION ${PROJECT_VERSION}
TARGET_LINK_LIBRARIES (
ProcessProposal
+ CORSIKAparticles
+ CORSIKAutilities
+ CORSIKAunits
+ CORSIKAthirdparty
+ CORSIKAgeometry
CORSIKAenvironment
${PROPOSAL_LIBRARY}
)
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index a19b6bf5ff82a9ecb96689d47156e111bbaddee8..7faaa9dbb8b75f28acd77bc69748d6d79b37ad6e 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -2,6 +2,7 @@
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/process/proposal/Interaction.h>
+#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
@@ -14,12 +15,11 @@
#include <tuple>
using Component_PROPOSAL = PROPOSAL::Components::Component;
-namespace corsika::process::proposal {
+namespace corsika::process::proposal {
using namespace corsika::setup;
using namespace corsika::environment;
using namespace corsika::units::si;
- using SetupParticle = corsika::setup::Stack::StackIterator;
template <>
std::unordered_map<particles::Code, PROPOSAL::ParticleDef>
@@ -65,43 +65,42 @@ namespace corsika::process::proposal {
template <>
template <>
corsika::process::EProcessReturn Interaction<SetupEnvironment>::DoInteraction(
- SetupParticle& vP) {
- auto calc = GetCalculator(vP); // [CrossSections]
- std::uniform_real_distribution<double> distr(0., 1.);
- auto [type, comp_ptr, v] = std::get<INTERACTION>(calc->second)
- ->TypeInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
- auto rnd = std::vector<double>();
- for (size_t i = 0;
- i < std::get<SECONDARIES>(calc->second).RequiredRandomNumbers(type); ++i)
- rnd.push_back(distr(fRNG));
- double primary_energy = vP.GetEnergy() / 1_GeV;
- auto point =
- PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm, vP.GetPosition().GetY() / 1_cm,
- vP.GetPosition().GetZ() / 1_cm);
- auto p = vP.GetMomentum().GetComponents();
- auto direction = PROPOSAL::Vector3D(p[0] / 1_GeV, p[1] / 1_GeV, p[2] / 1_GeV);
- auto loss =
- make_tuple(static_cast<int>(type), point, direction, v * primary_energy, 0.);
- auto sec = std::get<SECONDARIES>(calc->second)
- .CalculateSecondaries(primary_energy, loss, *comp_ptr, rnd);
- for (auto& s : sec) {
- auto energy = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
- auto vec = corsika::geometry::QuantityVector(
- std::get<PROPOSAL::Loss::DIRECTION>(s).GetX() * energy,
- std::get<PROPOSAL::Loss::DIRECTION>(s).GetY() * energy,
- std::get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * energy);
- auto momentum = corsika::stack::MomentumVector(
- corsika::geometry::RootCoordinateSystem::GetInstance()
- .GetRootCoordinateSystem(),
- vec);
- particles::Code sec_code = corsika::particles::ConvertFromPDG(
- static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
- corsika::units::si::HEPEnergyType sec_energy = energy;
- stack::MomentumVector sec_momentum = momentum;
- geometry::Point sec_position = vP.GetPosition();
- corsika::units::si::TimeType sec_time = vP.GetTime();
- vP.AddSecondary(
- make_tuple(sec_code, sec_energy, sec_momentum, sec_position, sec_time));
+ setup::StackView::StackIterator& vP) {
+ if (CanInteract(vP.GetPID())) {
+ auto calc = GetCalculator(vP); // [CrossSections]
+ std::uniform_real_distribution<double> distr(0., 1.);
+ auto [type, comp_ptr, v] =
+ std::get<INTERACTION>(calc->second)
+ ->TypeInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
+ auto rnd = std::vector<double>();
+ for (size_t i = 0;
+ i < std::get<SECONDARIES>(calc->second).RequiredRandomNumbers(type); ++i)
+ rnd.push_back(distr(fRNG));
+ double primary_energy = vP.GetEnergy() / 1_GeV;
+ auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
+ vP.GetPosition().GetY() / 1_cm,
+ vP.GetPosition().GetZ() / 1_cm);
+ auto p = vP.GetMomentum().GetComponents();
+ auto direction = PROPOSAL::Vector3D(p[0] / 1_GeV, p[1] / 1_GeV, p[2] / 1_GeV);
+ auto loss =
+ make_tuple(static_cast<int>(type), point, direction, v * primary_energy, 0.);
+ auto sec = std::get<SECONDARIES>(calc->second)
+ .CalculateSecondaries(primary_energy, loss, *comp_ptr, rnd);
+ for (auto& s : sec) {
+ auto energy = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
+ auto vec = corsika::geometry::QuantityVector(
+ std::get<PROPOSAL::Loss::DIRECTION>(s).GetX() * energy,
+ std::get<PROPOSAL::Loss::DIRECTION>(s).GetY() * energy,
+ std::get<PROPOSAL::Loss::DIRECTION>(s).GetZ() * energy);
+ auto momentum = corsika::stack::MomentumVector(
+ corsika::geometry::RootCoordinateSystem::GetInstance()
+ .GetRootCoordinateSystem(),
+ vec);
+ particles::Code sec_code = corsika::particles::ConvertFromPDG(
+ static_cast<particles::PDGCode>(get<PROPOSAL::Loss::TYPE>(s)));
+ vP.AddSecondary(
+ make_tuple(sec_code, energy, momentum, vP.GetPosition(), vP.GetTime()));
+ }
}
return process::EProcessReturn::eOk;
}
@@ -109,13 +108,16 @@ namespace corsika::process::proposal {
template <>
template <>
corsika::units::si::GrammageType Interaction<SetupEnvironment>::GetInteractionLength(
- SetupParticle const& vP) {
- auto calc = GetCalculator(vP); // [CrossSections]
- std::uniform_real_distribution<double> distr(0., 1.);
- auto energy = get<INTERACTION>(calc->second)
- ->EnergyInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
- return get<DISPLACEMENT>(calc->second)
- ->SolveTrackIntegral(vP.GetEnergy() / 1_GeV, energy) *
- 1_g / 1_cm / 1_cm;
+ setup::Stack::StackIterator& vP) {
+ if (CanInteract(vP.GetPID())) {
+ auto calc = GetCalculator(vP); // [CrossSections]
+ std::uniform_real_distribution<double> distr(0., 1.);
+ auto energy = get<INTERACTION>(calc->second)
+ ->EnergyInteraction(vP.GetEnergy() / 1_GeV, distr(fRNG));
+ return get<DISPLACEMENT>(calc->second)
+ ->SolveTrackIntegral(vP.GetEnergy() / 1_GeV, energy) *
+ 1_g / 1_cm / 1_cm;
+ }
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index 8692dd4e68f0860f361aad77a762873a31adc122..503827f356892d636303a62c33f0ed31d5cb882d 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -38,12 +38,10 @@ namespace corsika::process::proposal {
static std::unordered_map<particles::Code, PROPOSAL::ParticleDef> particles;
std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
- enum { SECONDARIES, INTERACTION, DISPLACEMENT };
-
corsika::random::RNG& fRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ corsika::random::RNGManager::GetInstance().GetRandomStream("p_rndm");
- auto IsTracked(particles::Code pcode) const noexcept {
+ bool CanInteract(particles::Code pcode) const noexcept {
auto search = particles.find(pcode);
if (search != particles.end()) return true;
return false;
@@ -54,6 +52,7 @@ namespace corsika::process::proposal {
unique_ptr<PROPOSAL::Displacement>>;
std::unordered_map<const NuclearComposition*, calculator_t> calculators;
+ enum { SECONDARIES, INTERACTION, DISPLACEMENT };
template <typename Particle>
auto GetCalculator(Particle& vP) {
auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
@@ -81,6 +80,5 @@ namespace corsika::process::proposal {
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
};
-
} // namespace corsika::process::proposal
#endif