Prevent potential divide-by-zero in ContinuousTabular::sample

This commit is contained in:
Paul Romano 2020-06-29 09:23:15 -05:00
parent 1b0400d9fe
commit 1319031a0f

View file

@ -238,12 +238,14 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
double E_l_k1 = distribution_[l].e_out[k+1];
double p_l_k1 = distribution_[l].p[k+1];
double frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k);
if (frac == 0.0) {
E_out = E_l_k + (r1 - c_k)/p_l_k;
} else {
E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k +
2.0*frac*(r1 - c_k))) - p_l_k)/frac;
if (E_l_k != E_l_k1) {
double frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k);
if (frac == 0.0) {
E_out = E_l_k + (r1 - c_k)/p_l_k;
} else {
E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k +
2.0*frac*(r1 - c_k))) - p_l_k)/frac;
}
}
} else {
throw std::runtime_error{"Unexpected interpolation for continuous energy "