Add Polynomial Interpolation to Speed Up Total Xsec

[tomeichlersmith]

While the details of the total cross section calculation are incredibly complicated, the resulting function of total cross section vs energy within a specific material is relatively smooth. This begs for a performance improvement by only doing the full calculation at a few intelligently selected points and then interpolating between them.

DMG4 does implement this idea; however, the points that are sampled are hard-coded and chosen at compile time. I'd like to implement a more dynamic scheme where the sample points are chosen while the process is running so that the interpolation is "specialized" to the specific running case.

[tomeichlersmith]

https://kluge.in-chemnitz.de/opensource/spline/

[tomeichlersmith]

Our cross section function is smooth and monotonic, so we can look at implementing https://en.wikipedia.org/wiki/Monotone_cubic_interpolation which calls the underlying cross section computation if the energy requested is outside of the bounds of interpolation.

[tomeichlersmith]

Since our function is expected to be so smooth, I'm going to avoid more general interpolation algorithms and instead specialize a simple interpolation algorithm. We expect the function to roughly follow $y \sim \sqrt{x-a}$ so I am going to use a quadratic-order interpolation. We can simplify the generalized polynomial interpolation to the $n=2$ case. In our specific case, we know that $a = 2 m{A'}$, so - from here on out - assume $x = E - 2 m{A'}$ so we don't have to worry about going below the kinematic minimum.

$$ p(x) = \left(\frac{x-x_1}{x_0-x_1}\right)\left(\frac{x-x_2}{x_0-x_2}\right) y_0 + \left(\frac{x-x_0}{x_1-x_0}\right)\left(\frac{x-x_2}{x_1-x_2}\right) y_1 + \left(\frac{x-x_0}{x_2-x_0}\right)\left(\frac{x-x_1}{x_2-x_1}\right) y_2 $$

We can use this formula for $x_0 <= x <= x_2$, so let's assume we start with the ordered set $(x_k)$ of size $n = 3$ that fit that criteria. If a new $x$ arrives that is outside that range, we can expand that range to include it by sampling below $x_0$ or above $x_n$ whichever is required.

How do we choose where to sample?!? Some ideas (focus on going lower since we expect that to be used more frequently, going higher can be deduced similarly.) Assume $x < x_0$ where $x_0$ is the lowest currently-sampled point. Changing a sampled point expands the set $(x_k)$ shifting indices.

• $x_0 \to x_0 / 2$ until $x_0 < x$
  • could allow the parameter $2$ to be chosen by user
• $x_0 \to \sqrt{x_0}$ until $x_0 < x$
  • works as long as all x are above 1
• $x_0 \to x / 2$ and $x_1 \to x$
  • Gets us a "free" perfect evaluation and maybe more accurate of an interpolation since we physically expect many particles to follow the same energy-loss path as defined by earlier particles in the run?