update

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])

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

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

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.
+

scripts/polyfit.pro

@@ -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