Smoothing interpolators: pair production (atomic nucleus and electron), log-log splines

[lynntf]

Current behavior for interpolation with pair production does not appear to use a correct linearization function for the data. Although this seems to be an accurate reproduction of the available XCOM software and what is explained in the XCOM documents, the interpolation is not behaving as it should.

Consider the output of the following using the current code (adapted from the XCOM example notebook):

import numpy as np
import xcom
import matplotlib.pyplot as plt

def plot(data):
    x = data["energy"]/1e6
    plt.plot(x, data["coherent"], label="Coherent")
    plt.plot(x, data["incoherent"], label="Incoherent")
    plt.plot(x, data["photoelectric"], label="Photoelectric")
    plt.plot(x, data["pair_atom"], label="Pair atom")
    plt.plot(x, data["pair_electron"], label="Pair electron")
    plt.plot(x, data["total"], label="Total")
    plt.xlabel("MeV")
    plt.ylabel("barn/atom")
    plt.xscale("log")
    plt.yscale("log")
    plt.legend()

# Interpolate between tabulated datapoints
plt.figure()
energies = np.logspace(3,11,500)
interpolated_data_92 = xcom.calculate_cross_section("U", energies)
plot(interpolated_data_92)

Notice the large fluctuations near the pair production energy thresholds.

• XCOM documents indicate that log-log interpolation for pair production should be done with $\left(1-\frac{E}{E'}\right)^3 \sigma_{\mathrm{PAIR}}(E) = \left(\frac{E' - E}{E'}\right)^3 \sigma_{\mathrm{PAIR}}(E)$ where $E'$ is the pair production threshold. This is what is implemented in the current nist-calculators/xcom code.
• XCOM source code indicates that the log-log interpolation for pair production is done with:

        TERM=(E(M)-EPAIR1)/E(M) 
    420 PR1(M)=LOG(PAIRAT(M)/(TERM**3)) 

  That is, $\left(1 - \frac{E'}{E}\right)^{-3} \sigma_{\mathrm{PAIR}}(E) = \left(\frac{E}{E - E'}\right)^3 \sigma_{\mathrm{PAIR}}(E)$.
• Unfortunately, I am not able to see the implementation used by the web interface, although I suspect that it is neither of these as the interpolated cross-sections appear smooth.

These are not equivalent, but similar. However, both seem to be incorrect for creating a linear log-log plot for the log-log interpolation.

The normalizing function appears to have been mistaken at some point and the corrected linearizing function implemented now is $\left(E (E - E')\right)^{-3} \sigma_{\mathrm{PAIR}}(E)$. Perhaps at some point this was written as $\sigma_{\mathrm{PAIR}}(E)/E^3/(E - E')^3$ and incorrectly interpreted as $\sigma_{\mathrm{PAIR}}(E)/\left(E^3/(E - E')^3\right)$, but it is not clear.

Additional changes were made to the boundary condition types for the interpolators (allowing for linear extrapolation in the log-log domain; bc_type='natural') along with some other linting changes and some more explicit code.

[lynntf]

Some other changes have been made:

• It appears that the material mixing had a bad index resulting in incorrect computations of mixed attenuation coefficients.
• The 'energy' data was being rescaled when it should not have been.
• Material mixing by formula (e.g., H2O 0.7 Bi4Ge3O12 0.3 as would be allowed on the XCOM web interface is not allowed)
• Not all of the oscillations of the CubicSpline interpolation were dealt with by changing the boundary conditions to 'natural', so the interpolator has been replaced with a PchipInterpolator. This may introduce problems when dμ/dE = 0 (Pchip is monotonic and may not handle the change from increasing to decreasing present in incoherent scattering).
• Linear extrapolation has been added (optional warnings when it is being used also added)

[lynntf]

Another note based on #7, computation of total_without_coherent is fixed here.