From 370d68ad27f7d53c69e353fbc0e65ab1e314cf3c Mon Sep 17 00:00:00 2001
From: Maximilian Sackel <maximilian.sackel@tu-dortmund.de>
Date: Wed, 30 Sep 2020 10:23:35 +0000
Subject: [PATCH] add type checking before processing particle with proposal
---
Documentation/Examples/vertical_EAS.cc | 14 +++++++-------
Processes/Proposal/ContinuousProcess.cc | 10 ++++++++--
Processes/Proposal/Interaction.cc | 5 ++++-
Processes/Proposal/ProposalProcessBase.cc | 2 +-
4 files changed, 20 insertions(+), 11 deletions(-)
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index d78f1f624..26b259c84 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -188,7 +188,7 @@ int main(int argc, char** argv) {
PROPOSAL::InterpolationDef::path_to_tables = "~/.local/share/PROPOSAL/tables/";
PROPOSAL::InterpolationDef::path_to_tables_readonly = "~/.local/share/PROPOSAL/tables/";
- process::particle_cut::ParticleCut cut{60_GeV, true, true};
+ process::particle_cut::ParticleCut cut{60_GeV, false, true};
process::proposal::Interaction proposal(env, cut);
process::proposal::ContinuousProcess em_continuous(env, cut);
process::interaction_counter::InteractionCounter proposalCounted(proposal);
@@ -205,9 +205,9 @@ int main(int argc, char** argv) {
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
- process::conex_source_cut::CONEXSourceCut conexSource(
- center, showerAxis, t, injectionHeight, E0,
- particles::GetPDG(particles::Code::Proton));
+ /* process::conex_source_cut::CONEXSourceCut conexSource( */
+ /* center, showerAxis, t, injectionHeight, E0, */
+ /* particles::GetPDG(particles::Code::Proton)); */
// assemble all processes into an ordered process list
@@ -218,9 +218,9 @@ int main(int argc, char** argv) {
//auto sequence = switchProcess << reset_particle_mass << decaySequence << conexSource
// << longprof << eLoss << cut << observationLevel;
- auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << cut << em_continuous
+ auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted << cut << em_continuous
<< longprof << observationLevel;
-
+
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
@@ -232,7 +232,7 @@ int main(int argc, char** argv) {
EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile
- conexSource.SolveCE();
+ /* conexSource.SolveCE(); */
cut.ShowResults();
const HEPEnergyType Efinal =
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 1977467c2..582a59882 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -23,7 +23,10 @@ namespace corsika::process::proposal {
void ContinuousProcess::BuildCalculator(particles::Code code,
environment::NuclearComposition const& comp) {
- auto c = cross[code](media.at(&comp), cut);
+ auto p_cross = cross.find(code);
+ if (p_cross == cross.end())
+ throw std::runtime_error("PROPOSAL could not find corresponding builder");
+ auto c = p_cross->second(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)] =
@@ -83,13 +86,16 @@ namespace corsika::process::proposal {
template <>
units::si::LengthType ContinuousProcess::MaxStepLength(
setup::Stack::ParticleType const& vP, setup::Trajectory const& vT) {
+ if (!CanInteract(vP.GetPID()))
+ return units::si::meter * std::numeric_limits<double>::infinity();
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/Interaction.cc b/Processes/Proposal/Interaction.cc
index 219712d65..2ac0e7b70 100644
--- a/Processes/Proposal/Interaction.cc
+++ b/Processes/Proposal/Interaction.cc
@@ -31,7 +31,10 @@ namespace corsika::process::proposal {
void Interaction::BuildCalculator(particles::Code code,
environment::NuclearComposition const& comp) {
- auto c = cross[code](media.at(&comp), cut);
+ auto p_cross = cross.find(code);
+ if (p_cross == cross.end())
+ throw std::runtime_error("PROPOSAL could not find corresponding builder");
+ auto c = p_cross->second(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)),
diff --git a/Processes/Proposal/ProposalProcessBase.cc b/Processes/Proposal/ProposalProcessBase.cc
index 335cf4186..815f032b5 100644
--- a/Processes/Proposal/ProposalProcessBase.cc
+++ b/Processes/Proposal/ProposalProcessBase.cc
@@ -49,7 +49,7 @@ namespace corsika::process::proposal {
}
PROPOSAL::InterpolationDef::order_of_interpolation = 2;
PROPOSAL::InterpolationDef::nodes_cross_section = 100;
- PROPOSAL::InterpolationDef::nodes_propagate = 100;
+ PROPOSAL::InterpolationDef::nodes_propagate = 1000;
}
size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const noexcept {
--
GitLab