Make a quantile-quantile plot comparing multiple sets of events and models.

    *arrays*
      X.
    *key_text*
      Text describing the quantile-quantile comparison quantity; will be
      shown on the plot legend.
    Returns:
      An :class:`omega.RectPlot` instance.

    *TODO*: nothing about this is Sherpa-specific. Same goes for some of the
    plotting routines in :mod:`pkwit.environments.casa.data`; might be
    reasonable to add a submodule for generic X-ray-y plotting routines.

    *TODO*: Some gross code duplication here.

def make_multi_qq_plots(arrays, key_text):
    """Make a quantile-quantile plot comparing multiple sets of events and models.

    *arrays*
      X.
    *key_text*
      Text describing the quantile-quantile comparison quantity; will be
      shown on the plot legend.
    Returns:
      An :class:`omega.RectPlot` instance.

    *TODO*: nothing about this is Sherpa-specific. Same goes for some of the
    plotting routines in :mod:`pkwit.environments.casa.data`; might be
    reasonable to add a submodule for generic X-ray-y plotting routines.

    *TODO*: Some gross code duplication here.

    """
    import omega as om

    p = om.RectPlot()
    p.addXY([0, 1.], [0, 1.], '1:1')

    for index, array in enumerate(arrays):
        kev, obs, mdl = array
        c_obs = np.cumsum(obs)
        c_mdl = np.cumsum(mdl)

        mx = 0.5 * (c_obs[-1] + c_mdl[-1])
        c_obs /= mx
        c_mdl /= mx

        p.addXY(c_mdl, c_obs, '%s #%d' % (key_text, index))

    # HACK: this range of numbers is chosen to give reasonable sampling for my
    # sources, which are typically quite soft.
    #
    # Note: this reuses the variables from the last loop iteration.
    locs = np.array([0, 0.05, 0.08, 0.11, 0.17, 0.3, 0.4, 0.7, 1]) * (kev.size - 2)
    c0 = 1.05
    c1 = 1.1

    for loc in locs:
        i0 = int(np.floor(loc))
        frac = loc - i0
        kevval = (1 - frac) * kev[i0] + frac * kev[i0+1]
        mdlval = (1 - frac) * c_mdl[i0] + frac * c_mdl[i0+1]
        obsval = (1 - frac) * c_obs[i0] + frac * c_obs[i0+1]
        p.addXY([mdlval, mdlval], [c0, c1], '%.2f keV' % kevval, dsn=2)
        p.addXY([c0, c1], [obsval, obsval], None, dsn=2)

    p.setLabels('Cumulative rescaled model', 'Cumulative rescaled data')
    p.defaultKeyOverlay.vAlign = 0.3
    return p