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

sim

Browse Source
This commit is contained in:
Roman Krivonos
2024-07-08 18:57:46 +03:00
parent d7ddc5dc43
commit ac0cea1a3a
941 changed files with 566503 additions and 93421 deletions
Download Patch File Download Diff File Expand all files Collapse all files

79
scripts/03_grxe_spec.py
View File

@@ -7,6 +7,8 @@ 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
@@ -35,9 +37,11 @@ from ridge.config import *
inkey="ALL"
n_bins = 80
sigma=3
plotme=False
nsim=20000
simfrac=10
"""
ebands0={#'E02':[0.0,0.0],
@@ -123,6 +127,9 @@ ebands_sim={'B01':[],
'B21':[],
}
mcrab=u.def_unit('mCrab')
ctss=u.def_unit('cts/s')
u.add_enabled_units([mcrab,ctss])
#skey='Geminga'
if len(sys.argv) > 1:
@@ -136,7 +143,7 @@ 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(ignored_rev)
ign=ignored_rev.tolist()
@@ -148,7 +155,6 @@ else:
sys.exit()
"""
nsim=1000
for skey in skeys:
if not skey in skyreg.keys():
@@ -167,10 +173,10 @@ for skey in skeys:
#with fits.open(proddir+fn) as data:
# df = pd.DataFrame(data[1].data)
dat = Table.read(proddir+fn)
dat = Table.read(proddir+fn, unit_parse_strict='silent')
df = dat.to_pandas()
print(df.columns)
#print(df.columns)
#sys.exit()
#df = df.query("REV == @ign")
@@ -181,26 +187,28 @@ for skey in skeys:
skyreg[skey]['lat'] - skyreg[skey]['wlat']/2,
skyreg[skey]['lat'] + skyreg[skey]['wlat']/2)
print(query)
#sys.exit()
df = df.query(query)
print("{}, {}: {} N={}".format(skey, enkey, query, df.shape[0]))
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))
with open("{}fits/{}.{}.livetime".format(specdir,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
#plt.scatter(df['LON'],df['LAT'])
#plt.show()
print("*** Data frame size {} ***".format(df.size))
sg_mean,sg_sem = get_spec(df, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey)
#print("*** {} {} Data Frame size {} ***".format(skey, enkey, df.size))
sg_mean,sg_sem = get_spec(df, sigma=2, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=False)
ebands0[enkey]=[sg_mean,sg_sem]
nsel = int(df.shape[0]/10)
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)
@@ -208,24 +216,14 @@ for skey in skeys:
#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])
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()
###
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]))
flux=ebands0[enkey][0]
err=ebands0[enkey][1]
print("[DATA] {}: {} {:.6f} {:.6f}".format(skey,enkey,flux,err))
fp.write("0 {} {:.6f} {:.6f} 0.0\n".format(bands[enkey],flux,err))
subprocess.run(["perl", "do_pha_igr_v3_mCrab.pl", fspec])
###
@@ -236,35 +234,30 @@ for skey in skeys:
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)
fp.write("0 {} {:.6f} {:.6f} 0.0\n".format(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)
y = norm.pdf(bins, mu, sg)
l = plt.plot(bins, y, 'r--', linewidth=2)
#plot
plt.xlabel('Flux')
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(r'$\mathrm{Histogram\ of\ IQ:}\ \mu=%.3f,\ \sigma=%.3f$' %(mu, sg))
plt.title("{} {:.2f} {:.2f}".format(enkey, mu, sg))
plt.title("[BOOT] {}: {:.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])

1
scripts/do_pha_igr_v3_mCrab.pl
View File

@@ -13,6 +13,7 @@ $mean_rms=0.73;
$inp=$ARGV[0];
print "*** READING $inp ***\n";
$out=$inp;
open INPP,">$inp.dat";
@spec=`cat $inp`;