docs: def _plot_data_to_ax(
      data_all,
      ax1,
      e_unit=None,
      sed=True,
      ylabel=None,
      ulim_opts={},
      errorbar_opts={},
  ):
      """ Plots data errorbars and upper limits onto ax.
      X label is left to plot_data and plot_fit because they depend on whether
      residuals are plotted.
      """
  
      if e_unit is None:
          e_unit = data_all["energy"].unit
  
      f_unit, sedf = sed_conversion(
          data_all["energy"], data_all["flux"].unit, sed
      )
  
      if "group" not in data_all.keys():
          data_all["group"] = np.zeros(len(data_all))
  
      groups = np.unique(data_all["group"])
  
      for g in groups:
          data = data_all[np.where(data_all["group"] == g)]
          _, sedfg = sed_conversion(data["energy"], data["flux"].unit, sed)
  
          # wrap around color and marker cycles
          color = color_cycle[int(g) % len(color_cycle)]
          marker = marker_cycle[int(g) % len(marker_cycle)]
  
          ul = data["ul"]
          notul = ~ul
  
          # Hack to show y errors compatible with 0 in loglog plot
          yerr_lo = data["flux_error_lo"][notul]
          y = data["flux"][notul].to(yerr_lo.unit)
          bad_err = np.where((y - yerr_lo) <= 0.0)
          yerr_lo[bad_err] = y[bad_err] * (1.0 - 1e-7)
          yerr = u.Quantity((yerr_lo, data["flux_error_hi"][notul]))
          xerr = u.Quantity((data["energy_error_lo"], data["energy_error_hi"]))
  
          opts = dict(
              zorder=100,
              marker=marker,
              ls="",
              elinewidth=2,
              capsize=0,
              mec=color,
              mew=0.1,
              ms=5,
              color=color,
          )
          opts.update(**errorbar_opts)
  
          ax1.errorbar(
              data["energy"][notul].to(e_unit).value,
              (data["flux"][notul] * sedfg[notul]).to(f_unit).value,
              yerr=(yerr * sedfg[notul]).to(f_unit).value,
              xerr=xerr[:, notul].to(e_unit).value,
              **opts
          )
  
          if np.any(ul):
              if "elinewidth" in errorbar_opts:
                  ulim_opts["elinewidth"] = errorbar_opts["elinewidth"]
  
              _plot_ulims(
                  ax1,
                  data["energy"][ul].to(e_unit).value,
                  (data["flux"][ul] * sedfg[ul]).to(f_unit).value,
                  (xerr[:, ul]).to(e_unit).value,
                  color,
                  **ulim_opts
              )
  
      ax1.set_xscale("log")
      ax1.set_yscale("log")
      xmin = 10 ** np.floor(
          np.log10(
              np.min(data["energy"] - data["energy_error_lo"]).to(e_unit).value
          )
      )
      xmax = 10 ** np.ceil(
          np.log10(
              np.max(data["energy"] + data["energy_error_hi"]).to(e_unit).value
          )
      )
      ax1.set_xlim(xmin, xmax)
      # avoid autoscaling to errorbars to 0
      notul = ~data_all["ul"]
      if np.any(data_all["flux_error_lo"][notul] >= data_all["flux"][notul]):
          elo = (data_all["flux"][notul] * sedf[notul]).to(f_unit).value - (
              data_all["flux_error_lo"][notul] * sedf[notul]
          ).to(f_unit).value
          gooderr = np.where(
              data_all["flux_error_lo"][notul] < data_all["flux"][notul]
          )
          ymin = 10 ** np.floor(np.log10(np.min(elo[gooderr])))
          ax1.set_ylim(bottom=ymin)
  
      if ylabel is None:
          if sed:
              ax1.set_ylabel(
                  r"$E^2\mathrm{{d}}N/\mathrm{{d}}E$"
                  " [{0}]".format(u.Unit(f_unit).to_string("latex_inline"))
              )
          else:
              ax1.set_ylabel(
                  r"$\mathrm{{d}}N/\mathrm{{d}}E$"
                  " [{0}]".format(u.Unit(f_unit).to_string("latex_inline"))
              )
      else:
          ax1.set_ylabel(ylabel)