def compare_SED(castroData1, castroData2, ylims, TS_thresh=4.0,
                errSigma=1.0, specVals=[]):
    """ Compare two SEDs 

    castroData1: A CastroData object, with the 
                 log-likelihood v. normalization for each energy bin
    castroData2: A CastroData object, with the 
                 log-likelihood v. normalization for each energy bin
    ylims      : y-axis limits
    TS_thresh  : TS value above with to plot a point, 
                 rather than an upper limit
    errSigma   : Number of sigma to use for error bars
    specVals   : List of spectra to add to plot
    returns fig,ax which are matplotlib figure and axes objects
    """
    import matplotlib.pyplot as plt

    xmin = min(castroData1.refSpec.ebins[0], castroData2.refSpec.ebins[0])
    xmax = max(castroData1.refSpec.ebins[-1], castroData2.refSpec.ebins[-1])
    ymin = ylims[0]
    ymax = ylims[1]

    fig = plt.figure()
    ax = fig.add_subplot(111)

    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_xlim((xmin, xmax))
    ax.set_ylim((ymin, ymax))

    ax.set_xlabel("Energy [GeV]")
    ax.set_ylabel(NORM_LABEL[castroData1.norm_type])

    plotSED_OnAxis(ax, castroData1, TS_thresh, errSigma,
                   colorLim='blue', colorPoint='blue')
    plotSED_OnAxis(ax, castroData2, TS_thresh, errSigma,
                   colorLim='red', colorPoint='red')

    for spec in specVals:
        ax.loglog(castroData1.refSpec.eref, spec)

    return fig, ax