code: def compare_SED(castroData1, castroData2, ylims, TS_thresh=4.0,
                  errSigma=1.0, specVals=[]):
      
      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
docs: 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