generated from erosita/uds
minor
This commit is contained in:
Roman Krivonos
@@ -26,9 +26,10 @@ from numpy import arange
|
||||
from ridge.utils import *
|
||||
from ridge.config import *
|
||||
|
||||
plotme=False
|
||||
enkey = sys.argv[1]
|
||||
#outkey = sys.argv[2]
|
||||
outkey = "ALL"
|
||||
outkey = sys.argv[2]
|
||||
#outkey = "ALL"
|
||||
|
||||
fn="detcnts.{}.fits".format(enkey)
|
||||
|
||||
@@ -60,6 +61,7 @@ rev = data.field('rev')
|
||||
mjd = data.field('mjd')
|
||||
clean = data.field('clean')
|
||||
phase = data.field('phase')
|
||||
texp = data.field('exposure')
|
||||
|
||||
obsid0=[]
|
||||
rev0=[]
|
||||
@@ -79,6 +81,7 @@ lon0=[]
|
||||
lat0=[]
|
||||
base0=[]
|
||||
c0=[]
|
||||
texp0=[]
|
||||
|
||||
d = fits.getdata(datadir+fn)
|
||||
df = pd.DataFrame(np.array(d).byteswap().newbyteorder())
|
||||
@@ -139,6 +142,7 @@ for i, row in df.iterrows():
|
||||
base0.append(abs(orbit - int(np.min([r1,r2]))))
|
||||
clean0.append(clean[i])
|
||||
mjd0.append(mjd[i])
|
||||
texp0.append(texp[i])
|
||||
model0.append(m)
|
||||
resid0.append(clean[i]-m)
|
||||
grxe0.append(1000*(clean[i]-m)/p(orbit))
|
||||
@@ -156,6 +160,57 @@ for i, row in df.iterrows():
|
||||
lat0.append(row['LAT'])
|
||||
|
||||
|
||||
rev_min=np.min(rev0)
|
||||
rev_max=np.max(rev0)
|
||||
orbits=np.array(rev0)
|
||||
resid_arr=np.array(resid0)
|
||||
distr_val=[]
|
||||
distr_rev=[]
|
||||
for r in range(rev_min,rev_max):
|
||||
ind = np.nonzero(orbits == r)
|
||||
if not (len(resid_arr[ind])>10):
|
||||
continue
|
||||
#print(r,np.mean(resid_arr[ind]), len(resid_arr[ind]))
|
||||
distr_val.append(np.mean(resid_arr[ind]))
|
||||
distr_rev.append(r)
|
||||
|
||||
dval=np.array(distr_val)
|
||||
drev=np.array(distr_rev)
|
||||
|
||||
sigma=2
|
||||
n_bins=20
|
||||
filtered_data = sigma_clip(dval, sigma=sigma, maxiters=10, return_bounds=True)
|
||||
filtered_arr=filtered_data[0]
|
||||
filtered_min=filtered_data[1]
|
||||
filtered_max=filtered_data[2]
|
||||
|
||||
|
||||
print("length orig: {} taken: {} filtered: {}".format(len(dval),len(dval[filtered_arr.mask==False]),len(dval[filtered_arr.mask==True])))
|
||||
|
||||
sg_mean, sg_med, sg_std = sigma_clipped_stats(distr_val, sigma=sigma, maxiters=10)
|
||||
sg_sem = sem(dval[filtered_arr.mask==False])
|
||||
print("Sigma clipping: mean {:.2f} med {:.2f} std {:.2f} ".format(sg_mean, sg_med, sg_std))
|
||||
|
||||
|
||||
print(len(drev[filtered_arr.mask==True]))
|
||||
|
||||
if(plotme):
|
||||
k=1.2
|
||||
plt.hist(dval, bins=n_bins, range=[filtered_min*k, filtered_max*k])
|
||||
plt.hist(dval[filtered_arr.mask], bins=n_bins, range=[filtered_min*k, filtered_max*k])
|
||||
plt.axvline(sg_mean, color="black")
|
||||
#plt.axvline(sg_mean+sg_sem, color="black", linestyle="dashed")
|
||||
#plt.axvline(sg_mean-sg_sem, color="black", linestyle="dashed")
|
||||
plt.axvline(sg_mean+sg_std, color="blue", linestyle="dashed")
|
||||
plt.axvline(sg_mean-sg_std, color="blue", linestyle="dashed")
|
||||
plt.xlabel("Resid")
|
||||
plt.show()
|
||||
|
||||
with open(proddir+fn.replace(".fits",".ignored_rev.resid.pkl"), 'wb') as fp:
|
||||
pickle.dump(drev[filtered_arr.mask==True], fp, protocol=pickle.HIGHEST_PROTOCOL)
|
||||
print("Removed REVs:",ignored_rev)
|
||||
|
||||
|
||||
indices = sorted(
|
||||
range(len(mjd0)),
|
||||
key=lambda index: mjd0[index]
|
||||
@@ -170,6 +225,7 @@ coldefs = fits.ColDefs([
|
||||
fits.Column(name='REV', format='J', unit='', array=[rev0[index] for index in indices]),
|
||||
fits.Column(name='BASE', format='J', unit='', array=[base0[index] for index in indices]), # nearest BKG
|
||||
fits.Column(name='MJD', format='D', unit='', array=[mjd0[index] for index in indices]),
|
||||
fits.Column(name='TEXP', format='D', unit='', array=[texp0[index] for index in indices]),
|
||||
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]),
|
||||
|
||||
@@ -1,77 +1,71 @@
|
||||
./02_grxe_resid.py E01 ALL
|
||||
./02_grxe_resid.py E02 ALL
|
||||
./02_grxe_resid.py E03 ALL
|
||||
./02_grxe_resid.py E04 ALL
|
||||
./02_grxe_resid.py E05 ALL
|
||||
./02_grxe_resid.py E06 ALL
|
||||
./02_grxe_resid.py E07 ALL
|
||||
./02_grxe_resid.py E08 ALL
|
||||
./02_grxe_resid.py E09 ALL
|
||||
./02_grxe_resid.py E10 ALL
|
||||
./02_grxe_resid.py E11 ALL
|
||||
./02_grxe_resid.py E12 ALL
|
||||
./02_grxe_resid.py E13 ALL
|
||||
./02_grxe_resid.py E14 ALL
|
||||
./02_grxe_resid.py E15 ALL
|
||||
#./02_grxe_resid.py E01 ALL
|
||||
#./02_grxe_resid.py E02 ALL
|
||||
#./02_grxe_resid.py E03 ALL
|
||||
#./02_grxe_resid.py E04 ALL
|
||||
#./02_grxe_resid.py E05 ALL
|
||||
#./02_grxe_resid.py E06 ALL
|
||||
#./02_grxe_resid.py E07 ALL
|
||||
#./02_grxe_resid.py E08 ALL
|
||||
#./02_grxe_resid.py E09 ALL
|
||||
#./02_grxe_resid.py E10 ALL
|
||||
#./02_grxe_resid.py E11 ALL
|
||||
#./02_grxe_resid.py E12 ALL
|
||||
#./02_grxe_resid.py E13 ALL
|
||||
#./02_grxe_resid.py E14 ALL
|
||||
#./02_grxe_resid.py E15 ALL
|
||||
|
||||
./02_grxe_resid.py A01 ALL
|
||||
./02_grxe_resid.py A02 ALL
|
||||
./02_grxe_resid.py A03 ALL
|
||||
./02_grxe_resid.py A04 ALL
|
||||
./02_grxe_resid.py A05 ALL
|
||||
./02_grxe_resid.py A06 ALL
|
||||
./02_grxe_resid.py A07 ALL
|
||||
./02_grxe_resid.py A08 ALL
|
||||
./02_grxe_resid.py A09 ALL
|
||||
./02_grxe_resid.py A10 ALL
|
||||
./02_grxe_resid.py A11 ALL
|
||||
#./02_grxe_resid.py A01 ALL
|
||||
#./02_grxe_resid.py A02 ALL
|
||||
#./02_grxe_resid.py A03 ALL
|
||||
#./02_grxe_resid.py A04 ALL
|
||||
#./02_grxe_resid.py A05 ALL
|
||||
#./02_grxe_resid.py A06 ALL
|
||||
#./02_grxe_resid.py A07 ALL
|
||||
#./02_grxe_resid.py A08 ALL
|
||||
#./02_grxe_resid.py A09 ALL
|
||||
#./02_grxe_resid.py A10 ALL
|
||||
#./02_grxe_resid.py A11 ALL
|
||||
|
||||
./02_grxe_resid.py B01 ALL
|
||||
./02_grxe_resid.py B02 ALL
|
||||
./02_grxe_resid.py B03 ALL
|
||||
./02_grxe_resid.py B04 ALL
|
||||
./02_grxe_resid.py B05 ALL
|
||||
./02_grxe_resid.py B06 ALL
|
||||
./02_grxe_resid.py B07 ALL
|
||||
./02_grxe_resid.py B08 ALL
|
||||
./02_grxe_resid.py B09 ALL
|
||||
./02_grxe_resid.py B10 ALL
|
||||
./02_grxe_resid.py B11 ALL
|
||||
./02_grxe_resid.py B12 ALL
|
||||
./02_grxe_resid.py B13 ALL
|
||||
./02_grxe_resid.py B14 ALL
|
||||
./02_grxe_resid.py B15 ALL
|
||||
./02_grxe_resid.py B16 ALL
|
||||
./02_grxe_resid.py B17 ALL
|
||||
./02_grxe_resid.py B18 ALL
|
||||
./02_grxe_resid.py B19 ALL
|
||||
./02_grxe_resid.py B20 ALL
|
||||
./02_grxe_resid.py B21 ALL
|
||||
|
||||
|
||||
#./02_grxe_resid.py E01 GAL
|
||||
#./02_grxe_resid.py E02 GAL
|
||||
#./02_grxe_resid.py E03 GAL
|
||||
#./02_grxe_resid.py E04 GAL
|
||||
#./02_grxe_resid.py E05 GAL
|
||||
#./02_grxe_resid.py E06 GAL
|
||||
#./02_grxe_resid.py E07 GAL
|
||||
#./02_grxe_resid.py E08 GAL
|
||||
#./02_grxe_resid.py E09 GAL
|
||||
#./02_grxe_resid.py E10 GAL
|
||||
#./02_grxe_resid.py E11 GAL
|
||||
#./02_grxe_resid.py E12 GAL
|
||||
|
||||
#./02_grxe_resid.py E01 BKG
|
||||
#./02_grxe_resid.py E02 BKG
|
||||
#./02_grxe_resid.py E03 BKG
|
||||
#./02_grxe_resid.py E04 BKG
|
||||
#./02_grxe_resid.py E05 BKG
|
||||
#./02_grxe_resid.py E06 BKG
|
||||
#./02_grxe_resid.py E07 BKG
|
||||
#./02_grxe_resid.py E08 BKG
|
||||
#./02_grxe_resid.py E09 BKG
|
||||
#./02_grxe_resid.py E10 BKG
|
||||
#./02_grxe_resid.py E11 BKG
|
||||
#./02_grxe_resid.py E12 BKG
|
||||
#./02_grxe_resid.py B01 ALL
|
||||
#./02_grxe_resid.py B02 ALL
|
||||
#./02_grxe_resid.py B03 ALL
|
||||
#./02_grxe_resid.py B04 ALL
|
||||
#./02_grxe_resid.py B05 ALL
|
||||
#./02_grxe_resid.py B06 ALL
|
||||
#./02_grxe_resid.py B07 ALL
|
||||
#./02_grxe_resid.py B08 ALL
|
||||
#./02_grxe_resid.py B09 ALL
|
||||
#./02_grxe_resid.py B10 ALL
|
||||
#./02_grxe_resid.py B11 ALL
|
||||
#./02_grxe_resid.py B12 ALL
|
||||
#./02_grxe_resid.py B13 ALL
|
||||
#./02_grxe_resid.py B14 ALL
|
||||
#./02_grxe_resid.py B15 ALL
|
||||
#./02_grxe_resid.py B16 ALL
|
||||
#./02_grxe_resid.py B17 ALL
|
||||
#./02_grxe_resid.py B18 ALL
|
||||
#./02_grxe_resid.py B19 ALL
|
||||
#./02_grxe_resid.py B20 ALL
|
||||
#./02_grxe_resid.py B21 ALL
|
||||
|
||||
./02_grxe_resid.py B01 BKG
|
||||
./02_grxe_resid.py B02 BKG
|
||||
./02_grxe_resid.py B03 BKG
|
||||
./02_grxe_resid.py B04 BKG
|
||||
./02_grxe_resid.py B05 BKG
|
||||
./02_grxe_resid.py B06 BKG
|
||||
./02_grxe_resid.py B07 BKG
|
||||
./02_grxe_resid.py B08 BKG
|
||||
./02_grxe_resid.py B09 BKG
|
||||
./02_grxe_resid.py B10 BKG
|
||||
./02_grxe_resid.py B11 BKG
|
||||
./02_grxe_resid.py B12 BKG
|
||||
./02_grxe_resid.py B13 BKG
|
||||
./02_grxe_resid.py B14 BKG
|
||||
./02_grxe_resid.py B15 BKG
|
||||
./02_grxe_resid.py B16 BKG
|
||||
./02_grxe_resid.py B17 BKG
|
||||
./02_grxe_resid.py B18 BKG
|
||||
./02_grxe_resid.py B19 BKG
|
||||
./02_grxe_resid.py B20 BKG
|
||||
./02_grxe_resid.py B21 BKG
|
||||
|
||||
@@ -25,7 +25,6 @@ from astropy.stats import sigma_clipped_stats
|
||||
from scipy.stats import norm
|
||||
from scipy.stats import describe
|
||||
from scipy.stats import sem
|
||||
|
||||
import subprocess
|
||||
|
||||
from numpy import absolute
|
||||
@@ -77,6 +76,53 @@ ebands0={'B01':[0.0,0.0],
|
||||
'B21':[0.0,0.0],
|
||||
}
|
||||
|
||||
"""
|
||||
ebands_sim={'B01':[[],[]],
|
||||
'B02':[[],[]],
|
||||
'B03':[[],[]],
|
||||
'B04':[[],[]],
|
||||
'B05':[[],[]],
|
||||
'B06':[[],[]],
|
||||
'B07':[[],[]],
|
||||
'B08':[[],[]],
|
||||
'B09':[[],[]],
|
||||
'B10':[[],[]],
|
||||
'B11':[[],[]],
|
||||
'B12':[[],[]],
|
||||
'B13':[[],[]],
|
||||
'B14':[[],[]],
|
||||
'B15':[[],[]],
|
||||
'B16':[[],[]],
|
||||
'B17':[[],[]],
|
||||
'B18':[[],[]],
|
||||
'B19':[[],[]],
|
||||
'B20':[[],[]],
|
||||
'B21':[[],[]],
|
||||
}
|
||||
"""
|
||||
ebands_sim={'B01':[],
|
||||
'B02':[],
|
||||
'B03':[],
|
||||
'B04':[],
|
||||
'B05':[],
|
||||
'B06':[],
|
||||
'B07':[],
|
||||
'B08':[],
|
||||
'B09':[],
|
||||
'B10':[],
|
||||
'B11':[],
|
||||
'B12':[],
|
||||
'B13':[],
|
||||
'B14':[],
|
||||
'B15':[],
|
||||
'B16':[],
|
||||
'B17':[],
|
||||
'B18':[],
|
||||
'B19':[],
|
||||
'B20':[],
|
||||
'B21':[],
|
||||
}
|
||||
|
||||
|
||||
#skey='Geminga'
|
||||
if len(sys.argv) > 1:
|
||||
@@ -88,55 +134,61 @@ else:
|
||||
if not os.path.exists(specdir):
|
||||
os.makedirs(specdir)
|
||||
|
||||
with open(proddir+'detcnts.B21.ignored_rev.resid.pkl', 'rb') as fp:
|
||||
ignored_rev = pickle.load(fp)
|
||||
print(ignored_rev)
|
||||
|
||||
ign=ignored_rev.tolist()
|
||||
sys.exit()
|
||||
nsim=1000
|
||||
|
||||
for skey in skeys:
|
||||
if not skey in skyreg.keys():
|
||||
print("{} not found in {}".format(skey,list(skyreg.keys())))
|
||||
sys.exit()
|
||||
for enkey in ebands0.keys():
|
||||
bkg_fn="detcnts.{}.BKG.resid.fits".format(enkey,inkey)
|
||||
syserr, bkg_sem = get_syserror(proddir+bkg_fn)
|
||||
|
||||
fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
|
||||
d = fits.getdata(proddir+fn)
|
||||
df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
|
||||
#print(df.columns)
|
||||
|
||||
df = df.query('LON > {} & LON < {} & LAT > {} & LAT < {}'.format(skyreg[skey]['lon'] - skyreg[skey]['wlon']/2,
|
||||
skyreg[skey]['lon'] + skyreg[skey]['wlon']/2,
|
||||
skyreg[skey]['lat'] - skyreg[skey]['wlat']/2,
|
||||
skyreg[skey]['lat'] + skyreg[skey]['wlat']/2))
|
||||
#df = df.query("REV == @ign")
|
||||
|
||||
#print("Number of ScWs: {}".format(df.shape[0]))
|
||||
|
||||
df = df.query("LON > {} & LON < {} & LAT > {} & LAT < {} & REV != @ign".format(
|
||||
skyreg[skey]['lon'] - skyreg[skey]['wlon']/2,
|
||||
skyreg[skey]['lon'] + skyreg[skey]['wlon']/2,
|
||||
skyreg[skey]['lat'] - skyreg[skey]['wlat']/2,
|
||||
skyreg[skey]['lat'] + skyreg[skey]['wlat']/2)
|
||||
)
|
||||
|
||||
|
||||
t = Table.from_pandas(df)
|
||||
t.write("{}fits/{}.{}.fits".format(specdir,skey,enkey),overwrite=True)
|
||||
|
||||
texp = np.array(df['TEXP'])
|
||||
print("{} Number of ScWs: {}, {:.1f} Ms".format(skey,df.shape[0],np.sum(texp)/1e6))
|
||||
if not (df.shape[0]>0):
|
||||
continue
|
||||
|
||||
#plt.scatter(df['LON'],df['LAT'])
|
||||
#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,grxe_err_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]
|
||||
filtered_min = filtered_data[1]
|
||||
filtered_max = filtered_data[2]
|
||||
|
||||
#print("length orig: {} taken: {} filtered: {}".format(len(grxe),len(grxe[filtered_grxe.mask==False]),len(grxe[filtered_grxe.mask==True])))
|
||||
|
||||
sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=10)
|
||||
sg_sem = sem(grxe[filtered_grxe.mask==False])
|
||||
print("{}: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f} N={}".format(enkey, sg_mean, sg_med, sg_std, sg_sem, len(grxe[filtered_grxe.mask==False])))
|
||||
|
||||
#sg_sem*=1.5
|
||||
if(sg_mean<0.0):
|
||||
sg_mean=1e-6
|
||||
#sg_sem*=2
|
||||
|
||||
print("*** Data frame size {} ***".format(df.size))
|
||||
sg_mean,sg_sem = get_spec(df, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey)
|
||||
ebands0[enkey]=[sg_mean,sg_sem]
|
||||
|
||||
nsel = int(df.shape[0]/10)
|
||||
for n in range(nsim):
|
||||
df0=df.sample(nsel)
|
||||
sg_mean,sg_sem = get_spec(df0, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey)
|
||||
ebands_sim[enkey].append(sg_mean)
|
||||
#ebands_sim[enkey][1].append(sg_sem)
|
||||
|
||||
|
||||
|
||||
if(plotme):
|
||||
k=1.2
|
||||
plt.hist(grxe, bins=n_bins, range=[filtered_min*k, filtered_max*k])
|
||||
@@ -149,12 +201,53 @@ for skey in skeys:
|
||||
plt.xlabel("mCrab")
|
||||
plt.show()
|
||||
|
||||
###
|
||||
fspec="{}{}.spec".format(specdir,skey)
|
||||
with open(fspec, 'w') as fp:
|
||||
for enkey in ebands0.keys():
|
||||
fp.write("0 {} {:.6f} {:.6f} 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])
|
||||
|
||||
###
|
||||
|
||||
|
||||
|
||||
fspec="{}{}.sim.spec".format(specdir,skey)
|
||||
with open(fspec, 'w') as fp:
|
||||
for enkey in ebands_sim.keys():
|
||||
data=ebands_sim[enkey]
|
||||
#print(type(data))
|
||||
|
||||
(mu, sg) = norm.fit(data)
|
||||
#n, bins, patches = plt.hist(data, 60, density=True, facecolor='green', alpha=0.75)
|
||||
|
||||
if(plotme):
|
||||
|
||||
# add a 'best fit' line
|
||||
y = norm.pdf( bins, mu, sg)
|
||||
l = plt.plot(bins, y, 'r--', linewidth=2)
|
||||
|
||||
#plot
|
||||
plt.xlabel('Flux')
|
||||
plt.ylabel('Probability')
|
||||
#plt.title(r'$\mathrm{Histogram\ of\ IQ:}\ \mu=%.3f,\ \sigma=%.3f$' %(mu, sg))
|
||||
plt.title("{} {:.2f} {:.2f}".format(enkey, mu, sg))
|
||||
plt.grid(True)
|
||||
plt.show()
|
||||
|
||||
print(mu,sg)
|
||||
|
||||
filtered_data = sigma_clip(data, sigma=sigma, maxiters=10, return_bounds=True)
|
||||
filtered_arr=filtered_data[0]
|
||||
filtered_min=filtered_data[1]
|
||||
filtered_max=filtered_data[2]
|
||||
sg_mean, sg_med, sg_std = sigma_clipped_stats(data, sigma=sigma, maxiters=10)
|
||||
sg_sem = sem(data)
|
||||
|
||||
fp.write("0 {} {:.6f} {:.6f} 0.0\n".format(bands[enkey],sg_mean,sg_std))
|
||||
subprocess.run(["perl", "do_pha_igr_v3_mCrab.pl", fspec])
|
||||
|
||||
|
||||
try:
|
||||
for remfile in ["cols","cols1","cols2","header",]:
|
||||
os.remove(remfile)
|
||||
|
||||
98
scripts/03_grxe_syserr.py
Executable file
98
scripts/03_grxe_syserr.py
Executable file
@@ -0,0 +1,98 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
__author__ = "Roman Krivonos"
|
||||
__copyright__ = "Space Research Institute (IKI)"
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from astropy.io import fits
|
||||
from astropy.table import Table, Column
|
||||
import matplotlib.pyplot as plt
|
||||
import math, sys
|
||||
import pickle
|
||||
|
||||
from sklearn.linear_model import LinearRegression
|
||||
from sklearn.linear_model import HuberRegressor
|
||||
from sklearn.linear_model import RANSACRegressor
|
||||
from sklearn.linear_model import TheilSenRegressor
|
||||
from sklearn.model_selection import cross_val_score
|
||||
from sklearn.model_selection import RepeatedKFold
|
||||
|
||||
#from statsmodels.robust.scale import huber
|
||||
from astropy.stats import sigma_clip
|
||||
from astropy.stats import sigma_clipped_stats
|
||||
|
||||
from scipy.stats import norm
|
||||
from scipy.stats import describe
|
||||
from scipy.stats import sem
|
||||
|
||||
|
||||
from numpy import absolute
|
||||
from numpy import arange
|
||||
|
||||
from ridge.utils import *
|
||||
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("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'])
|
||||
|
||||
grxe = (clean-model)/clean*100
|
||||
|
||||
|
||||
mean = np.mean(grxe)
|
||||
std = np.std(grxe)
|
||||
print("\n*** Unfiltered:")
|
||||
print("mean {:.2f} std {:.2f}".format(mean,std))
|
||||
print("min {:.2f}".format(grxe.min()))
|
||||
print("max {:.2f}".format(grxe.max()))
|
||||
print("mean {:.2f}".format(grxe.mean()))
|
||||
print("median {:.2f}".format(np.median(grxe)))
|
||||
print("var {:.2f}".format(grxe.var()))
|
||||
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))
|
||||
|
||||
print("\n***")
|
||||
|
||||
filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=10, return_bounds=True)
|
||||
filtered_grxe=filtered_data[0]
|
||||
filtered_min=filtered_data[1]
|
||||
filtered_max=filtered_data[2]
|
||||
|
||||
print("length orig: {} taken: {} filtered: {}".format(len(grxe),len(grxe[filtered_grxe.mask==False]),len(grxe[filtered_grxe.mask==True])))
|
||||
|
||||
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])
|
||||
plt.axvline(sg_mean, color="black")
|
||||
plt.axvline(sg_mean+sg_sem, color="black", linestyle="dashed")
|
||||
plt.axvline(sg_mean-sg_sem, color="black", linestyle="dashed")
|
||||
plt.axvline(sg_mean+sg_std, color="blue", linestyle="dashed")
|
||||
plt.axvline(sg_mean-sg_std, color="blue", linestyle="dashed")
|
||||
plt.xlabel("mCrab")
|
||||
plt.show()
|
||||
58
scripts/polyfit.pro
Normal file
58
scripts/polyfit.pro
Normal file
@@ -0,0 +1,58 @@
|
||||
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
|
||||
Reference in New Issue
Block a user