Make a plot of multiple spectral models and data.

    *model_plot*
      A model plot object returned by Sherpa from a call like
      ``ui.get_model_plot()`` or ``ui.get_bkg_model_plot()``.
    *data_plots*
      An iterable of data plot objects returned by Sherpa from calls like
      ``ui.get_source_plot(id)`` or ``ui.get_bkg_plot(id)``.
    *desc*
      Text describing the origin of the data; will be shown in the plot legend
      (with "Model" and "Data #<number>" appended).
    *xmin_clamp*
      The smallest "x" (energy axis) value that will be plotted; default is 0.01.
      This is needed to allow the plot to be shown on a logarithmic scale if
      the energy axes of the model go all the way to 0.
    *min_valid_x*
      Either None, or the smallest "x" (energy axis) value in which the model and
      data are valid; this could correspond to a range specified in the "notice"
      command during analysis. If specified, a gray band will be added to the plot
      showing the invalidated regions.
    *max_valid_x*
      Like *min_valid_x* but for the largest "x" (energy axis) value in which the
      model and data are valid.
    Returns:
      A tuple ``(plot, xlow, xhigh)``, where *plot* an OmegaPlot RectPlot
      instance, *xlow* is the left edge of the plot bounds, and *xhigh* is the
      right edge of the plot bounds.

    TODO: not happy about the code duplication with :func:`make_spectrum_plot`
    but here we are.

def make_multi_spectrum_plots(model_plot, plotids, data_getter, desc, xmin_clamp=0.01,
                              min_valid_x=None, max_valid_x=None):
    """Make a plot of multiple spectral models and data.

    *model_plot*
      A model plot object returned by Sherpa from a call like
      ``ui.get_model_plot()`` or ``ui.get_bkg_model_plot()``.
    *data_plots*
      An iterable of data plot objects returned by Sherpa from calls like
      ``ui.get_source_plot(id)`` or ``ui.get_bkg_plot(id)``.
    *desc*
      Text describing the origin of the data; will be shown in the plot legend
      (with "Model" and "Data #<number>" appended).
    *xmin_clamp*
      The smallest "x" (energy axis) value that will be plotted; default is 0.01.
      This is needed to allow the plot to be shown on a logarithmic scale if
      the energy axes of the model go all the way to 0.
    *min_valid_x*
      Either None, or the smallest "x" (energy axis) value in which the model and
      data are valid; this could correspond to a range specified in the "notice"
      command during analysis. If specified, a gray band will be added to the plot
      showing the invalidated regions.
    *max_valid_x*
      Like *min_valid_x* but for the largest "x" (energy axis) value in which the
      model and data are valid.
    Returns:
      A tuple ``(plot, xlow, xhigh)``, where *plot* an OmegaPlot RectPlot
      instance, *xlow* is the left edge of the plot bounds, and *xhigh* is the
      right edge of the plot bounds.

    TODO: not happy about the code duplication with :func:`make_spectrum_plot`
    but here we are.

    """
    import omega as om
    from omega.stamps import DataThemedStamp, WithYErrorBars

    model_x = np.concatenate((model_plot.xlo, [model_plot.xhi[-1]]))
    model_x[0] = max(model_x[0], xmin_clamp)
    model_y = np.concatenate((model_plot.y, [0.]))

    # Sigh, sometimes Sherpa gives us bad values.
    is_bad = ~np.isfinite(model_y)
    if is_bad.sum():
        from .cli import warn
        warn('bad Sherpa model Y value(s) at: %r', np.where(is_bad)[0])
        model_y[is_bad] = 0

    p = om.RectPlot()
    data_csps = []
    data_lines = []
    xlow = xhigh = None

    for index, plotid in enumerate(plotids):
        data_plot = data_getter(plotid)
        data_left_edges = data_plot.x - 0.5 * data_plot.xerr
        data_left_edges[0] = max(data_left_edges[0], xmin_clamp)
        data_hist_x = np.concatenate((data_left_edges, [data_plot.x[-1] + 0.5 * data_plot.xerr[-1]]))
        data_hist_y = np.concatenate((data_plot.y, [0.]))

        if xlow is None:
            xlow = model_x[0]
            xhigh = model_x[-1]
        else:
            xlow = min(xlow, model_x[0])
            xhigh = max(xhigh, model_x[-1])

        csp = om.rect.ContinuousSteppedPainter(keyText=None)
        csp.setFloats(data_hist_x, data_hist_y)
        data_csps.append(csp)

        inner_stamp = DataThemedStamp(None)
        stamp = WithYErrorBars(inner_stamp)

        lines = om.rect.XYDataPainter(
            lines = False,
            pointStamp = stamp,
            keyText = '%s Data #%d' % (desc, index)
        )
        lines.setFloats(data_plot.x, data_plot.y,
                        data_plot.y + data_plot.yerr,
                        data_plot.y - data_plot.yerr)
        inner_stamp.setHolder(lines)
        data_lines.append(lines)

    log_bounds_pad_factor = 0.9
    xlow *= log_bounds_pad_factor
    xhigh /= log_bounds_pad_factor

    if min_valid_x is not None:
        p.add(om.rect.XBand(1e-3 * xlow, min_valid_x, keyText=None), zheight=-1, dsn=1)
    if max_valid_x is not None:
        p.add(om.rect.XBand(max_valid_x, xhigh * 1e3, keyText=None), zheight=-1, dsn=1)

    model_csp = om.rect.ContinuousSteppedPainter(keyText=desc + ' Model')
    model_csp.setFloats(model_x, model_y)
    p.add(model_csp)

    for index, (data_csp, lines) in enumerate(zip(data_csps, data_lines)):
        p.add(data_csp, dsn=index + 1)
        p.add(lines, dsn=index + 1)

    p.setLabels(data_plot.xlabel, data_plot.ylabel) # data_plot = last one from the for loop
    p.setLinLogAxes(True, False)
    p.setBounds (xlow, xhigh)
    return p, xlow, xhigh