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

submit

Browse Source
This commit is contained in:
Roman Krivonos
2024-09-27 19:46:24 +03:00
parent 462133ce2e
commit 0839938594
149 changed files with 157799 additions and 157347 deletions
Download Patch File Download Diff File Expand all files Collapse all files

24
scripts/02_grxe_resid_plot.py
View File

@@ -62,7 +62,7 @@ t = Table.from_pandas(df)
t.write("{}/{}.{}.resid_filtered_rev.fits".format(proddir,inkey,enkey),overwrite=True)
fresid1="{}/{}.{}.resid_filtered_spec.fits".format(proddir,inkey,enkey)
sg_mean,sg_sem = get_spec(df, sigma=sigma, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, gaussfit=True, fout=fresid1)
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"
@@ -79,7 +79,7 @@ t = Table.from_pandas(df)
t.write("{}/{}.{}.resid_filtered_rev.fits".format(proddir,inkey,enkey),overwrite=True)
fresid2="{}/{}.{}.resid_filtered_spec.fits".format(proddir,inkey,enkey)
sg_mean,sg_sem = get_spec(df, sigma=sigma, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, gaussfit=True, fout=fresid2)
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"
@@ -96,7 +96,7 @@ t = Table.from_pandas(df)
t.write("{}/{}.{}.resid_filtered_rev.fits".format(proddir,inkey,enkey),overwrite=True)
fresid3="{}/{}.{}.resid_filtered_spec.fits".format(proddir,inkey,enkey)
sg_mean,sg_sem = get_spec(df, sigma=sigma, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, gaussfit=True, fout=fresid3)
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=fresid3)
###
### Plot light curve
@@ -218,6 +218,22 @@ n, bins, patches = ax3.hist(data, nbins, density=True, facecolor='green', alpha=
# add a 'best fit' line
y = norm.pdf(bins, mu, sg)
l = ax3.plot(bins, y, 'r--', linewidth=2)
area = np.sum(n * np.diff(bins))
xdata = bins[:-1]+np.diff(bins)/2
ydata = n
print("Initial Gaiss fit: mu={:.2f} sigma={:.2f}".format(mu,sg))
y_peak = norm.pdf(mu, mu, sg)
params=[
y_peak*area, # height
mu, # mu
sg, # sigma
0.0, # const 1
0.0, # const 2
]
pfit, perr = fit_leastsq(params, xdata, ydata, const_gaussian_ff)
#plt.plot(bins, const_gaussian_ff(bins, pfit), c='black' )
#plot
ax3.axvline(mu, color="black", linewidth=2)
ax3.axvline(mu+sg, color="blue", linestyle="dashed", linewidth=2)
@@ -232,7 +248,7 @@ 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.yscale('linear')
plt.savefig(proddir+'bkgmodel_histogram.png', bbox_inches='tight')
plt.close(fig)

22
scripts/03_grxe_flux.py
View File

@@ -39,7 +39,7 @@ inkey="ALL"
sigma=3
plotme=True
plotme=False
ebands0={
'E01':[0.0,0.0], # 25-60 keV
@@ -50,6 +50,7 @@ ebands0={
if len(sys.argv) > 1:
skeys = [sys.argv[1]]
#simfrac = int(sys.argv[2])
else:
skeys = list(skyreg.keys())
@@ -62,6 +63,9 @@ with open(ignored_rev_file, 'rb') as fp:
ign=ignored_rev.tolist()
for skey in skeys:
if not skey in skyreg.keys():
print("{} not found in {}".format(skey,list(skyreg.keys())))
@@ -100,15 +104,17 @@ for skey in skeys:
if not (df.shape[0]>0):
continue
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]
sg_mean,sg_sem,skew_val,skew_err = get_spec(df, sigma=4, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, bootstrap=False, gaussfit=True)
ebands0[enkey]=[sg_mean,sg_sem]
nsel = int(df.shape[0]/simfrac)
nsel = int(df.shape[0]*simfrac/100)
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)
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)
###
@@ -122,13 +128,13 @@ for skey in skeys:
###
fspec="{}{}.sim.dat".format(fluxdir,skey)
fspec="{}{}.sim{:02d}.dat".format(fluxdir,skey,simfrac)
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))
fp.write("{:02d} {} {} {} {:.6f} {:.6f}\n".format(simfrac,skey,enkey,bands[enkey],mu,sg))
print("[BOOT] {}: {} {:.6f} {:.6f}".format(skey,enkey,mu,sg))
if(plotme):

23
scripts/03_grxe_flux.sh Normal file
View File

@@ -0,0 +1,23 @@
#./03_grxe_flux.py GC06 1
#./03_grxe_flux.py GC06 2
#./03_grxe_flux.py GC06 3
#./03_grxe_flux.py GC06 4
#./03_grxe_flux.py GC06 5
#./03_grxe_flux.py GC06 10
#./03_grxe_flux.py GC06 15
#./03_grxe_flux.py GC06 20
#./03_grxe_flux.py GC06 25
#./03_grxe_flux.py GC06 30
#./03_grxe_flux.py GC06 35
#./03_grxe_flux.py GC06 40
#./03_grxe_flux.py GC06 45
#./03_grxe_flux.py GC06 50
#./03_grxe_flux.py GC06 55
#./03_grxe_flux.py GC06 60
#./03_grxe_flux.py GC06 65
#./03_grxe_flux.py GC06 70
#./03_grxe_flux.py GC06 75
#./03_grxe_flux.py GC06 80
#./03_grxe_flux.py GC06 85
#./03_grxe_flux.py GC06 90
#./03_grxe_flux.py GC06 95

11
scripts/03_grxe_galprof.py
View File

@@ -108,15 +108,18 @@ for i in range(lon_nbin):
if (df0.shape[0] < nscw_min):
continue
sg_mean,sg_sem = get_spec(df0, sigma=sigma, grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
sg_mean,sg_sem,skew_val,skew_err = get_spec(df0, sigma=sigma, grxe_err_cut=grxe_err_cut, enkey=enkey, plotme=True, nscw_min=nscw_min, gaussfit=True)
nsel = int(df0.shape[0]/simfrac)
nsel = int(df0.shape[0]*simfrac/100)
print("lon {:.2f} ".format(glon[i]),"nsel=",nsel,df0.shape[0])
"""
for n in range(nsim):
df1=df0.sample(nsel)
sg_mean1,sg_sem1 = get_spec(df1, grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
sg_mean1,sg_sem1,skew_val,skew_err = get_spec(df1, grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
mean_sim[dkey].append(sg_mean1)
"""
#print('sg_sem',sg_sem)
mean_map[i] = sg_mean

18
scripts/03_grxe_galprof_plot.py
View File

@@ -50,9 +50,9 @@ ax1.set_xlim(185,-185)
ax2.set_xlim(185,-185)
ax3.set_xlim(185,-185)
ax1.set_ylim(-0.3,1.8)
ax2.set_ylim(-0.3,1.8)
ax3.set_ylim(-0.3,1.8)
ax1.set_ylim(-0.3,1.6)
ax2.set_ylim(-0.3,1.)
ax3.set_ylim(-0.3,1.)
ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax2.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
@@ -68,21 +68,21 @@ ax3.yaxis.set_minor_locator(ticker.MultipleLocator(0.2))
lon=(df1['LON1']+df1['LON2'])/2
ax1.plot(df_cobe[0],df_cobe[2]/3, color='red', linewidth=4)
ax1.errorbar(lon, df1['GRXE_SIM_FLUX']/scale,
yerr=df1['GRXE_SIM_ERROR']/scale,
ax1.errorbar(lon, df1['GRXE_FLUX']/scale,
yerr=df1['GRXE_ERROR']/scale,
xerr=(df1['LON2']-df1['LON1'])/2,
fmt='o', color='black')
ax1.grid(visible=True)
lon=(df2['LON1']+df2['LON2'])/2
ax2.errorbar(lon, df2['GRXE_SIM_FLUX']/scale,
yerr=df2['GRXE_SIM_ERROR']/scale,
ax2.errorbar(lon, df2['GRXE_FLUX']/scale,
yerr=df2['GRXE_ERROR']/scale,
xerr=(df2['LON2']-df2['LON1'])/2,fmt='o' , color='black')
ax2.grid(visible=True)
lon=(df3['LON1']+df3['LON2'])/2
ax3.errorbar(lon, df3['GRXE_SIM_FLUX']/scale,
yerr=df3['GRXE_SIM_ERROR']/scale,
ax3.errorbar(lon, df3['GRXE_FLUX']/scale,
yerr=df3['GRXE_ERROR']/scale,
xerr=(df3['LON2']-df3['LON1'])/2, fmt='o' , color='black')
ax3.grid(visible=True)

2
scripts/03_grxe_map.py
View File

@@ -147,7 +147,7 @@ for i in range(sx):
sg_mean,sg_sem = get_spec(df0, sigma=sigma, grxe_err_cut=grxe_err_cut, enkey=enkey, nscw_min=nscw_min)
nsel = int(df0.shape[0]/simfrac)
nsel = int(df0.shape[0]*simfrac/100)
#print("nsel=",nsel,df0.shape[0])
for n in range(nsim):
df1=df0.sample(nsel)

40
scripts/03_grxe_spec.py
View File

@@ -78,6 +78,29 @@ ebands0={'B01':[0.0,0.0],
'B21':[0.0,0.0],
}
skew0={'B01':[0.0,0.0],
'B02':[0.0,0.0],
'B03':[0.0,0.0],
'B04':[0.0,0.0],
'B05':[0.0,0.0],
'B06':[0.0,0.0],
'B07':[0.0,0.0],
'B08':[0.0,0.0],
'B09':[0.0,0.0],
'B10':[0.0,0.0],
'B11':[0.0,0.0],
'B12':[0.0,0.0],
'B13':[0.0,0.0],
'B14':[0.0,0.0],
'B15':[0.0,0.0],
'B16':[0.0,0.0],
'B17':[0.0,0.0],
'B18':[0.0,0.0],
'B19':[0.0,0.0],
'B20':[0.0,0.0],
'B21':[0.0,0.0],
}
mcrab=u.def_unit('mCrab')
ctss=u.def_unit('cts/s')
u.add_enabled_units([mcrab,ctss])
@@ -158,14 +181,15 @@ for skey in skeys:
#plt.scatter(df['LON'],df['LAT'])
#plt.show()
#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)
sg_mean,sg_sem,skew_val,skew_err = get_spec(df, sigma=3, grxe_err_cut=grxe_err_cut, skey=skey, enkey=enkey, plotme=True, bootstrap=False, gaussfit=True)
ebands0[enkey]=[sg_mean,sg_sem]
skew0[enkey]=[skew_val,skew_err]
nsel = int(df.shape[0]/simfrac)
nsel = int(df.shape[0]*simfrac/100)
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)
#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)
@@ -181,7 +205,13 @@ for skey in skeys:
###
fspec="{}{}.skew".format(specdir,skey)
with open(fspec, 'w') as fp:
fp.write("read serr 4\n")
count=1
for enkey in skew0.keys():
fp.write("{} {} {:.6f} {:.6f}\n".format(count,bands[enkey],skew0[enkey][0],skew0[enkey][1]))
count+=1
fspec="{}{}.sim.spec".format(specdir,skey)
with open(fspec, 'w') as fp: