IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 77404e91 authored by Felix Riehn's avatar Felix Riehn Committed by ralfulrich
Browse files

removed explicit threshold from proposal

parent 284a7906
No related branches found
No related tags found
No related merge requests found
......@@ -28,11 +28,13 @@ namespace corsika::proposal {
auto p_cross = cross.find(code);
if (p_cross == cross.end())
throw std::runtime_error("PROPOSAL could not find corresponding builder");
// interpolate the crosssection for given media and energy cut. These may
// take some minutes if you have to build the tables and cannot read the
// from disk
auto c = p_cross->second(media.at(comp.getHash()), emCut_);
auto const emCut = get_energy_threshold(
code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Build displacement integral and scattering object and interpolate them too and
// saved in the calc map by a key build out of a hash of composed of the component and
......@@ -45,9 +47,9 @@ namespace corsika::proposal {
}
template <>
ContinuousProcess::ContinuousProcess(setup::Environment const& _env,
HEPEnergyType _emCut)
: ProposalProcessBase(_env, _emCut) {}
ContinuousProcess::ContinuousProcess(setup::Environment const& _env
)
: ProposalProcessBase(_env) {}
template <>
void ContinuousProcess::scatter(setup::Stack::particle_type& vP,
......@@ -115,21 +117,24 @@ namespace corsika::proposal {
template <>
LengthType ContinuousProcess::getMaxStepLength(setup::Stack::particle_type const& vP,
setup::Trajectory const& vT) {
if (!canInteract(vP.getPID())) return meter * std::numeric_limits<double>::infinity();
auto const code = vP.getPID();
if (!canInteract(code)) return meter * std::numeric_limits<double>::infinity();
// Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a
// hyper parameter which must be adjusted.
//
// in any case: never go below 0.99*emCut_ This needs to be
// slightly smaller than emCut_ since, either this Step is limited
auto const emCut = get_energy_threshold(
code); //! energy thresholds globally defined for individual particles
// in any case: never go below 0.99*emCut This needs to be
// slightly smaller than emCut since, either this Step is limited
// by energy_lim, then the particle is stopped in a very short
// range (before doing anythin else) and is then removed
// instantly. The exact position where it reaches emCut is not
// important, the important fact is that its E_kin is zero
// afterwards.
//
auto energy_lim = std::max(0.9 * vP.getEnergy(), 0.99 * emCut_);
auto energy_lim = std::max(0.9 * vP.getEnergy(), 0.99 * emCut);
// solving the track integral for giving energy lim
auto c = getCalculator(vP, calc);
......
......@@ -24,8 +24,8 @@
namespace corsika::proposal {
template <>
Interaction::Interaction(setup::Environment const& _env, HEPEnergyType _emCut)
: ProposalProcessBase(_env, _emCut) {}
Interaction::Interaction(setup::Environment const& _env)
: ProposalProcessBase(_env) {}
void Interaction::buildCalculator(Code code, NuclearComposition const& comp) {
// search crosssection builder for given particle
......@@ -36,7 +36,10 @@ namespace corsika::proposal {
// interpolate the crosssection for given media and energy cut. These may
// take some minutes if you have to build the tables and cannot read the
// from disk
auto c = p_cross->second(media.at(comp.getHash()), emCut_);
auto const emCut = get_energy_threshold(
code); //! energy thresholds globally defined for individual particles
auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Look which interactions take place and build the corresponding
// interaction and secondarie builder. The interaction integral will
......
......@@ -30,10 +30,9 @@ namespace corsika::proposal {
return false;
}
ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env,
HEPEnergyType _emCut)
: emCut_(_emCut)
, RNG_(RNGManager::getInstance().getRandomStream("proposal")) {
ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env
)
: RNG_(RNGManager::getInstance().getRandomStream("proposal")) {
_env.getUniverse()->walk([&](auto& vtn) {
if (vtn.hasModelProperties()) {
const auto& prop = vtn.getModelProperties();
......
......@@ -52,7 +52,7 @@ namespace corsika::proposal {
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
ContinuousProcess(TEnvironment const&, HEPEnergyType _emCut);
ContinuousProcess(TEnvironment const&);
//!
//! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into
......
......@@ -48,7 +48,7 @@ namespace corsika::proposal {
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
Interaction(TEnvironment const& env, HEPEnergyType emCut);
Interaction(TEnvironment const& env);
//!
//! Calculate the rates for the different targets and interactions. Sample a
......
......@@ -44,7 +44,11 @@ namespace corsika::proposal {
//!
template <typename T>
static auto cross_builder =
[](PROPOSAL::Medium& m, corsika::units::si::HEPEnergyType emCut) {
[](PROPOSAL::Medium& m,
corsika::units::si::HEPEnergyType
emCut) { //!< Stochastic losses smaller than the given cut
//!< will be handeled continuously.
using namespace corsika::units::si;
auto p_cut =
std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, true);
......@@ -73,8 +77,6 @@ namespace corsika::proposal {
//!
class ProposalProcessBase {
protected:
HEPEnergyType emCut_; //!< Stochastic losses smaller than the given cut
//!< will be handeled continuously.
RNGManager::prng_type RNG_; //!< random number generator used by proposal
std::unordered_map<std::size_t, PROPOSAL::Medium>
......@@ -85,7 +87,7 @@ namespace corsika::proposal {
//! Store cut and nuclear composition of the whole universe in media which are
//! required for creating crosssections by proposal.
//!
ProposalProcessBase(corsika::setup::Environment const& _env, HEPEnergyType _emCut);
ProposalProcessBase(corsika::setup::Environment const& _env);
//!
//! Checks if a particle can be processed by proposal
......
......@@ -143,8 +143,8 @@ int main(int argc, char** argv) {
// PROPOSAL processs proposal{...};
ParticleCut cut(10_GeV, 10_GeV, 100_PeV, 100_PeV, true);
corsika::proposal::Interaction proposal(env, cut.getElectronECut());
corsika::proposal::ContinuousProcess em_continuous(env, cut.getElectronECut());
corsika::proposal::Interaction proposal(env);
corsika::proposal::ContinuousProcess em_continuous(env);
InteractionCounter proposalCounted(proposal);
TrackWriter trackWriter("tracks.dat");
......
......@@ -222,8 +222,8 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
corsika::proposal::Interaction proposal(env, cut.getElectronECut());
corsika::proposal::ContinuousProcess em_continuous(env, cut.getElectronECut());
corsika::proposal::Interaction proposal(env);
corsika::proposal::ContinuousProcess em_continuous(env);
InteractionCounter proposalCounted(proposal);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
......
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