diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 81ba9825eb4243efa06be491db552546aceee7fe..1977467c2f91ab5513c2d491e2f468ec4ff67330 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -8,7 +8,6 @@
#include <PROPOSAL/PROPOSAL.h>
#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/NuclearComposition.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/proposal/ContinuousProcess.h>
#include <corsika/process/proposal/Interaction.h>
@@ -18,60 +17,38 @@
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
-using SetupParticle = corsika::setup::Stack::ParticleType;
-using SetupTrack = corsika::setup::Trajectory;
-
namespace corsika::process::proposal {
- using namespace corsika::setup;
- using namespace corsika::environment;
+
using namespace corsika::units::si;
- bool ContinuousProcess::CanInteract(particles::Code pcode) const noexcept {
- if (std::find(tracked_particles.begin(), tracked_particles.end(), pcode) !=
- tracked_particles.end())
- return true;
- return false;
+ void ContinuousProcess::BuildCalculator(particles::Code code,
+ environment::NuclearComposition const& comp) {
+ auto c = cross[code](media.at(&comp), cut);
+ auto disp = PROPOSAL::make_displacement(c, true);
+ auto scatter = PROPOSAL::make_scattering("highland", particle[code], media.at(&comp));
+ calc[std::make_pair(&comp, code)] =
+ std::make_tuple(std::move(disp), std::move(scatter));
}
template <>
- ContinuousProcess::ContinuousProcess(SetupEnvironment const& _env,
- CORSIKA_ParticleCut& _cut)
- : cut(_cut)
- , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
- auto all_compositions = std::vector<const NuclearComposition*>();
- _env.GetUniverse()->walk([&](auto& vtn) {
- if (vtn.HasModelProperties())
- all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
- });
- for (auto& ncarg : all_compositions) {
- auto comp_vec = std::vector<PROPOSAL::Components::Component>();
- auto frac_iter = ncarg->GetFractions().cbegin();
- for (auto& pcode : ncarg->GetComponents()) {
- comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
- *frac_iter);
- ++frac_iter;
- }
- media[ncarg] = PROPOSAL::Medium(
- "Modified Air", PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
- PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
- PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
- }
- }
+ ContinuousProcess::ContinuousProcess(setup::SetupEnvironment const& _env,
+ particle_cut::ParticleCut& _cut)
+ : ProposalProcessBase(_env, _cut) {}
template <>
- HEPEnergyType ContinuousProcess::TotalEnergyLoss(SetupParticle const& vP,
+ HEPEnergyType ContinuousProcess::TotalEnergyLoss(setup::Stack::ParticleType const& vP,
GrammageType const& vDX) {
- auto calc = GetCalculator(vP);
- return vP.GetEnergy() - get<DISPLACEMENT>(calc->second)
- ->UpperLimitTrackIntegral(vP.GetEnergy() / 1_MeV,
- vDX / 1_g * 1_cm * 1_cm) *
+ auto c = GetCalculator(vP, calc);
+ return vP.GetEnergy() - get<DISPLACEMENT>(c->second)->UpperLimitTrackIntegral(
+ vP.GetEnergy() / 1_MeV, vDX / 1_g * 1_cm * 1_cm) *
1_MeV;
}
template <>
- void ContinuousProcess::Scatter(SetupParticle& vP, HEPEnergyType const& loss,
+ void ContinuousProcess::Scatter(setup::Stack::ParticleType& vP,
+ HEPEnergyType const& loss,
GrammageType const& grammage) {
- auto calc = GetCalculator(vP);
+ auto c = GetCalculator(vP, calc);
auto d = vP.GetDirection().GetComponents();
auto direction = PROPOSAL::Vector3D(d.GetX().magnitude(), d.GetY().magnitude(),
d.GetZ().magnitude());
@@ -79,10 +56,10 @@ namespace corsika::process::proposal {
std::uniform_real_distribution<double> distr(0., 1.);
auto rnd = array<double, 4>();
for (auto& it : rnd) it = distr(fRNG);
- auto [mean_dir, final_dir] =
- get<SCATTERING>(calc->second)
- ->Scatter(grammage / 1_g * square(1_cm), vP.GetEnergy() / 1_MeV, E_f / 1_MeV,
- direction, rnd);
+ auto [mean_dir, final_dir] = get<SCATTERING>(c->second)->Scatter(
+ grammage / 1_g * square(1_cm), vP.GetEnergy() / 1_MeV, E_f / 1_MeV, direction,
+ rnd);
+ (void)mean_dir;
auto vec = corsika::geometry::QuantityVector(
final_dir.GetX() * E_f, final_dir.GetY() * E_f, final_dir.GetZ() * E_f);
vP.SetMomentum(corsika::stack::MomentumVector(
@@ -91,8 +68,8 @@ namespace corsika::process::proposal {
}
template <>
- EProcessReturn ContinuousProcess::DoContinuous(SetupParticle& vP,
- SetupTrack const& vT) {
+ EProcessReturn ContinuousProcess::DoContinuous(setup::Stack::ParticleType& vP,
+ setup::Trajectory const& vT) {
if (!CanInteract(vP.GetPID())) return process::EProcessReturn::eOk;
auto dX = vP.GetNode()->GetModelProperties().IntegratedGrammage(vT, vT.GetLength());
auto energy_loss = TotalEnergyLoss(vP, dX);
@@ -104,15 +81,15 @@ namespace corsika::process::proposal {
}
template <>
- units::si::LengthType ContinuousProcess::MaxStepLength(SetupParticle const& vP,
- SetupTrack const& vT) {
- auto energy_lim = 0.9 * vP.GetEnergy() / 1_MeV;
- auto calc = GetCalculator(vP);
- auto grammage = get<DISPLACEMENT>(calc->second)
- ->SolveTrackIntegral(vP.GetEnergy() / 1_MeV, energy_lim) *
+ units::si::LengthType ContinuousProcess::MaxStepLength(
+ setup::Stack::ParticleType const& vP, setup::Trajectory const& vT) {
+ auto energy_lim = 0.9 * vP.GetEnergy();
+ if (cut.GetECut() > energy_lim) energy_lim = cut.GetECut();
+ auto c = GetCalculator(vP, calc);
+ auto grammage = get<DISPLACEMENT>(c->second)->SolveTrackIntegral(
+ vP.GetEnergy() / 1_MeV, energy_lim / 1_MeV) *
1_g / square(1_cm);
- return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage) *
- 1.0001;
+ return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage) * 1.0001;
}
} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index cbb5a337c74c727fc9fc1574f4f666a6f3e4ed4a..c50ac8d541bb4e3073dec4476e8a27cbd2abeb3a 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -13,102 +13,70 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/proposal/ProposalProcessBase.h>
#include <corsika/random/RNGManager.h>
#include <corsika/random/UniformRealDistribution.h>
#include <unordered_map>
-using std::unordered_map;
-
-using namespace corsika::environment;
-using namespace corsika::units::si;
-
-using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
-
namespace corsika::process::proposal {
- class ContinuousProcess
- : public corsika::process::ContinuousProcess<ContinuousProcess> {
- CORSIKA_ParticleCut& cut;
- corsika::random::RNG& fRNG;
- static constexpr std::array<particles::Code, 7> tracked_particles{
- particles::Code::Gamma, particles::Code::Electron, particles::Code::Positron,
- particles::Code::MuMinus, particles::Code::MuPlus, particles::Code::TauPlus,
- particles::Code::TauMinus,
- };
- unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
-
- bool CanInteract(particles::Code pcode) const noexcept;
+ using namespace corsika::units::si;
- using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
- using calc_t =
- tuple<unique_ptr<PROPOSAL::Displacement>, unique_ptr<PROPOSAL::Scattering>>;
-
- struct disp_hash {
- size_t operator()(const calc_key_t& p) const {
- return std::hash<const NuclearComposition*>{}(p.first) ^
- std::hash<particles::Code>{}(p.second);
- }
- };
+ //!
+ //! Electro-magnetic and gamma continous losses produced by proposal. It makes
+ //! use of interpolation tables which are runtime intensive calculation, but can be
+ //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable.
+ //!
+ class ContinuousProcess : public process::ContinuousProcess<ContinuousProcess>,
+ ProposalProcessBase {
enum { DISPLACEMENT, SCATTERING };
- unordered_map<calc_key_t, calc_t, disp_hash> calc;
-
- template <typename Particle>
- auto BuildCalculator(particles::Code code, Particle p_def,
- NuclearComposition const& comp) {
- auto cross = GetStdCrossSections(
- p_def, media.at(&comp),
- make_shared<const PROPOSAL::EnergyCutSettings>(cut.GetECut() / 1_MeV, 1, false),
- true);
- auto [insert_it, success] =
- calc.insert({std::make_pair(&comp, code),
- std::make_tuple(PROPOSAL::make_displacement(cross, true),
- PROPOSAL::make_scattering("highland", p_def,
- media.at(&comp)))});
- return insert_it;
- }
+ using calc_t = std::tuple<std::unique_ptr<PROPOSAL::Displacement>,
+ std::unique_ptr<PROPOSAL::Scattering>>;
- auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
- if (corsika_code == particles::Code::Gamma)
- return BuildCalculator(particles::Code::Gamma, PROPOSAL::GammaDef(), comp);
- if (corsika_code == particles::Code::Electron)
- return BuildCalculator(particles::Code::Electron, PROPOSAL::EMinusDef(), comp);
- if (corsika_code == particles::Code::Positron)
- return BuildCalculator(particles::Code::Positron, PROPOSAL::EPlusDef(), comp);
- if (corsika_code == particles::Code::MuMinus)
- return BuildCalculator(particles::Code::MuMinus, PROPOSAL::MuMinusDef(), comp);
- if (corsika_code == particles::Code::MuPlus)
- return BuildCalculator(particles::Code::MuPlus, PROPOSAL::MuPlusDef(), comp);
- if (corsika_code == particles::Code::TauMinus)
- return BuildCalculator(particles::Code::TauMinus, PROPOSAL::TauMinusDef(), comp);
- if (corsika_code == particles::Code::TauPlus)
- return BuildCalculator(particles::Code::TauPlus, PROPOSAL::TauPlusDef(), comp);
- throw std::runtime_error("PROPOSAL could not find corresponding builder");
- }
+ std::unordered_map<calc_key_t, calc_t, hash>
+ calc; //!< Stores the displacement and scattering calculators.
- template <typename Particle>
- auto GetCalculator(Particle& vP) {
- auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
- auto calc_it = calc.find(std::make_pair(&comp, vP.GetPID()));
- if (calc_it != calc.end()) return calc_it;
- return BuildCalculator(vP.GetPID(), comp);
- }
+ //!
+ //! Build the displacement and scattering calculators and add it to calc.
+ //!
+ void BuildCalculator(particles::Code, environment::NuclearComposition const&) final;
public:
+ //! Produces the continuous loss calculator for leptons based on nuclear
+ //! compositions and stochastic description limited by the particle cut.
+ //!
template <typename TEnvironment>
- ContinuousProcess(TEnvironment const&, CORSIKA_ParticleCut&);
+ ContinuousProcess(TEnvironment const&, particle_cut::ParticleCut&);
+ //!
+ //! Calculates the energy of the particle which has passed the given
+ //! grammage.
+ //!
template <typename Particle>
corsika::units::si::HEPEnergyType TotalEnergyLoss(
Particle const&, corsika::units::si::GrammageType const&);
+ //!
+ //! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into
+ //! account. Displacment of the track due to multiple scattering is not possible
+ //! because of the constant referernce. The final direction will be updated anyway.
+ //!
template <typename Particle>
void Scatter(Particle&, corsika::units::si::HEPEnergyType const&,
corsika::units::si::GrammageType const&);
+ //!
+ //! Produces the loss and deflection after given distance for the particle.
+ //! If the particle if below the given energy threshold where it will be
+ //! considered stochastically, it will be absorbed.
+ //!
template <typename Particle, typename Track>
EProcessReturn DoContinuous(Particle&, Track const&);
+ //!
+ //! Calculates maximal step length of process.
+ //!
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle const&, Track const&);
};
diff --git a/Processes/Proposal/Interaction.cc b/Processes/Proposal/Interaction.cc
index 71ddfcfee72e230b2a1343338e03d8a9b2884480..219712d65c6f950be8d616de50e4270ae9255e03 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -8,6 +8,7 @@
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/proposal/Interaction.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
@@ -19,55 +20,34 @@
#include <random>
#include <tuple>
-using Component_PROPOSAL = PROPOSAL::Components::Component;
-
namespace corsika::process::proposal {
- using namespace corsika::setup;
using namespace corsika::environment;
using namespace corsika::units::si;
- bool Interaction::CanInteract(particles::Code pcode) const noexcept {
- if (std::find(tracked_particles.begin(), tracked_particles.end(), pcode) !=
- tracked_particles.end())
- return true;
- return false;
- }
-
template <>
- Interaction::Interaction(SetupEnvironment const& _env, CORSIKA_ParticleCut& _cut)
- : cut(_cut)
- , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
- auto all_compositions = std::vector<const NuclearComposition*>();
- _env.GetUniverse()->walk([&](auto& vtn) {
- if (vtn.HasModelProperties())
- all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
- });
- for (auto& ncarg : all_compositions) {
- auto comp_vec = std::vector<Component_PROPOSAL>();
- auto frac_iter = ncarg->GetFractions().cbegin();
- for (auto& pcode : ncarg->GetComponents()) {
- comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
- *frac_iter);
- ++frac_iter;
- }
- media[ncarg] = PROPOSAL::Medium(
- "Modified Air", PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
- PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
- PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
- }
+ Interaction::Interaction(setup::SetupEnvironment const& _env,
+ particle_cut::ParticleCut& _cut)
+ : ProposalProcessBase(_env, _cut) {}
+
+ void Interaction::BuildCalculator(particles::Code code,
+ environment::NuclearComposition const& comp) {
+ auto c = cross[code](media.at(&comp), cut);
+ auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
+ calc[std::make_pair(&comp, code)] = std::make_tuple(
+ PROPOSAL::make_secondaries(inter_types, particle[code], media.at(&comp)),
+ PROPOSAL::make_interaction(c, true));
}
template <>
corsika::process::EProcessReturn Interaction::DoInteraction(
setup::StackView::StackIterator& vP) {
if (CanInteract(vP.GetPID())) {
- auto calc = GetCalculator(vP); // [CrossSections]
+ auto c = GetCalculator(vP, calc); // [CrossSections]
std::uniform_real_distribution<double> distr(0., 1.);
- auto [type, comp_ptr, v] =
- get<INTERACTION>(calc->second)
- ->TypeInteraction(vP.GetEnergy() / 1_MeV, distr(fRNG));
- auto rnd =
- vector<double>(get<SECONDARIES>(calc->second)->RequiredRandomNumbers(type));
+ auto rates = get<INTERACTION>(c->second)->Rates(vP.GetEnergy() / 1_MeV);
+ auto [type, comp_ptr, v] = get<INTERACTION>(c->second)->SampleLoss(
+ vP.GetEnergy() / 1_MeV, rates, distr(fRNG));
+ auto rnd = vector<double>(get<SECONDARIES>(c->second)->RequiredRandomNumbers(type));
for (auto& it : rnd) it = distr(fRNG);
auto point = PROPOSAL::Vector3D(vP.GetPosition().GetX() / 1_cm,
vP.GetPosition().GetY() / 1_cm,
@@ -77,8 +57,8 @@ namespace corsika::process::proposal {
d.GetZ().magnitude());
auto loss = make_tuple(static_cast<int>(type), point, direction,
v * vP.GetEnergy() / 1_MeV, 0.);
- auto sec = get<SECONDARIES>(calc->second)
- ->CalculateSecondaries(vP.GetEnergy() / 1_MeV, loss, *comp_ptr, rnd);
+ auto sec = get<SECONDARIES>(c->second)->CalculateSecondaries(vP.GetEnergy() / 1_MeV,
+ loss, *comp_ptr, rnd);
for (auto& s : sec) {
auto E = get<PROPOSAL::Loss::ENERGY>(s) * 1_MeV;
auto vec = corsika::geometry::QuantityVector(
@@ -101,8 +81,8 @@ namespace corsika::process::proposal {
corsika::units::si::GrammageType Interaction::GetInteractionLength(
setup::Stack::StackIterator const& vP) {
if (CanInteract(vP.GetPID())) {
- auto calc = GetCalculator(vP);
- return get<INTERACTION>(calc->second)->MeanFreePath(vP.GetEnergy() / 1_MeV) * 1_g /
+ auto c = GetCalculator(vP, calc);
+ return get<INTERACTION>(c->second)->MeanFreePath(vP.GetEnergy() / 1_MeV) * 1_g /
(1_cm * 1_cm);
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
diff --git a/Processes/Proposal/Interaction.h b/Processes/Proposal/Interaction.h
index 6e1b2d57941c204a1a5a1b210c0e62f54a7295c8..990cf11156cf107acd91ff8c4f1c57871f6b8943 100644
--- a/Processes/Proposal/Interaction.h
+++ b/Processes/Proposal/Interaction.h
@@ -14,89 +14,29 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/process/proposal/ProposalProcessBase.h>
#include <corsika/random/RNGManager.h>
#include <corsika/random/UniformRealDistribution.h>
#include <array>
-using namespace corsika::environment;
-using namespace corsika::units::si;
-
-using CORSIKA_ParticleCut = corsika::process::particle_cut::ParticleCut;
-using std::make_pair;
-using std::make_tuple;
-
namespace corsika::process::proposal {
- class Interaction : public corsika::process::InteractionProcess<Interaction> {
- CORSIKA_ParticleCut& cut;
- corsika::random::RNG& fRNG;
- static constexpr std::array<particles::Code, 7> tracked_particles{
- particles::Code::Gamma, particles::Code::Electron, particles::Code::Positron,
- particles::Code::MuMinus, particles::Code::MuPlus, particles::Code::TauPlus,
- particles::Code::TauMinus,
- };
- std::unordered_map<const NuclearComposition*, PROPOSAL::Medium> media;
+ using namespace corsika::units::si;
- bool CanInteract(particles::Code pcode) const noexcept;
+ class Interaction : public InteractionProcess<Interaction>, ProposalProcessBase {
using calculator_t = tuple<unique_ptr<PROPOSAL::SecondariesCalculator>,
unique_ptr<PROPOSAL::Interaction>>;
- using calc_key_t = std::pair<const NuclearComposition*, particles::Code>;
-
- struct interaction_hash {
- size_t operator()(const calc_key_t& p) const {
- return std::hash<const NuclearComposition*>{}(p.first) ^
- std::hash<particles::Code>{}(p.second);
- }
- };
- std::unordered_map<calc_key_t, calculator_t, interaction_hash> calculators;
+ std::unordered_map<calc_key_t, calculator_t, hash> calc;
- template <typename Particle>
- auto BuildCalculator(particles::Code code, Particle p_def,
- NuclearComposition const& comp) {
- auto cross = GetStdCrossSections(
- p_def, media.at(&comp),
- make_shared<const PROPOSAL::EnergyCutSettings>(cut.GetECut() / 1_MeV, 1, false),
- true);
- auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(cross);
- auto [insert_it, success] = calculators.insert(
- {make_pair(&comp, code),
- make_tuple(PROPOSAL::make_secondaries(inter_types, p_def, media.at(&comp)),
- PROPOSAL::make_interaction(cross, true))});
- return insert_it;
- }
-
- auto BuildCalculator(particles::Code corsika_code, NuclearComposition const& comp) {
- if (corsika_code == particles::Code::Gamma)
- return BuildCalculator(particles::Code::Gamma, PROPOSAL::GammaDef(), comp);
- if (corsika_code == particles::Code::Electron)
- return BuildCalculator(particles::Code::Electron, PROPOSAL::EMinusDef(), comp);
- if (corsika_code == particles::Code::Positron)
- return BuildCalculator(particles::Code::Positron, PROPOSAL::EPlusDef(), comp);
- if (corsika_code == particles::Code::MuMinus)
- return BuildCalculator(particles::Code::MuMinus, PROPOSAL::MuMinusDef(), comp);
- if (corsika_code == particles::Code::MuPlus)
- return BuildCalculator(particles::Code::MuPlus, PROPOSAL::MuPlusDef(), comp);
- if (corsika_code == particles::Code::TauMinus)
- return BuildCalculator(particles::Code::TauMinus, PROPOSAL::TauMinusDef(), comp);
- if (corsika_code == particles::Code::TauPlus)
- return BuildCalculator(particles::Code::TauPlus, PROPOSAL::TauPlusDef(), comp);
- throw std::runtime_error("PROPOSAL could not find corresponding builder");
- }
+ void BuildCalculator(particles::Code, environment::NuclearComposition const&) final;
enum { SECONDARIES, INTERACTION };
- template <typename Particle>
- auto GetCalculator(Particle& vP) {
- auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
- auto calc_it = calculators.find(make_pair(&comp, vP.GetPID()));
- if (calc_it != calculators.end()) return calc_it;
- return BuildCalculator(vP.GetPID(), comp);
- }
public:
template <typename TEnvironment>
- Interaction(TEnvironment const& env, CORSIKA_ParticleCut& cut);
+ Interaction(TEnvironment const& env, particle_cut::ParticleCut& cut);
template <typename Particle>
corsika::process::EProcessReturn DoInteraction(Particle&);
diff --git a/Processes/Proposal/PROPOSAL b/Processes/Proposal/PROPOSAL
index 17bc2e3a6fb777c51993dca12831c45a8eb40ff8..92b9b9ca20826725cb805d135a1a5d5b29a7e06a 160000
--- a/Processes/Proposal/PROPOSAL
+++ b/Processes/Proposal/PROPOSAL
@@ -1 +1 @@
-Subproject commit 17bc2e3a6fb777c51993dca12831c45a8eb40ff8
+Subproject commit 92b9b9ca20826725cb805d135a1a5d5b29a7e06a
diff --git a/Processes/Proposal/ProposalProcessBase.cc b/Processes/Proposal/ProposalProcessBase.cc
new file mode 100644
index 0000000000000000000000000000000000000000..335cf41860160dbbd73df8e79dcd94535a306eef
--- /dev/null
+++ b/Processes/Proposal/ProposalProcessBase.cc
@@ -0,0 +1,60 @@
+/*
+ * (c) Copyright 2020 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.
+ */
+
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/proposal/ProposalProcessBase.h>
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+#include <limits>
+#include <memory>
+#include <random>
+#include <tuple>
+
+namespace corsika::process::proposal {
+ bool ProposalProcessBase::CanInteract(particles::Code pcode) const {
+ if (std::find(begin(tracked), end(tracked), pcode) != end(tracked)) return true;
+ return false;
+ }
+
+ ProposalProcessBase::ProposalProcessBase(setup::SetupEnvironment const& _env,
+ particle_cut::ParticleCut& _cut)
+ : cut(_cut)
+ , fRNG(corsika::random::RNGManager::GetInstance().GetRandomStream("proposal")) {
+ auto all_compositions = std::vector<const environment::NuclearComposition*>();
+ _env.GetUniverse()->walk([&](auto& vtn) {
+ if (vtn.HasModelProperties())
+ all_compositions.push_back(&vtn.GetModelProperties().GetNuclearComposition());
+ });
+ for (auto& ncarg : all_compositions) {
+ auto comp_vec = std::vector<PROPOSAL::Components::Component>();
+ auto frac_iter = ncarg->GetFractions().cbegin();
+ for (auto& pcode : ncarg->GetComponents()) {
+ comp_vec.emplace_back(GetName(pcode), GetNucleusZ(pcode), GetNucleusA(pcode),
+ *frac_iter);
+ ++frac_iter;
+ }
+ media[ncarg] = PROPOSAL::Medium(
+ "Modified Air", PROPOSAL::Air().GetI(), PROPOSAL::Air().GetC(),
+ PROPOSAL::Air().GetA(), PROPOSAL::Air().GetM(), PROPOSAL::Air().GetX0(),
+ PROPOSAL::Air().GetX1(), PROPOSAL::Air().GetD0(), 1.0, comp_vec);
+ }
+ PROPOSAL::InterpolationDef::order_of_interpolation = 2;
+ PROPOSAL::InterpolationDef::nodes_cross_section = 100;
+ PROPOSAL::InterpolationDef::nodes_propagate = 100;
+ }
+
+ size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const noexcept {
+ return std::hash<const environment::NuclearComposition*>{}(p.first) ^
+ std::hash<particles::Code>{}(p.second);
+ }
+
+} // namespace corsika::process::proposal
diff --git a/Processes/Proposal/ProposalProcessBase.h b/Processes/Proposal/ProposalProcessBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2f503e775f4584d0a1809ee0c8b3d6bc57eed5a
--- /dev/null
+++ b/Processes/Proposal/ProposalProcessBase.h
@@ -0,0 +1,122 @@
+#pragma once
+
+#include <PROPOSAL/PROPOSAL.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/particle_cut/ParticleCut.h>
+#include <corsika/random/RNGManager.h>
+#include <array>
+
+namespace corsika::process::proposal {
+
+ using namespace corsika::units::si;
+ using namespace std::placeholders;
+
+ //!
+ //! Particles which can be handled by proposal. That means they can be
+ //! propagated and decayed if they decays.
+ //!
+ static constexpr std::array<particles::Code, 7> tracked{
+ particles::Code::Gamma, particles::Code::Electron, particles::Code::Positron,
+ particles::Code::MuMinus, particles::Code::MuPlus, particles::Code::TauPlus,
+ particles::Code::TauMinus,
+ };
+
+ //!
+ //! Internal map from particle codes to particle properties required for
+ //! crosssections, decay and scattering algorithms. In the future the
+ //! particles may be created by reading out the Corsica constants.
+ //!
+ static std::map<particles::Code, PROPOSAL::ParticleDef> particle = {
+ {particles::Code::Gamma, PROPOSAL::GammaDef()},
+ {particles::Code::Electron, PROPOSAL::EMinusDef()},
+ {particles::Code::Positron, PROPOSAL::EPlusDef()},
+ {particles::Code::MuMinus, PROPOSAL::MuMinusDef()},
+ {particles::Code::MuPlus, PROPOSAL::MuPlusDef()},
+ {particles::Code::TauMinus, PROPOSAL::TauMinusDef()},
+ {particles::Code::TauPlus, PROPOSAL::TauPlusDef()}};
+
+ //!
+ //! Crosssection factories for different particle types.
+ //!
+ template <typename T>
+ static auto cross_builder = [](PROPOSAL::Medium& m, particle_cut::ParticleCut& cut) {
+ auto p_cut = std::make_shared<const PROPOSAL::EnergyCutSettings>(
+ cut.GetECut() / 1_MeV, 1, true);
+ return PROPOSAL::DefaultCrossSections<T>::template Get<std::false_type>(T(), m, p_cut,
+ true);
+ };
+
+ //!
+ //! PROPOSAL default crosssections are maped to corresponding corsika particle
+ //! code.
+ //!
+ static std::map<particles::Code, std::function<PROPOSAL::crosssection_list_t<
+ PROPOSAL::ParticleDef, PROPOSAL::Medium>(
+ PROPOSAL::Medium&, particle_cut::ParticleCut&)>>
+ cross = {{particles::Code::Gamma, cross_builder<PROPOSAL::GammaDef>},
+ {particles::Code::Electron, cross_builder<PROPOSAL::EMinusDef>},
+ {particles::Code::Positron, cross_builder<PROPOSAL::EPlusDef>},
+ {particles::Code::MuMinus, cross_builder<PROPOSAL::MuMinusDef>},
+ {particles::Code::MuPlus, cross_builder<PROPOSAL::MuPlusDef>},
+ {particles::Code::TauMinus, cross_builder<PROPOSAL::TauMinusDef>},
+ {particles::Code::TauPlus, cross_builder<PROPOSAL::TauPlusDef>}};
+
+ //!
+ //! PROPOSAL base process which handels mapping of particle codes to
+ //! stored interpolation tables.
+ //!
+ class ProposalProcessBase {
+ protected:
+ particle_cut::ParticleCut& cut; //!< Stochastic losses smaller than the given cut
+ //!< will be handeled continuously.
+ corsika::random::RNG& fRNG; //!< random number generator used by proposal
+
+ std::unordered_map<const environment::NuclearComposition*, PROPOSAL::Medium>
+ media; //!< maps nuclear composition from univers to media to produce
+ //!< crosssections, which requires further ionization constants.
+
+ //!
+ //! Store cut and nuclear composition of the whole universe in media which are
+ //! required for creating crosssections by proposal.
+ //!
+ ProposalProcessBase(setup::SetupEnvironment const& _env,
+ particle_cut::ParticleCut& _cut);
+
+ //!
+ //! Checks if a particle can be processed by proposal
+ //!
+ bool CanInteract(particles::Code pcode) const;
+
+ using calc_key_t = std::pair<const environment::NuclearComposition*, particles::Code>;
+
+ //!
+ //! Hash to store interpolation tables related to a pair of particle and nuclear
+ //! composition.
+ //!
+ struct hash {
+ size_t operator()(const calc_key_t& p) const noexcept;
+ };
+
+ //!
+ //! Builds the calculator to the corresponding class
+ //!
+ virtual void BuildCalculator(particles::Code,
+ environment::NuclearComposition const&) = 0;
+
+ //!
+ //! Searches the particle dependet calculator dependent of actuall medium composition
+ //! and particle type. If no calculator is found, the corresponding new calculator is
+ //! built and then returned.
+ //!
+ template <typename Particle, typename Calculators>
+ auto GetCalculator(Particle& vP, Calculators& calc) {
+ auto& comp = vP.GetNode()->GetModelProperties().GetNuclearComposition();
+ auto calc_it = calc.find(std::make_pair(&comp, vP.GetPID()));
+ if (calc_it != calc.end()) return calc_it;
+ BuildCalculator(vP.GetPID(), comp);
+ return GetCalculator(vP, calc);
+ }
+ };
+} // namespace corsika::process::proposal