docs: def _spectrum(self, photon_energy):
          """
          Compute differential spectrum from pp interactions using the
          parametrization of Kafexhiu, E., Aharonian, F., Taylor, A.~M., and
          Vila, G.~S.\ 2014, `arXiv:1406.7369
          <http://www.arxiv.org/abs/1406.7369>`_.
  
          Parameters
          ----------
          photon_energy : :class:`~astropy.units.Quantity` instance
              Photon energy array.
          """
  
          # Load LUT if available, otherwise use self._diffsigma
          if self.useLUT:
              LUT_base = "PionDecayKafexhiu14_LUT_"
              if self.nuclear_enhancement:
                  LUT_base += "NucEnh_"
              LUT_fname = LUT_base + "{0}.npz".format(self.hiEmodel)
              # only reload LUT if it has changed or hasn't been loaded yet
              try:
                  if os.path.basename(self.diffsigma.fname) != LUT_fname:
                      self._loadLUT(LUT_fname)
              except AttributeError:
                  self._loadLUT(LUT_fname)
          else:
              self.diffsigma = self._diffsigma
  
          Egamma = _validate_ene(photon_energy).to("GeV")
          Ep = self._Ep * u.GeV
          J = self._J * u.Unit("1/GeV")
  
          specpp = []
          for Eg in Egamma:
              diffsigma = self.diffsigma(Ep.value, Eg.value) * u.Unit("cm2/GeV")
              specpp.append(trapz_loglog(diffsigma * J, Ep))
  
          self.specpp = u.Quantity(specpp)
  
          self.specpp *= self.nh * c.cgs
  
          return self.specpp.to("1/(s eV)")