numerical issue in O16 CC cross sections

[JostMigenda]

Discovered this issue while working on a PR for SNOwGLoBES: The oxygen-16 CC cross sections are implemented using a fit from arXiv:1809.08398, that has four "excitation energy" parameters. If the neutrino energy is barely higher than one of these excitation energies, the cross section has an unphysical spike (because the fit function diverges as the difference between neutrino energy and excitation energy goes to zero).

Sample code to demonstrate:

from sntools.interaction_channels import o16e
for i in range(1,5):
    eG, a, b, c = o16e.fit_parameters[i]
    eNu = 29.350001 # excitation energy of the 4th group: 29.35 MeV
    print(f"*** {i}: eNu={eNu}, eG={eG} ***")
    if eNu - eG > 0:
        d = log10(eNu**0.25 - eG**0.25)
        log_sigma = a + b * d + c * d**2
        print(f"d = {d}\nlog_sigma = {log_sigma}\npartial_dSigma_dE(eNu, eNu-eG, {i}) = {o16e.partial_dSigma_dE(eNu, eNu - eG,i)}\n")

In reality, the difference between eNu and eG needs to be above the detector threshold (e.g. 0.8 MeV in water Cherenkov detectors). Enforcing that may be sufficient to fix this?