docs: def _calc_specpp_hiE(self, Egamma):
          """
          Spectrum computed as in Eq. 42 for Egamma >= 0.1 TeV
          """
          # Fixed quad with n=40 is about 15 times faster and is always within
          # 0.5% of the result of adaptive quad for Egamma>0.1
          # WARNING: It also produces artifacts for steep distributions (e.g.
          # Maxwellian) at ~500 GeV. Reverting to adaptative quadrature
          # from scipy.integrate import fixed_quad
          # result=c*fixed_quad(self._photon_integrand, 0., 1., args = [Egamma,
          # ], n = 40)[0]
          from scipy.integrate import quad
  
          Egamma = Egamma.to("TeV").value
          specpp = (
              c.cgs.value
              * quad(
                  self._photon_integrand,
                  0.0,
                  1.0,
                  args=Egamma,
                  epsrel=1e-3,
                  epsabs=0,
              )[0]
          )
  
          return specpp * u.Unit("1/(s TeV)")