IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 370d68ad authored by Maximilian Sackel's avatar Maximilian Sackel
Browse files

add type checking before processing particle with proposal

parent 7b37877b
No related branches found
No related tags found
No related merge requests found
...@@ -188,7 +188,7 @@ int main(int argc, char** argv) { ...@@ -188,7 +188,7 @@ int main(int argc, char** argv) {
PROPOSAL::InterpolationDef::path_to_tables = "~/.local/share/PROPOSAL/tables/"; PROPOSAL::InterpolationDef::path_to_tables = "~/.local/share/PROPOSAL/tables/";
PROPOSAL::InterpolationDef::path_to_tables_readonly = "~/.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::Interaction proposal(env, cut);
process::proposal::ContinuousProcess em_continuous(env, cut); process::proposal::ContinuousProcess em_continuous(env, cut);
process::interaction_counter::InteractionCounter proposalCounted(proposal); process::interaction_counter::InteractionCounter proposalCounted(proposal);
...@@ -205,9 +205,9 @@ int main(int argc, char** argv) { ...@@ -205,9 +205,9 @@ int main(int argc, char** argv) {
process::UrQMD::UrQMD urqmd; process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd}; process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
process::conex_source_cut::CONEXSourceCut conexSource( /* process::conex_source_cut::CONEXSourceCut conexSource( */
center, showerAxis, t, injectionHeight, E0, /* center, showerAxis, t, injectionHeight, E0, */
particles::GetPDG(particles::Code::Proton)); /* particles::GetPDG(particles::Code::Proton)); */
// assemble all processes into an ordered process list // assemble all processes into an ordered process list
...@@ -218,9 +218,9 @@ int main(int argc, char** argv) { ...@@ -218,9 +218,9 @@ int main(int argc, char** argv) {
//auto sequence = switchProcess << reset_particle_mass << decaySequence << conexSource //auto sequence = switchProcess << reset_particle_mass << decaySequence << conexSource
// << longprof << eLoss << cut << observationLevel; // << 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; << longprof << observationLevel;
// define air shower object, run simulation // define air shower object, run simulation
tracking_line::TrackingLine tracking; tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack); cascade::Cascade EAS(env, tracking, sequence, stack);
...@@ -232,7 +232,7 @@ int main(int argc, char** argv) { ...@@ -232,7 +232,7 @@ int main(int argc, char** argv) {
EAS.Run(); EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile eLoss.PrintProfile(); // print longitudinal profile
conexSource.SolveCE(); /* conexSource.SolveCE(); */
cut.ShowResults(); cut.ShowResults();
const HEPEnergyType Efinal = const HEPEnergyType Efinal =
......
...@@ -23,7 +23,10 @@ namespace corsika::process::proposal { ...@@ -23,7 +23,10 @@ namespace corsika::process::proposal {
void ContinuousProcess::BuildCalculator(particles::Code code, void ContinuousProcess::BuildCalculator(particles::Code code,
environment::NuclearComposition const& comp) { 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 disp = PROPOSAL::make_displacement(c, true);
auto scatter = PROPOSAL::make_scattering("highland", particle[code], media.at(&comp)); auto scatter = PROPOSAL::make_scattering("highland", particle[code], media.at(&comp));
calc[std::make_pair(&comp, code)] = calc[std::make_pair(&comp, code)] =
...@@ -83,13 +86,16 @@ namespace corsika::process::proposal { ...@@ -83,13 +86,16 @@ namespace corsika::process::proposal {
template <> template <>
units::si::LengthType ContinuousProcess::MaxStepLength( units::si::LengthType ContinuousProcess::MaxStepLength(
setup::Stack::ParticleType const& vP, setup::Trajectory const& vT) { 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(); auto energy_lim = 0.9 * vP.GetEnergy();
if (cut.GetECut() > energy_lim) energy_lim = cut.GetECut(); if (cut.GetECut() > energy_lim) energy_lim = cut.GetECut();
auto c = GetCalculator(vP, calc); auto c = GetCalculator(vP, calc);
auto grammage = get<DISPLACEMENT>(c->second)->SolveTrackIntegral( auto grammage = get<DISPLACEMENT>(c->second)->SolveTrackIntegral(
vP.GetEnergy() / 1_MeV, energy_lim / 1_MeV) * vP.GetEnergy() / 1_MeV, energy_lim / 1_MeV) *
1_g / square(1_cm); 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 } // namespace corsika::process::proposal
...@@ -31,7 +31,10 @@ namespace corsika::process::proposal { ...@@ -31,7 +31,10 @@ namespace corsika::process::proposal {
void Interaction::BuildCalculator(particles::Code code, void Interaction::BuildCalculator(particles::Code code,
environment::NuclearComposition const& comp) { 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); auto inter_types = PROPOSAL::CrossSectionVector::GetInteractionTypes(c);
calc[std::make_pair(&comp, code)] = std::make_tuple( calc[std::make_pair(&comp, code)] = std::make_tuple(
PROPOSAL::make_secondaries(inter_types, particle[code], media.at(&comp)), PROPOSAL::make_secondaries(inter_types, particle[code], media.at(&comp)),
......
...@@ -49,7 +49,7 @@ namespace corsika::process::proposal { ...@@ -49,7 +49,7 @@ namespace corsika::process::proposal {
} }
PROPOSAL::InterpolationDef::order_of_interpolation = 2; PROPOSAL::InterpolationDef::order_of_interpolation = 2;
PROPOSAL::InterpolationDef::nodes_cross_section = 100; 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 { size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const noexcept {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment