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

clean

Browse Source
This commit is contained in:
Roman Krivonos
2024-11-01 11:40:05 +03:00
parent 770ee50a62
commit bb3b99abd3
1115 changed files with 62 additions and 1376385 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
scripts/01_crabmodel_plot_poly.py
View File

@@ -6,7 +6,7 @@ __copyright__ = "Space Research Institute (IKI)"
from astropy.table import Table, Column
import matplotlib.pyplot as plt
import sys
import sys, os
from ridge.utils import *
from ridge.config import *
@@ -14,15 +14,15 @@ from ridge.config import *
scale = 1e-3
fn="detcnts.E01.crabmodel.fits"
fn="detcnts.A01.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df1 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.E14.crabmodel.fits"
fn="detcnts.A02.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df2 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.E13.crabmodel.fits"
fn="detcnts.A03.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df3 = dat.to_pandas().sort_values(by=['REV'])
@@ -69,5 +69,10 @@ ax2.set_ylabel('Count rate, x10$^{-3}$ cts s$^{-1}$ pix$^{-1}$',fontsize=14, fon
#plt.xscale('linear')
#plt.yscale('linear')
plt.savefig(proddir+'crabmodel_poly.png', bbox_inches='tight')
if not os.path.exists(figdir):
os.makedirs(figdir)
filename=figdir+'crabmodel_poly.png'
plt.savefig(filename, bbox_inches='tight')
plt.close(fig)
print("Result is saved as {}".format(filename))

20
scripts/01_crabmodel_plot_sys.py
View File

@@ -6,7 +6,7 @@ __copyright__ = "Space Research Institute (IKI)"
from astropy.table import Table, Column
import matplotlib.pyplot as plt
import sys
import sys, os
from scipy.stats import norm
from ridge.utils import *
@@ -16,15 +16,15 @@ from ridge.config import *
scale = 1e-3
nbins=30
fn="detcnts.E01.crabmodel.fits"
fn="detcnts.A01.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df1 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.E14.crabmodel.fits"
fn="detcnts.A02.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df2 = dat.to_pandas().sort_values(by=['REV'])
fn="detcnts.E13.crabmodel.fits"
fn="detcnts.A03.crabmodel.fits"
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df3 = dat.to_pandas().sort_values(by=['REV'])
@@ -47,7 +47,7 @@ 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)
#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)
@@ -60,7 +60,7 @@ ax1.grid(visible=True)
data=(df2['B_EST']-df2['B_POLY'])/df2['B_POLY']*1000
(mu, sg) = norm.fit(data)
print(mu, sg)
#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)
@@ -73,7 +73,7 @@ ax2.grid(visible=True)
data=(df3['B_EST']-df3['B_POLY'])/df3['B_POLY']*1000
(mu, sg) = norm.fit(data)
print(mu, sg)
#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
@@ -95,5 +95,9 @@ plt.xlabel('Flux, mCrab',fontsize=14, fontweight='normal')
#plt.xscale('linear')
#plt.yscale('linear')
plt.savefig(proddir+'crabmodel_sys.png', bbox_inches='tight')
if not os.path.exists(figdir):
os.makedirs(figdir)
filename=figdir+'crabmodel_sys.png'
plt.savefig(filename, bbox_inches='tight')
plt.close(fig)
print("Result is saved as {}".format(filename))

23
scripts/02_grxe_resid_plot.py
View File

@@ -48,7 +48,7 @@ with open(ignored_rev_file, 'rb') as fp:
ign=ignored_rev.tolist()
enkey="E01"
enkey="A01"
fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df = dat.to_pandas()
@@ -65,7 +65,7 @@ fresid1="{}/{}.{}.resid_filtered_spec.fits".format(proddir,inkey,enkey)
sg_mean,sg_sem,skew_val,skew_err = get_spec(df, sigma=sigma, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, gaussfit=True, fout=fresid1)
enkey="E14"
enkey="A02"
fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df = dat.to_pandas()
@@ -82,7 +82,7 @@ fresid2="{}/{}.{}.resid_filtered_spec.fits".format(proddir,inkey,enkey)
sg_mean,sg_sem, skew_val, skew_err = get_spec(df, sigma=sigma, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, gaussfit=True, fout=fresid2)
enkey="E13"
enkey="A03"
fn="detcnts.{}.{}.resid.fits".format(enkey,inkey)
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df = dat.to_pandas()
@@ -250,8 +250,10 @@ plt.xlabel('Flux, mCrab',fontsize=14, fontweight='normal')
#plt.xscale('linear')
plt.yscale('linear')
plt.savefig(proddir+'bkgmodel_histogram.png', bbox_inches='tight')
filename=figdir+'bkgmodel_histogram.png'
plt.savefig(filename, bbox_inches='tight')
plt.close(fig)
print("\nResult is saved as {}".format(filename))
###
@@ -277,7 +279,7 @@ ax3.ticklabel_format(style='sci', axis='y', scilimits=(-4,4))
data=(df1['CLEAN']-df1['MODEL'])/df1['CLEAN']*100
(mu, sg) = norm.fit(data)
print(mu, sg)
#print(mu, sg)
txt="25-60 keV, $\\sigma=${:.1f}%".format(sg)
ax1.set_title(txt)
n, bins, patches = ax1.hist(data, nbins, density=True, facecolor='green', alpha=0.75)
@@ -292,7 +294,7 @@ ax1.grid(visible=True)
data=(df2['CLEAN']-df2['MODEL'])/df2['CLEAN']*100
(mu, sg) = norm.fit(data)
print(mu, sg)
#print(mu, sg)
txt="60-80 keV, $\\sigma=${:.1f}%".format(sg)
ax2.set_title(txt)
n, bins, patches = ax2.hist(data, nbins, density=True, facecolor='green', alpha=0.75)
@@ -308,7 +310,7 @@ ax2.grid(visible=True)
data=(df3['CLEAN']-df3['MODEL'])/df3['CLEAN']*100
(mu, sg) = norm.fit(data)
print(mu, sg)
#print(mu, sg)
txt="80-200 keV, $\\sigma=${:.1f}%".format(sg)
ax3.set_title(txt)
@@ -333,5 +335,10 @@ plt.xlabel('Residuals, %',fontsize=14, fontweight='normal')
#plt.xscale('linear')
#plt.yscale('linear')
plt.savefig(proddir+'bkgmodel_systematic.png', bbox_inches='tight')
if not os.path.exists(figdir):
os.makedirs(figdir)
filename=figdir+'bkgmodel_systematic.png'
plt.savefig(filename, bbox_inches='tight')
plt.close(fig)
print("Result is saved as {}".format(filename))

20
scripts/03_grxe_flux.py
View File

@@ -38,13 +38,14 @@ from ridge.config import *
inkey="ALL"
sigma=3
plotme=False
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
'A01':[0.0,0.0], # 25-60 keV
'A02':[0.0,0.0], # 60-80 keV
'A03':[0.0,0.0], # 80-200 keV
}
@@ -57,6 +58,10 @@ else:
if not os.path.exists(fluxdir):
os.makedirs(fluxdir)
fitsdir = "{}fits/".format(fluxdir)
if not os.path.exists(fitsdir):
os.makedirs(fitsdir)
with open(ignored_rev_file, 'rb') as fp:
ignored_rev = pickle.load(fp)
print("{} orbits ignored".format(len(ignored_rev)))
@@ -95,10 +100,10 @@ for skey in skeys:
print("{}, {}: {} N={}".format(skey, enkey, query, df.shape[0]))
t = Table.from_pandas(df)
t.write("{}fits/{}.{}.fits".format(fluxdir,skey,enkey),overwrite=True)
t.write("{}{}.{}.fits".format(fitsdir,skey,enkey),overwrite=True)
texp = np.array(df['TEXP'])
with open("{}fits/{}.{}.livetime".format(fluxdir,skey,enkey), 'w') as fp:
with open("{}{}.{}.livetime".format(fitsdir,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):
@@ -113,10 +118,7 @@ for skey in skeys:
for n in range(nsim):
df0=df.sample(nsel)
sg_mean,sg_sem,skew_val,skew_err = get_spec(df0, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey)
ebands_sim[enkey].append(sg_mean)
ebands_sim[enkey].append(sg_mean)
###
fspec="{}{}.dat".format(fluxdir,skey)
with open(fspec, 'w') as fp:

12
scripts/README.md
View File

@@ -14,10 +14,22 @@ Calibrates IBIS/ISGRI background model
Calibrates Crab detector count rate model.
### 01_crabmodel_plot_poly.py
Plot long-term Crab detector count rate approximated by cubic polynomial function (Fig. B7 in paper).
### 01_crabmodel_plot_sys.py
Plots the normalized distribution of the residuals between the IBIS/ISGRI Crab count rate and the corresponding polynomial fit (Fig. B8 in paper). The distribution is approximated with a Gaussian function.
### 02_grxe_resid.py
Calculates difference between detector count rate and that predicted by background model in mCrab units for `GAL` (Galaxy), `BKG` (Extragalactic), or `ALL` (all-sky) regions.
### 02_grxe_resid_plot.py
Plots normalized distribution of the relative residuals of the background model obtained in three energy bands for BKG region (Figs. A5 and A6 in paper).
### 02_grxe_map.py
Makes the map of the residuals in mCrab units (not covered in the paper).