docs: def _calc_specpp_loE(self, Egamma):
          """
          Delta-functional approximation for low energies Egamma < 0.1 TeV
          """
          from scipy.integrate import quad
  
          Egamma = Egamma.to("TeV").value
          Epimin = Egamma + self._m_pi ** 2 / (4 * Egamma)
  
          result = (
              2
              * quad(
                  self._delta_integrand, Epimin, np.inf, epsrel=1e-3, epsabs=0
              )[0]
          )
  
          return result * u.Unit("1/(s TeV)")