diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index d78f1f62478d631f3aaeb44431df3321bfff8d4a..26b259c84434178308394fe3e302e04ae3fe9662 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 1977467c2f91ab5513c2d491e2f468ec4ff67330..582a59882b70a3aa3aa35fbd4be03379b0368d59 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 219712d65c6f950be8d616de50e4270ae9255e03..2ac0e7b70b4f256a10dd2b79f058d1fff03943ab 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 335cf41860160dbbd73df8e79dcd94535a306eef..815f032b533de2aca85587dd036d7d3a43b20417 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 {