more filtering

data/xspec/GC06.xcm

@@ -0,0 +1,12 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0
+model  atable{../../data/polarmodel.fits} + powerlaw
+        0.73486          1          0          0          3          3
+    2.27115e-09       0.01          0          0      1e+20      1e+24
+        1.23053       0.01         -3         -2          9         10
+      0.0103601       0.01          0          0      1e+20      1e+24
+bayes off

data/xspec/GC10.xcm

@@ -0,0 +1,12 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0
+model  atable{../../data/polarmodel.fits} + powerlaw
+       0.713934          1          0          0          3          3
+    2.27522e-09       0.01          0          0      1e+20      1e+24
+        1.80029       0.01         -3         -2          9         10
+       0.138919       0.01          0          0      1e+20      1e+24
+bayes off

data/xspec/LON+20.xcm

@@ -0,0 +1,10 @@
+method leven 10 0.01
+abund angr
+xsect vern
+cosmo 70 0 0.73
+xset delta 0.01
+systematic 0
+model  atable{../../data/polarmodel.fits}
+       0.724678          1          0          0          3          3
+    1.03624e-09       0.01          0          0      1e+20      1e+24
+bayes off

ridge/ridge/config.py

@@ -45,6 +45,11 @@ SEM означает standatd error on mean (~RMS/sqrt(n))
 """
 sem_cut=98 
 
+"""
+Измерение GRXE ((data-model)/crab) со значением ошибки выше этого персентиля будет отброшено
+"""
+grxe_err_cut=90 
+
 """ Координаты Краба """
 crab_ra  = 83.6287
 crab_dec = 22.014

scripts/02_grxe_resid.py

@@ -43,8 +43,9 @@ with open(proddir+fn.replace(".fits",".ignored_rev.pkl"), 'rb') as fp:
 with open(proddir+fn.replace(".fits",".crabmodel.pkl"), 'rb') as fp:
     crabmodel, z = pickle.load(fp)
     p = np.poly1d(z)
+    crabmodel_keys = list(crabmodel.keys())
     #print(crabmodel)
-
+#sys.exit()
     
 crab_rev_max = np.max(list(crabmodel.keys()))
 print("Crab is defined untill orbit {}".format(crab_rev_max))
@@ -66,11 +67,13 @@ clean0=[]
 model0=[]
 resid0=[] # residuals in cts/s
 grxe0=[] # mCrab
+grxe_err0=[] # mCrab
 crab0=[] # Crab count rate
 mjd0=[]
 a0=[]
 b0=[]
-err0=[]
+model_err0=[]
+crab_err0=[]
 lon0=[]
 lat0=[]
 base0=[]
@@ -114,11 +117,22 @@ for i, row in df.iterrows():
     a = bgdmodel[orbit]['a']
     b = bgdmodel[orbit]['b']
     c = bgdmodel[orbit]['c']
+
+    # Crab error is defined only for Crab resolutions, so we interpolate between
+    if(orbit in crabmodel_keys):
+        crab_err = crabmodel[orbit]['err']
+    else:
+        left,right = find_nearest(crabmodel_keys, orbit)
+        crab_err   = np.interp(orbit, [left,right], [crabmodel[left]['err'], crabmodel[right]['err']])
+        #print()
+        #print(orbit, left, right)
+        #print(orbit, 'err', crabmodel[left]['err'], crab_err, crabmodel[right]['err'])
+        #sys.exit()
     err = bgdmodel[orbit]['err']
     m = a*row['PHASE']+b
     r1 = bgdmodel[orbit]['r1'] # nearest left orbit used for calibration 
     r2 = bgdmodel[orbit]['r2'] # nearest right orbit used for calibration 
-
+    
     c0.append(c)
     base0.append(abs(orbit - int(np.min([r1,r2]))))
     clean0.append(clean[i])
@@ -126,12 +140,14 @@ for i, row in df.iterrows():
     model0.append(m)
     resid0.append(clean[i]-m)
     grxe0.append(1000*(clean[i]-m)/p(orbit))
+    grxe_err0.append(1000*np.sqrt(err**2 + crab_err**2)/p(orbit))
     crab0.append(p(orbit))
         
     a0.append(a)
     b0.append(b)
     c0.append(c)
-    err0.append(err)
+    model_err0.append(err)
+    crab_err0.append(crab_err)
     phase0.append(row['PHASE'])
     rev0.append(orbit)
     lon0.append(row['LON'])
@@ -155,13 +171,15 @@ coldefs = fits.ColDefs([
     fits.Column(name='PHASE', format='D', unit='', array=[phase0[index] for index in indices]),
     fits.Column(name='CLEAN', format='D', unit='cts/s', array=[clean0[index] for index in indices]),
     fits.Column(name='MODEL', format='D', unit='cts/s', array=[model0[index] for index in indices]),
+    fits.Column(name='MODEL_ERR', format='D', unit='', array=[model_err0[index] for index in indices]),
     fits.Column(name='RESID', format='D', unit='cts/s', array=[resid0[index] for index in indices]),
     fits.Column(name='GRXE', format='D', unit='mCrab', array=[grxe0[index] for index in indices]),
+    fits.Column(name='GRXE_ERR', format='D', unit='mCrab', array=[grxe_err0[index] for index in indices]),
     fits.Column(name='CRAB', format='D', unit='cts/s', array=[crab0[index] for index in indices]),
+    fits.Column(name='CRAB_ERR', format='D', unit='', array=[crab_err0[index] for index in indices]),
     fits.Column(name='A', format='D', unit='', array=[a0[index] for index in indices]),
     fits.Column(name='B', format='D', unit='', array=[b0[index] for index in indices]),
     fits.Column(name='C', format='D', unit='', array=[c0[index] for index in indices]),
-    fits.Column(name='ERR', format='D', unit='', array=[err0[index] for index in indices]),
 ])
 
 fout = fn.replace(".fits",".{}.resid.fits".format(outkey))

scripts/03_grxe_spec.py

@@ -85,6 +85,13 @@ for skey in skeys:
         #plt.show()
 
         grxe = np.array(df['GRXE'])
+        grxe_err = np.array(df['GRXE_ERR'])
+
+        perc = np.percentile(grxe_err, grxe_err_cut, axis=0, keepdims=False)
+        print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(skey,enkey,sem_cut,perc))
+        idx=np.where(grxe_err < perc)
+        grxe=grxe[idx]
+        grxe_err=grxe_err[idx]
 
         filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=10, return_bounds=True)
         filtered_grxe = filtered_data[0]
@@ -99,7 +106,7 @@ for skey in skeys:
 
         #sg_sem*=1.5
         #if(sg_mean<0.0):
-        #    sg_mean=1e-9
+        #    sg_mean=1e-6
         #    sg_sem*=2
 
         ebands0[enkey]=[sg_mean,sg_sem]   
@@ -119,7 +126,7 @@ for skey in skeys:
     fspec="{}{}.spec".format(specdir,skey)
     with open(fspec, 'w') as fp:
         for enkey in ebands0.keys():
-            fp.write("0 {} {:.2f} {:.2f} 0.0\n".format(bands[enkey],ebands0[enkey][0],ebands0[enkey][1]))
+            fp.write("0 {} {:.6f} {:.6f} 0.0\n".format(bands[enkey],ebands0[enkey][0],ebands0[enkey][1]))
         
     subprocess.run(["perl", "do_pha_igr_v3_mCrab.pl", fspec])
     try:

scripts/03_plot_skyreg.py

@@ -71,6 +71,7 @@ plt.show()
 
 rev = np.array(df['REV'])
 grxe = np.array(df['GRXE'])
+grxe_err = np.array(df['GRXE_ERR'])
 
 
 
@@ -86,6 +87,35 @@ sstr = '%-14s mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f'
 n, (smin, smax), sm, sv, ss, sk = describe(grxe)
 print(sstr % ('sample:', sm, sv, ss, sk))
     
+# Calculate the percentiles across the x dimension
+errmax=np.max(grxe_err)
+perc = np.percentile(grxe_err, grxe_err_cut, axis=0, keepdims=False)
+print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(skey,enkey,sem_cut,perc))
+idx=np.where(grxe_err < perc)
+
+plt.hist(grxe_err, bins=n_bins, range=[0,errmax], color="red")
+rev=rev[idx]
+grxe=grxe[idx]
+grxe_err=grxe_err[idx]
+plt.hist(grxe_err, bins=n_bins, range=[0,errmax], color="grey")
+
+plt.xlabel("GRXE ERROR, mCrab")
+plt.title("Distribution of errors {}".format(skey))
+plt.show()
+
+grxe_sign=np.divide(grxe,grxe_err)
+plt.hist(grxe_sign, bins=n_bins)
+plt.xlabel("GRXE S/N")
+plt.title("Distribution of significance {}".format(skey))
+plt.show()
+
+A=np.sum(grxe/grxe_err**2)
+B=np.sum(1.0/grxe_err**2)
+wgt_mean=A/B
+wgt_mean_err=np.sqrt(1.0/B)
+
+print("Weighted mean: {:.2f}+/-{:.2f} normal mean: {:.2f}".format(wgt_mean,wgt_mean_err,np.mean(grxe)))
+
 #filtered_grxe = sigma_clip(grxe, sigma=sigma, maxiters=10)
 filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=10, return_bounds=True)
 filtered_grxe = filtered_data[0]