integral/ridge
1
0
Fork 0
generated from erosita/uds
Code Issues Pull Requests Packages Projects Releases Wiki Activity

July

Browse Source
This commit is contained in:
Roman Krivonos
2024-07-12 13:06:10 +03:00
parent ac0cea1a3a
commit 96c41060a6
1221 changed files with 954203 additions and 80565 deletions
Download Patch File Download Diff File Expand all files Collapse all files

80
scripts/00_stats.py Executable file
View File

@@ -0,0 +1,80 @@
#!/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
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 numpy import absolute
from numpy import arange
from astropy.table import Table, Column
from ridge.utils import *
from ridge.config import *
#enkey = sys.argv[1]
enkey="E01"
inkey="ALL"
fn="detcnts.{}.fits".format(enkey)
print("Reading {}".format(datadir+fn))
with fits.open(datadir+fn) as data:
df = pd.DataFrame(data[1].data)
print(df.columns)
df_bkg = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & ( abs(LAT) > {} | abs(LON) > {}) & PHASE > {} & PHASE < {}'.format(revmin,revmax,bmax,lmax,phmin,phmax))
df_gal = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & abs(LAT) < {} & abs(LON) < {} & PHASE > {} & PHASE < {}'.format(revmin,revmax,bmax,lmax,phmin,phmax))
df_tot = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & PHASE > {} & PHASE < {}'.format(revmin,revmax,phmin,phmax))
print("\n Initial data set")
print(" GAL {} ScWs, {:.1f} Ms".format(df_gal.shape[0], np.sum(df_gal['EXPOSURE'])/1e6))
print(" BKG {} ScWs, {:.1f} Ms".format(df_bkg.shape[0], np.sum(df_bkg['EXPOSURE'])/1e6))
print("Total {} ScWs, {:.1f} Ms".format(df_tot.shape[0], np.sum(df_tot['EXPOSURE'])/1e6))
fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
print("Reading {}".format(proddir+fn))
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df = dat.to_pandas()
df_bkg = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & ( abs(LAT) > {} | abs(LON) > {}) & PHASE > {} & PHASE < {}'.format(revmin,revmax,bmax,lmax,phmin,phmax))
df_gal = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & abs(LAT) < {} & abs(LON) < {} & PHASE > {} & PHASE < {}'.format(revmin,revmax,bmax,lmax,phmin,phmax))
df_tot = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & PHASE > {} & PHASE < {}'.format(revmin,revmax,phmin,phmax))
print("\n Residuals data set")
print(" GAL {} ScWs, {:.1f} Ms".format(df_gal.shape[0], np.sum(df_gal['TEXP'])/1e6))
print(" BKG {} ScWs, {:.1f} Ms".format(df_bkg.shape[0], np.sum(df_bkg['TEXP'])/1e6))
print("Total {} ScWs, {:.1f} Ms".format(df_tot.shape[0], np.sum(df_tot['TEXP'])/1e6))
# read ignored orbits
with open(proddir+'detcnts.B21.ignored_rev.resid.pkl', 'rb') as fp:
ignored_rev = pickle.load(fp)
print("{} orbits ignored".format(len(ignored_rev)))
ign=ignored_rev.tolist()
query = "REV != @ign"
df_bkg = df_bkg.query(query)
df_gal = df_gal.query(query)
df_tot = df_tot.query(query)
print("\n After removing ignored orbits:")
print(" GAL {} ScWs, {:.1f} Ms".format(df_gal.shape[0], np.sum(df_gal['TEXP'])/1e6))
print(" BKG {} ScWs, {:.1f} Ms".format(df_bkg.shape[0], np.sum(df_bkg['TEXP'])/1e6))
print("Total {} ScWs, {:.1f} Ms".format(df_tot.shape[0], np.sum(df_tot['TEXP'])/1e6))

2
scripts/01_bgdmodel.py
View File

@@ -49,7 +49,7 @@ nrev=0
bgdmodel={}
ignored_scw=[]
ignored_rev=[] # ignore orbits with pecular slope over phase
ignored_rev=[334,1760] # ignore orbits with pecular slope over phase
if not os.path.exists(proddir):
os.makedirs(proddir)

6
scripts/01_bgdmodel.sh
View File

@@ -1,4 +1,4 @@
#./01_bgdmodel.py E01
./01_bgdmodel.py E01
#./01_bgdmodel.py E02
#./01_bgdmodel.py E03
#./01_bgdmodel.py E04
@@ -10,8 +10,8 @@
#./01_bgdmodel.py E10
#./01_bgdmodel.py E11
#./01_bgdmodel.py E12
#./01_bgdmodel.py E13
#./01_bgdmodel.py E14
./01_bgdmodel.py E13
./01_bgdmodel.py E14
#./01_bgdmodel.py E15
#./01_bgdmodel.py A01

14
scripts/01_crabmodel.py
View File

@@ -41,7 +41,7 @@ enkey = sys.argv[1]
sigma=3
# for these bands, slope is taken from stacked profile
fixed_slope = ['E04','E05','E06','E07','E08','E09','E10','E11','E12','E15',
fixed_slope = ['E04','E05','E06','E07','E08','E09','E10','E11','E12','E13','E14','E15',
'A02','A03','A04','A05','A06','A07','A08','A09','A10','A11',
'B02','B03','B04','B05','B06','B07','B08','B09','B10','B11','B12','B13','B14','B15','B16','B17','B18','B19','B20','B21']
@@ -49,7 +49,7 @@ fixed_slope = ['E04','E05','E06','E07','E08','E09','E10','E11','E12','E15',
free_slope = ['E01', 'E02', 'E03', 'A01','B01']
# for these bands, slope is fixed at constant (or if positive, which is not allowed)
const_slope = ['E10','E11','E12','A10','A11','B18','B19','B20','B21']
const_slope = ['E10','E11','E12','E13','E14','E15','A10','A11','B18','B19','B20','B21']
# for stacked profile, skip orbits>800 for energy channels <30 keV
skip800 = ['E02','E03','A01','B01']
@@ -63,10 +63,6 @@ fn="detcnts.{}.fits".format(enkey)
with open(proddir+fn.replace(".fits",".ignored_scw.pkl"), 'rb') as fp:
ignored_scw = pickle.load(fp)
#d = fits.getdata(datadir+fn)
#df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
#print(df.columns)
with fits.open(datadir+fn) as data:
df = pd.DataFrame(data[1].data)
@@ -115,17 +111,19 @@ for rev in range(revmin,revmax):
if(rev > 800):
continue
#if(rev > 1000):
# continue
df0 = df.query('SRC > 0.0 & REV == {} & PHASE > {} & PHASE < {} & CRAB_SEP < {}'.format(rev,phmin,phmax,crab_sep_max))
nobs=len(df0)
if not (nobs> crab_nmax):
continue
print(rev,nobs)
for n in df0['CRAB_SEP'].values:
totx.append(n)
for n in df0['SRC'].values:
toty.append(n)
#if(np.min(df0['SRC'])<1e-3):
# print(rev)
# sys.exit()

6
scripts/01_crabmodel.sh
View File

@@ -1,4 +1,4 @@
#./01_crabmodel.py E01
./01_crabmodel.py E01
#./01_crabmodel.py E02
#./01_crabmodel.py E03
#./01_crabmodel.py E04
@@ -10,8 +10,8 @@
#./01_crabmodel.py E10
#./01_crabmodel.py E11
#./01_crabmodel.py E12
#./01_crabmodel.py E13
#./01_crabmodel.py E14
./01_crabmodel.py E13
./01_crabmodel.py E14
#./01_crabmodel.py E15
#./01_crabmodel.py A01

73
scripts/01_crabmodel_plot_poly.py Executable file
View File

@@ -0,0 +1,73 @@
#!/usr/bin/env python
__author__ = "Roman Krivonos"
__copyright__ = "Space Research Institute (IKI)"
from astropy.table import Table, Column
import matplotlib.pyplot as plt
import sys
from ridge.utils import *
from ridge.config import *
scale = 1e-3
fn="detcnts.E01.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df1 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.E14.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df2 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.E13.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df3 = dat.to_pandas().sort_values(by=['REV'])
fig, (ax1, ax2, ax3) = plt.subplots(3, sharex=True, figsize=(9, 7), dpi=100)
#fig.suptitle('Vertically stacked subplots')
#plt.figure(figsize=(8, 6), dpi=80)
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(1)
ax2.spines[axis].set_linewidth(1)
ax3.spines[axis].set_linewidth(1)
ax1.tick_params(axis="both", width=1, labelsize=14)
ax2.tick_params(axis="both", width=1, labelsize=14)
ax3.tick_params(axis="both", width=1, labelsize=14)
ax1.set_title("25-60 keV")
ax2.set_title("60-80 keV")
ax3.set_title("80-200 keV")
ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax3.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax1.errorbar(df1['REV'], df1['B_EST']/scale, yerr=df1['ERR']/scale, fmt='o' )
ax1.plot(df1['REV'], df1['B_POLY']/scale, color='r', linewidth=2)
ax1.grid(visible=True)
ax2.errorbar(df2['REV'], df2['B_EST']/scale, yerr=df2['ERR']/scale, fmt='o' )
ax2.plot(df2['REV'], df2['B_POLY']/scale, color='r', linewidth=2)
ax2.grid(visible=True)
ax3.errorbar(df3['REV'], df3['B_EST']/scale, yerr=df3['ERR']/scale, fmt='o' )
ax3.plot(df3['REV'], df3['B_POLY']/scale, color='r', linewidth=2)
ax3.grid(visible=True)
#x = np.arange(-10, 10, 0.001)
#plot normal distribution with mean 0 and standard deviation 1
#plt.plot(x, norm.pdf(x, 0, 1), color='red', linewidth=2)
plt.xlabel('Revolution',fontsize=14, fontweight='normal')
ax2.set_ylabel('Count rate, x10$^{-3}$ cts s$^{-1}$ pix$^{-1}$',fontsize=14, fontweight='normal')
#plt.xscale('linear')
#plt.yscale('linear')
plt.savefig(proddir+'crabmodel_poly.png', bbox_inches='tight')
plt.close(fig)

99
scripts/01_crabmodel_plot_sys.py Executable file
View File

@@ -0,0 +1,99 @@
#!/usr/bin/env python
__author__ = "Roman Krivonos"
__copyright__ = "Space Research Institute (IKI)"
from astropy.table import Table, Column
import matplotlib.pyplot as plt
import sys
from scipy.stats import norm
from ridge.utils import *
from ridge.config import *
scale = 1e-3
nbins=30
fn="detcnts.E01.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df1 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.E14.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df2 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.E13.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df3 = dat.to_pandas().sort_values(by=['REV'])
fig, (ax1, ax2, ax3) = plt.subplots(3, sharex=True, figsize=(9, 7), dpi=100)
#fig.suptitle('Vertically stacked subplots')
#plt.figure(figsize=(8, 6), dpi=80)
for axis in ['top','bottom','left','right']:
ax1.spines[axis].set_linewidth(1)
ax2.spines[axis].set_linewidth(1)
ax3.spines[axis].set_linewidth(1)
ax1.tick_params(axis="both", width=1, labelsize=14)
ax2.tick_params(axis="both", width=1, labelsize=14)
ax3.tick_params(axis="both", width=1, labelsize=14)
ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax3.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
data=(df1['B_EST']-df1['B_POLY'])/df1['B_POLY']*1000
(mu, sg) = norm.fit(data)
print(mu, sg)
ax1.set_title("25-60 keV, $\\sigma=${:.0f} mCrab".format(sg))
n, bins, patches = ax1.hist(data, nbins, density=True, facecolor='green', alpha=0.75)
y = norm.pdf(bins, mu, sg)
l = ax1.plot(bins, y, 'r--', linewidth=2)
#plot
ax1.axvline(mu, color="black", linewidth=2)
ax1.axvline(mu+sg, color="blue", linestyle="dashed", linewidth=2)
ax1.axvline(mu-sg, color="blue", linestyle="dashed", linewidth=2)
ax1.grid(visible=True)
data=(df2['B_EST']-df2['B_POLY'])/df2['B_POLY']*1000
(mu, sg) = norm.fit(data)
print(mu, sg)
ax2.set_title("60-80 keV, $\\sigma=${:.0f} mCrab".format(sg))
n, bins, patches = ax2.hist(data, nbins, density=True, facecolor='green', alpha=0.75)
y = norm.pdf(bins, mu, sg)
l = ax2.plot(bins, y, 'r--', linewidth=2)
#plot
ax2.axvline(mu, color="black", linewidth=2)
ax2.axvline(mu+sg, color="blue", linestyle="dashed", linewidth=2)
ax2.axvline(mu-sg, color="blue", linestyle="dashed", linewidth=2)
ax2.grid(visible=True)
data=(df3['B_EST']-df3['B_POLY'])/df3['B_POLY']*1000
(mu, sg) = norm.fit(data)
print(mu, sg)
ax3.set_title("80-200 keV, $\\sigma=${:.0f} mCrab".format(sg))
n, bins, patches = ax3.hist(data, nbins, density=True, facecolor='green', alpha=0.75)
# add a 'best fit' line
y = norm.pdf(bins, mu, sg)
l = ax3.plot(bins, y, 'r--', linewidth=2)
#plot
ax3.axvline(mu, color="black", linewidth=2)
ax3.axvline(mu+sg, color="blue", linestyle="dashed", linewidth=2)
ax3.axvline(mu-sg, color="blue", linestyle="dashed", linewidth=2)
ax3.grid(visible=True)
#x = np.arange(-10, 10, 0.001)
#plot normal distribution with mean 0 and standard deviation 1
#plt.plot(x, norm.pdf(x, 0, 1), color='red', linewidth=2)
plt.xlabel('Flux, mCrab',fontsize=14, fontweight='normal')
#ax2.set_ylabel('No, x10$^{-3}$ cts s$^{-1}$ pix$^{-1}$',fontsize=14, fontweight='normal')
#plt.xscale('linear')
#plt.yscale('linear')
plt.savefig(proddir+'crabmodel_sys.png', bbox_inches='tight')
plt.close(fig)

39
scripts/02_grxe_resid.py
View File

@@ -91,34 +91,43 @@ with fits.open(datadir+fn) as data:
# BKG
if(outkey == 'BKG'):
df = df.query('CLEAN > 0.0 & ( abs(LAT) > {} | abs(LON) > {}) & PHASE > {} & PHASE < {}'.format(bmax,lmax,phmin,phmax))
df = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & ( abs(LAT) > {} | abs(LON) > {}) & PHASE > {} & PHASE < {}'.format(revmin,revmax,bmax,lmax,phmin,phmax))
print("N={}".format(df.shape[0]))
# GAL
if(outkey=='GAL'):
df = df.query('CLEAN > 0.0 & abs(LAT) < {} & abs(LON) < {} & PHASE > {} & PHASE < {}'.format(bmax,lmax,phmin,phmax))
df = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & abs(LAT) < {} & abs(LON) < {} & PHASE > {} & PHASE < {}'.format(revmin,revmax,bmax,lmax,phmin,phmax))
print("N={}".format(df.shape[0]))
# ALL
if(outkey=='ALL'):
df = df.query('CLEAN > 0.0 & PHASE > {} & PHASE < {}'.format(phmin,phmax))
df = df.query('REV >= {} & REV< {} & CLEAN > 0.0 & PHASE > {} & PHASE < {}'.format(revmin,revmax,phmin,phmax))
print("N={}".format(df.shape[0]))
"""
GAL 70226 ScWs, 111.7 Ms
BKG 61214 ScWs, 114.3 Ms
Total 131440 ScWs, 225.9 Ms
"""
for i, row in df.iterrows():
orbit=row['REV']
obsid=row['OBSID']#.decode("UTF-8")
if not (orbit > revmin and orbit < revmax):
print("Skip orbit",orbit,row['OBSID'])
if not (orbit >= revmin and orbit < revmax):
#print("Skip orbit",orbit,row['OBSID'])
continue
if not (orbit < crab_rev_max):
if not (orbit <= crab_rev_max):
print("Skip orbit",orbit,obsid)
continue
if (obsid in ignored_scw):
print("Skip ScW",obsid)
#print("Skip ScW",obsid)
continue
if (orbit in ignored_rev):
print("Skip REV",orbit)
#print("Skip REV",orbit)
continue
a = bgdmodel[orbit]['a']
@@ -162,11 +171,13 @@ for i, row in df.iterrows():
lon0.append(row['LON'])
lat0.append(row['LAT'])
print("N={} ScWs, {:.1f} Ms".format(len(resid0),np.sum(texp0)/1e6))
sigma=3
rev_min=np.min(rev0)
rev_max=np.max(rev0)
orbits=np.array(rev0)
resid_arr=np.array(resid0)
resid_arr=np.array(grxe0)
distr_val=[]
distr_rev=[]
for r in range(rev_min,rev_max):
@@ -179,8 +190,10 @@ for r in range(rev_min,rev_max):
dval=np.array(distr_val)
drev=np.array(distr_rev)
sigma=2
n_bins=20
(mu, sg) = norm.fit(dval)
print(mu,sg)
n_bins=40
filtered_data = sigma_clip(dval, sigma=sigma, maxiters=10, return_bounds=True)
filtered_arr=filtered_data[0]
filtered_min=filtered_data[1]
@@ -206,7 +219,7 @@ if(plotme):
#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.xlabel("Residuals, mCrab")
plt.show()
with open(proddir+fn.replace(".fits",".ignored_rev.resid.pkl"), 'wb') as fp:

6
scripts/02_grxe_resid.sh
View File

@@ -1,4 +1,4 @@
#./02_grxe_resid.py E01 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
@@ -10,8 +10,8 @@
#./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 E13 ALL
./02_grxe_resid.py E14 ALL
#./02_grxe_resid.py E15 ALL
#./02_grxe_resid.py A01 ALL

64
scripts/02_grxe_resid_plot.py Executable file
View File

@@ -0,0 +1,64 @@
#!/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
from astropy import units as u
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
import subprocess
from numpy import absolute
from numpy import arange
from ridge.utils import *
from ridge.config import *
inkey="BKG"
skey="GRXE-BKG"
sigma=3
plotme=False
with open(proddir+'detcnts.B21.ignored_rev.resid.pkl', 'rb') as fp:
ignored_rev = pickle.load(fp)
print("{} orbits ignored".format(len(ignored_rev)))
ign=ignored_rev.tolist()
enkey="E01"
fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df = dat.to_pandas()
print("N={}".format(df.shape[0]))
query = "REV != @ign"
df = df.query(query)
print("{} N={}".format(query, df.shape[0]))
t = Table.from_pandas(df)
t.write("{}/BKG.E01.resid_filtered.fits".format(proddir),overwrite=True)
sg_mean,sg_sem = get_spec(df, sigma=sigma, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, fout="test.fits")

176
scripts/03_grxe_flux.py
View File

@@ -7,8 +7,10 @@ import numpy as np
import pandas as pd
from astropy.io import fits
from astropy.table import Table, Column
from astropy import units as u
import matplotlib.pyplot as plt
import math, sys
import math, sys, os
import pickle
from sklearn.linear_model import LinearRegression
@@ -25,7 +27,7 @@ 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
from numpy import arange
@@ -33,68 +35,120 @@ from numpy import arange
from ridge.utils import *
from ridge.config import *
key="GAL"
enkey = sys.argv[1]
fn="detcnts.{}.{}.resid.fits".format(enkey,key)
inkey="ALL"
d = fits.getdata(proddir+fn)
df=pd.DataFrame(np.array(d).byteswap().newbyteorder())
#print(df.columns)
#key="GC"
#sz=5
#lon0=0.0
#lat0=0.0
key="BKG"
sz=5
lon0=0.0
lat0=0.0
df = df.query('LON > {} & LON < {} & LAT > {} & LAT < {}'.format(lon0-sz,lon0+sz,lat0-sz,lat0+sz))
print("Number of ScWs: {}".format(df.shape[0]))
n_bins = 80
sigma=3
plotme=False
grxe = np.array(df['GRXE'])
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))
ebands0={
'E01':[0.0,0.0], # 25-60 keV
'E14':[0.0,0.0], # 60-80 keV
'E13':[0.0,0.0], # 80-200 keV
}
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 len(sys.argv) > 1:
skeys = [sys.argv[1]]
else:
skeys = list(skyreg.keys())
if not os.path.exists(fluxdir):
os.makedirs(fluxdir)
with open(proddir+'detcnts.B21.ignored_rev.resid.pkl', 'rb') as fp:
ignored_rev = pickle.load(fp)
print("{} orbits ignored".format(len(ignored_rev)))
ign=ignored_rev.tolist()
for skey in skeys:
if not skey in skyreg.keys():
print("{} not found in {}".format(skey,list(skyreg.keys())))
sys.exit()
# generate array for bootstrap
ebands_sim={}
for enkey in ebands0.keys():
ebands_sim[enkey] = []
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)
print("Reading {}".format(proddir+fn))
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df = dat.to_pandas()
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)
df = df.query(query)
print("{}, {}: {} N={}".format(skey, enkey, query, df.shape[0]))
t = Table.from_pandas(df)
t.write("{}fits/{}.{}.fits".format(fluxdir,skey,enkey),overwrite=True)
texp = np.array(df['TEXP'])
with open("{}fits/{}.{}.livetime".format(fluxdir,skey,enkey), 'w') as fp:
fp.write("{} {} ScWs: {} Texp: {:.2f} Ms\n".format(skey,enkey,df.shape[0],np.sum(texp)/1e6))
if not (df.shape[0]>0):
continue
#print("*** {} {} Data Frame size {} ***".format(skey, enkey, df.size))
sg_mean,sg_sem = get_spec(df, sigma=sigma, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=plotme)
ebands0[enkey]=[sg_mean,sg_sem]
nsel = int(df.shape[0]/simfrac)
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)
###
fspec="{}{}.dat".format(fluxdir,skey)
with open(fspec, 'w') as fp:
for enkey in ebands0.keys():
flux=ebands0[enkey][0]
err=ebands0[enkey][1]
print("[DATA] {}: {} {:.6f} {:.6f}".format(skey,enkey,flux,err))
fp.write("{} {} {} {:.6f} {:.6f}\n".format(skey,enkey,bands[enkey],flux,err))
###
fspec="{}{}.sim.dat".format(fluxdir,skey)
with open(fspec, 'w') as fp:
for enkey in ebands_sim.keys():
data=ebands_sim[enkey]
(mu, sg) = norm.fit(data)
fp.write("{} {} {} {:.6f} {:.6f}\n".format(skey,enkey,bands[enkey],mu,sg))
print("[BOOT] {}: {} {:.6f} {:.6f}".format(skey,enkey,mu,sg))
if(plotme):
n, bins, patches = plt.hist(data, 60, density=True, facecolor='green', alpha=0.75)
# add a 'best fit' line
y = norm.pdf(bins, mu, sg)
l = plt.plot(bins, y, 'r--', linewidth=2)
#plot
plt.axvline(mu, color="black")
plt.axvline(ebands0[enkey][0], color="black", linestyle="dashed")
#plt.axvline(mu+sg_sem, color="black", linestyle="dashed")
#plt.axvline(mu-sg_sem, color="black", linestyle="dashed")
plt.axvline(mu+sg, color="blue", linestyle="dashed")
plt.axvline(mu-sg, color="blue", linestyle="dashed")
plt.xlabel('Flux, mCrab')
plt.ylabel('Probability')
plt.title("[BOOT] {}: {:.2f} {:.2f}".format(enkey, mu, sg))
plt.grid(True)
plt.show()

66
scripts/03_grxe_spec.py
View File

@@ -40,8 +40,6 @@ inkey="ALL"
sigma=3
plotme=False
nsim=20000
simfrac=10
"""
ebands0={#'E02':[0.0,0.0],
@@ -80,53 +78,6 @@ 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':[],
}
mcrab=u.def_unit('mCrab')
ctss=u.def_unit('cts/s')
u.add_enabled_units([mcrab,ctss])
@@ -136,6 +87,8 @@ if len(sys.argv) > 1:
skeys = [sys.argv[1]]
else:
skeys = list(skyreg.keys())
# don't plot, if all regions are taken
plotme=False
if not os.path.exists(specdir):
@@ -143,23 +96,22 @@ if not os.path.exists(specdir):
with open(proddir+'detcnts.B21.ignored_rev.resid.pkl', 'rb') as fp:
ignored_rev = pickle.load(fp)
#print(ignored_rev)
print("{} orbits ignored".format(len(ignored_rev)))
ign=ignored_rev.tolist()
"""
if(1091 in ign):
print("Removed")
else:
print("Taken")
sys.exit()
"""
for skey in skeys:
if not skey in skyreg.keys():
print("{} not found in {}".format(skey,list(skyreg.keys())))
sys.exit()
# generate array for bootstrap
ebands_sim={}
for enkey in ebands0.keys():
ebands_sim[enkey] = []
for enkey in ebands0.keys():
#bkg_fn="detcnts.{}.BKG.resid.fits".format(enkey,inkey)
#syserr, bkg_sem = get_syserror(proddir+bkg_fn)

2
scripts/do_pha_igr_v3_mCrab.pl
View File

@@ -9,7 +9,7 @@ $alpha=2.1;
$norm=10;
$mean_rms=0.73;
#$mean_rms=0.73;
$inp=$ARGV[0];

BIN
scripts/test.fits Normal file
View File

Binary file not shown.