generated from erosita/uds
116 lines
3.7 KiB
Python
Executable File
116 lines
3.7 KiB
Python
Executable File
#!/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, os
|
|
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 *
|
|
|
|
inkey="ALL"
|
|
|
|
n_bins = 80
|
|
sigma=3
|
|
plotme=False
|
|
|
|
ebands0={#'E02':[0.0,0.0],
|
|
'E03':[0.0,0.0],
|
|
'E04':[0.0,0.0],
|
|
'E05':[0.0,0.0],
|
|
'E06':[0.0,0.0],
|
|
'E07':[0.0,0.0],
|
|
'E08':[0.0,0.0],
|
|
'E09':[0.0,0.0],
|
|
'E10':[0.0,0.0],
|
|
'E11':[0.0,0.0],
|
|
'E12':[0.0,0.0],}
|
|
|
|
#skey='Geminga'
|
|
skey = sys.argv[1]
|
|
|
|
if not skey in skyreg.keys():
|
|
print("{} not found in {}".format(skey,list(skyreg.keys())))
|
|
sys.exit()
|
|
|
|
for enkey in ebands0.keys():
|
|
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))
|
|
|
|
#print("Number of ScWs: {}".format(df.shape[0]))
|
|
if not (df.shape[0]>0):
|
|
continue
|
|
#plt.scatter(df['LON'],df['LAT'])
|
|
#plt.show()
|
|
|
|
grxe = np.array(df['GRXE'])
|
|
|
|
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-9
|
|
sg_sem*=2
|
|
|
|
ebands0[enkey]=[sg_mean,sg_sem]
|
|
|
|
if(plotme):
|
|
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()
|
|
|
|
|
|
if not os.path.exists(specdir):
|
|
os.makedirs(specdir)
|
|
|
|
with open("{}{}.spec".format(specdir,skey), '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]))
|