diff --git a/scripts/03_grxe_syserr.py b/scripts/03_grxe_syserr.py
index 60b7317..1ca3ec3 100755
--- a/scripts/03_grxe_syserr.py
+++ b/scripts/03_grxe_syserr.py
@@ -35,23 +35,15 @@ from ridge.config import *
 
 key="BKG"
 enkey = sys.argv[1]
+
 fn="detcnts.{}.{}.resid.fits".format(enkey,key)
-
-
-d = fits.getdata(proddir+fn)
-df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
-print(fn)
-print(df.columns)
-
-
+print("Reading {}".format(proddir+fn))
+dat = Table.read(proddir+fn, unit_parse_strict='silent')
+df = dat.to_pandas()
 
 print("Number of ScWs: {}".format(df.shape[0]))
 
-
 n_bins = 80
-sigma=3
-
-#grxe = np.array(df['GRXE'])
 
 clean = np.array(df['CLEAN'])
 model = np.array(df['MODEL'])
@@ -85,7 +77,6 @@ sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=10)
 sg_sem = sem(grxe[filtered_grxe.mask==False])
 print("Sigma clipping: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f}".format(sg_mean, sg_med, sg_std, sg_sem))
 
-
 k=1.2
 plt.hist(grxe, bins=n_bins, range=[filtered_min*k, filtered_max*k])
 plt.hist(grxe[filtered_grxe.mask], bins=n_bins, range=[filtered_min*k, filtered_max*k])
diff --git a/scripts/03_plot_bkg.py b/scripts/03_plot_bkg.py
index 84cf738..b6c70f0 100755
--- a/scripts/03_plot_bkg.py
+++ b/scripts/03_plot_bkg.py
@@ -40,7 +40,7 @@ fn='detcnts.{}.{}.resid.fits'.format(enkey,key)
 
 
 n_bins = 80
-sigma=3
+
 
 with fits.open(proddir+fn) as hdul:
     data=hdul[1].data
diff --git a/scripts/03_plot_skyreg.py b/scripts/03_plot_skyreg.py
index 3537a8a..b0419b7 100755
--- a/scripts/03_plot_skyreg.py
+++ b/scripts/03_plot_skyreg.py
@@ -41,10 +41,11 @@ if not skey in skyreg.keys():
     sys.exit()
 
 key = "ALL"
-fn='detcnts.{}.{}.resid.fits'.format(enkey,key)
 
-d = fits.getdata(proddir+fn)
-df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
+fn="detcnts.{}.{}.resid.fits".format(enkey,key)
+print("Reading {}".format(proddir+fn))
+dat = Table.read(proddir+fn, unit_parse_strict='silent')
+df = dat.to_pandas()
 
 
 n_bins = 60
diff --git a/scripts/README.md b/scripts/README.md
index e91f3e5..aec225b 100644
--- a/scripts/README.md
+++ b/scripts/README.md
@@ -44,7 +44,7 @@ GB A01  25.00   60.00 154.886019 3.087652
 GB A02  60.00   80.00 96.472798 9.218161
 GB A03  80.00  200.00 102.207302 13.291950
 ```
-wgere GB means Galactic Bulge region (see Table 1 in the paper).
+where GB means Galactic Bulge region (see Table 1 in the paper).
 
 ### 03_grxe_galprof.py/03_grxe_galprof.sh and 03_grxe_galprof_plot.py
 
@@ -71,3 +71,8 @@ xspec
 @load_gb_cutoffpl.xcm
 fit
 ```
+
+### 03_grxe_syserr.py/03_plot_bkg.py/03_plot_skyreg.py
+
+Prints statistics and plots distribution of residuals for extragalactic data, and for selected sky regions.
+
diff --git a/scripts/polyfit.pro b/scripts/polyfit.pro
deleted file mode 100644
index 00b8c0d..0000000
--- a/scripts/polyfit.pro
+++ /dev/null
@@ -1,58 +0,0 @@
-pro polyfit
-  green = GETCOLOR('green', 100)
-  red = GETCOLOR('red', 100)
-
-
-  enkey="B21"
-  fn = '../products/detcnts.'+enkey+'.ALL.resid.fits'
-  
-  table  = readfits(fn, header, EXTEN_NO=1, /SILENT)
-  
-  
-  ;;tstart = TBGET( header, table, 'TSTART', /NOSCALE)
-  ;;tstop  = TBGET( header, table, 'TSTOP', /NOSCALE)
-  rev   = TBGET( header, table, 'REV', /NOSCALE)
-  data   = TBGET( header, table, 'RESID', /NOSCALE)
-  time = LINDGEN(N_ELEMENTS(data))
-
-
-  rev_min=min(rev)
-  rev_max=max(rev)
-
-  plot_x=[]
-  plot_y=[]
-  
-  for r=rev_min, rev_max do begin
-     index = where(rev eq r, count)
-     if(count lt 10) then continue
-     mn = AVG(data(index))
-     push,plot_x,r
-     push,plot_y,mn
-     print,r,count,mn
-  endfor
-  
-  plot,plot_x,plot_y                ;, yrange=[-0.0001,0.0001]
-
-  return
-
-  
-  ;;COEFF = ROBUST_LINEFIT(time, rate, YFIT, SIG, COEF_SIG)
-  COEFF = ROBUST_POLY_FIT(time,data,7,YFIT,SIG, /DOUBLE)
-
-  print,SIG
-  print,COEFF
-  
-  ;mean_polyfit=AVG(YFIT)
-  ;print,'mean fit',mean_polyfit
-  ;print,'resid',SIG/mean_polyfit*100
-
-  poly = coeff;/mean_polyfit
-  ;;save,FILENAME = enkey+'_poly.sav', poly, mean_polyfit
-  
-  scale=1.2
-  ;;yfit0 = (poly[0] + poly[1]*time + poly[2]*time^2 + poly[3]*time^3)*scale
-  
-  plot,time,data;, yrange=[-0.0001,0.0001]
-  oplot,time,YFIT,color=green
-  ;;oplot,time,YFIT0,color=red
-end