Skip to content

Return discrete lines exactly in lin-lin tabular energy sampling #4069

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dallonby wants to merge 2 commits into
base: develop
Choose a base branch
from
+121 −100
Open

Return discrete lines exactly in lin-lin tabular energy sampling #4069

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dfd3434
Return discrete lines exactly in lin-lin tabular energy sampling
dallonby Aug 21, 2026
05fcd20
simplify code a bit
GuySten Aug 22, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 53 additions & 42 deletions src/distribution_energy.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -211,59 +211,70 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
}

double E_l_k = distribution_[l].e_out[k];
double p_l_k = distribution_[l].p[k];
double E_out = E_l_k;
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0 && k >= n_discrete) {
E_out = E_l_k + (r1 - c_k) / p_l_k;
}

} else if (distribution_[l].interpolation == Interpolation::lin_lin) {
// Linear-linear interpolation
double E_l_k1 = distribution_[l].e_out[k + 1];
double p_l_k1 = distribution_[l].p[k + 1];

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) {
if (k < n_discrete) {
// Discrete case
return E_l_k;
} else {
// Continuous case
double p_l_k = distribution_[l].p[k];
double E_out;
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 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;
E_out = E_l_k;
}
} else if (distribution_[l].interpolation == Interpolation::lin_lin) {
// Linear-linear interpolation
double E_l_k1 = distribution_[l].e_out[k + 1];
double p_l_k1 = distribution_[l].p[k + 1];

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 {
E_out = E_l_k;
}
} else {
throw std::runtime_error {
"Unexpected interpolation for continuous energy "
"distribution."};
}
} else {
throw std::runtime_error {"Unexpected interpolation for continuous energy "
"distribution."};
}

// Now interpolate between incident energy bins i and i + 1
if (!histogram_interp && n_energy_out > 1 && k >= n_discrete) {
// Interpolation for energy E1 and EK
n_energy_out = distribution_[i].e_out.size();
n_discrete = distribution_[i].n_discrete;
const double E_i_1 = distribution_[i].e_out[n_discrete];
const double E_i_K = distribution_[i].e_out[n_energy_out - 1];
// Now interpolate between incident energy bins i and i + 1
if (!histogram_interp && n_energy_out > 1) {
// Interpolation for energy E1 and EK
n_energy_out = distribution_[i].e_out.size();
n_discrete = distribution_[i].n_discrete;
const double E_i_1 = distribution_[i].e_out[n_discrete];
const double E_i_K = distribution_[i].e_out[n_energy_out - 1];

n_energy_out = distribution_[i + 1].e_out.size();
n_discrete = distribution_[i + 1].n_discrete;
const double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
n_energy_out = distribution_[i + 1].e_out.size();
n_discrete = distribution_[i + 1].n_discrete;
const double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];

const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
const double E_K = E_i_K + r * (E_i1_K - E_i_K);
const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
const double E_K = E_i_K + r * (E_i1_K - E_i_K);

if (l == i) {
return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
if (l == i) {
return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
} else {
return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
}
} else {
return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
return E_out;
}
} else {
return E_out;
}
}

Expand Down
55 changes: 29 additions & 26 deletions src/secondary_correlated.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -210,34 +210,37 @@ Distribution& CorrelatedAngleEnergy::sample_dist(
}

double E_l_k = distribution_[l].e_out[k];
double p_l_k = distribution_[l].p[k];
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0 && k >= n_discrete) {
E_out = E_l_k + (r1 - c_k) / p_l_k;
} else {
E_out = E_l_k;
}

} else if (distribution_[l].interpolation == Interpolation::lin_lin) {
// Linear-linear interpolation
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 (k < n_discrete) {
// Discrete case
E_out = E_l_k;
} else {
// Continuous case
double p_l_k = distribution_[l].p[k];
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0) {
E_out = E_l_k + (r1 - c_k) / p_l_k;
} else {
E_out = E_l_k;
}
} else if (distribution_[l].interpolation == Interpolation::lin_lin) {
// Linear-linear interpolation
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;
}
}
}

// Now interpolate between incident energy bins i and i + 1
if (k >= n_discrete) {
// Now interpolate between incident energy bins i and i + 1
if (l == i) {
E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
} else {
Expand Down
71 changes: 39 additions & 32 deletions src/secondary_kalbach.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,46 +169,53 @@ void KalbachMann::sample_params(
}

double E_l_k = distribution_[l].e_out[k];
double p_l_k = distribution_[l].p[k];
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0 && k >= n_discrete) {
E_out = E_l_k + (r1 - c_k) / p_l_k;
} else {
E_out = E_l_k;
}

// Determine Kalbach-Mann parameters
if (k < n_discrete) {
// Discrete case
E_out = E_l_k;
km_r = distribution_[l].r[k];
km_a = distribution_[l].a[k];

} else {
// Linear-linear interpolation
double E_l_k1 = distribution_[l].e_out[k + 1];
double p_l_k1 = distribution_[l].p[k + 1];
// Continuous case
double p_l_k = distribution_[l].p[k];
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0) {
E_out = E_l_k + (r1 - c_k) / p_l_k;
} else {
E_out = E_l_k;
}

// Determine Kalbach-Mann parameters
km_r = distribution_[l].r[k];
km_a = distribution_[l].a[k];

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;
}
// Linear-linear interpolation
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;
}

// Determine Kalbach-Mann parameters
km_r = distribution_[l].r[k] +
(E_out - E_l_k) / (E_l_k1 - E_l_k) *
(distribution_[l].r[k + 1] - distribution_[l].r[k]);
km_a = distribution_[l].a[k] +
(E_out - E_l_k) / (E_l_k1 - E_l_k) *
(distribution_[l].a[k + 1] - distribution_[l].a[k]);
}
// Determine Kalbach-Mann parameters
km_r = distribution_[l].r[k] +
(E_out - E_l_k) / (E_l_k1 - E_l_k) *
(distribution_[l].r[k + 1] - distribution_[l].r[k]);
km_a = distribution_[l].a[k] +
(E_out - E_l_k) / (E_l_k1 - E_l_k) *
(distribution_[l].a[k + 1] - distribution_[l].a[k]);
}

// Now interpolate between incident energy bins i and i + 1
if (k >= n_discrete) {
// Now interpolate between incident energy bins i and i + 1
if (l == i) {
E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
} else {
Expand Down
Loading