Skip to content

Commit

Permalink
Merge pull request #50 from mtakahiro/nirspec_res
Browse files Browse the repository at this point in the history 
bug fixed
  • Loading branch information
@mtakahiro
mtakahiro authored May 12, 2023
2 parents 9af5c35 + 9edd45e commit b973ecb
Show file tree
Hide file tree
Showing 4 changed files with 68 additions and 63 deletions.
26 changes: 13 additions & 13 deletions gsf/fitting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -326,10 +326,10 @@ def update_input(self, inputs, c:float=3e18, Mpc_cm:float=3.08568025e+24, m0set:
self.filts = [x.strip() for x in self.filts.split(',')]
except:
self.filts = []
for column in self.fd_cat.columns:
if column[0] == 'F':
self.filts.append(column[1:])
pass
if not isinstance(self.fd_cat, type(None)):
for column in self.fd_cat.columns:
if column[0] == 'F':
self.filts.append(column[1:])

self.band = {}
for ii in range(len(self.filts)):
Expand Down Expand Up @@ -990,14 +990,13 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, #zliml=0.01, zlim
# ax1.plot(data_model_sort[:,0], data_model_sort[:,1], 'b', linestyle='--', linewidth=0.5, label='')

# Observed data;
# spec;
if self.has_spectrum:
ey_max = 1000
con = (self.dict['ey']<ey_max) & (self.dict['NR']<self.NRbb_lim) & (self.dict['ey']>=0)
ax1.errorbar(self.dict['x'][con], self.dict['fy'][con], yerr=self.dict['ey'][con], color='gray', capsize=0, linewidth=0.5, linestyle='', zorder=4)
ax1.plot(self.dict['x'][con], self.dict['fy'][con], '.r', linestyle='', linewidth=0.5, label='Observed spectrum', zorder=4)
# bb;
if include_photometry:

if include_photometry and self.has_photometry:
con = (self.dict['NR']>=self.NRbb_lim) & (self.dict['ey']>=0)
ax1.errorbar(self.dict['x'][con], self.dict['fy'][con], yerr=self.dict['ey'][con], ms=15, marker='None',
color='orange', capsize=0, linewidth=0.5, ls='None', label='', zorder=4)
Expand Down Expand Up @@ -1147,13 +1146,13 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, #zliml=0.01, zlim

print('\n##############################################################')
print('Input redshift is %.3f per cent agreement.'%((1.-zsigma)*100))
print('Error is %.3f per cent.'%(zzsigma*100))
print('Estimated error is %.3f per cent.'%(zzsigma*100))
print('Input Cz0 is %.3f per cent agreement.'%((1.-C0sigma)*100))
print('Error is %.3f per cent.'%(eC0sigma*100))
print('Estimated error is %.3f per cent.'%(eC0sigma*100))
print('Input Cz1 is %.3f per cent agreement.'%((1.-C1sigma)*100))
print('Error is %.3f per cent.'%(eC1sigma*100))
print('Estimated error is %.3f per cent.'%(eC1sigma*100))
print('Input Cz2 is %.3f per cent agreement.'%((1.-C2sigma)*100))
print('Error is %.3f per cent.'%(eC2sigma*100))
print('Estimated error is %.3f per cent.'%(eC2sigma*100))
print('##############################################################\n')

if fzvis==1:
Expand Down Expand Up @@ -1209,8 +1208,9 @@ def get_zdist(self, f_interact=False, f_ascii=True, return_figure=False):
yy = np.arange(0,np.max(n),1)
xx = yy * 0 + self.z_cz[1]
ax1.plot(xx,yy,linestyle='-',linewidth=1,color='orangered',\
label='$z=%.5f_{-%.5f}^{+%.5f}$\n$C_{z0}=%.3f$\n$C_{z1}=%.3f$\n$C_{z2}=%.3f$'%\
(self.z_cz[1],self.z_cz[1]-self.z_cz[0],self.z_cz[2]-self.z_cz[1], self.Cz0, self.Cz1, self.Cz2))
# label='$z=%.5f_{-%.5f}^{+%.5f}$\n$C_{z0}=%.3f$\n$C_{z1}=%.3f$\n$C_{z2}=%.3f$'%\
label='$z=%.5f_{-%.5f}^{+%.5f}$'%\
(self.z_cz[1],self.z_cz[1]-self.z_cz[0],self.z_cz[2]-self.z_cz[1]))

if f_ascii:
file_ascii_out = self.DIR_OUT + 'zmc_' + self.ID + '.txt'
Expand Down
37 changes: 20 additions & 17 deletions gsf/maketmp_filt.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1260,9 +1260,11 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,
MB.data['spec_fir_obs'] = dict_spec_fir_obs

# BB phot
MB.has_photometry = False
fw = open(file_tmp,'w')
fw_rem = open(file_tmp2, 'w')
for ii in range(len(ltmpbb[0,:])):
MB.has_photometry = True
if SFILT[ii] in SKIPFILT:# data point to be skiped;
fw.write('%d %.5f %.5e %.5e %.1f %s\n'%(ii+ncolbb, ltmpbb[0,ii], 0.0, fbb[ii], FWFILT[ii]/2., SFILT[ii]))
fw_rem.write('%d %.5f %.5e %.5e %.1f %s\n'%(ii+ncolbb, ltmpbb[0,ii], fbb[ii], ebb[ii], FWFILT[ii]/2., SFILT[ii]))
Expand All @@ -1276,23 +1278,24 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,
fw_rem.close()

# register;
dat = ascii.read(file_tmp, format='no_header')
NRbb = dat['col1']
xbb = dat['col2']
fybb = dat['col3']
eybb = dat['col4']
exbb = dat['col5']
dict_bb_obs = {'NR':NRbb, 'x':xbb, 'fy':fybb, 'ey':eybb, 'ex':exbb}
MB.data['bb_obs'] = dict_bb_obs
if len(SKIPFILT)>0:#try:
dat = ascii.read(file_tmp2, format='no_header')
NR_ex = dat['col1']
x_ex = dat['col2']
fy_ex = dat['col3']
ey_ex = dat['col4']
ex_ex = dat['col5']
dict_bb_obs_removed = {'NR':NR_ex, 'x':x_ex, 'fy':fy_ex, 'ey':ey_ex, 'ex':ex_ex}
MB.data['bb_obs_removed'] = dict_bb_obs_removed
if MB.has_photometry:
dat = ascii.read(file_tmp, format='no_header')
NRbb = dat['col1']
xbb = dat['col2']
fybb = dat['col3']
eybb = dat['col4']
exbb = dat['col5']
dict_bb_obs = {'NR':NRbb, 'x':xbb, 'fy':fybb, 'ey':eybb, 'ex':exbb}
MB.data['bb_obs'] = dict_bb_obs
if len(SKIPFILT)>0:#try:
dat = ascii.read(file_tmp2, format='no_header')
NR_ex = dat['col1']
x_ex = dat['col2']
fy_ex = dat['col3']
ey_ex = dat['col4']
ex_ex = dat['col5']
dict_bb_obs_removed = {'NR':NR_ex, 'x':x_ex, 'fy':fy_ex, 'ey':ey_ex, 'ex':ex_ex}
MB.data['bb_obs_removed'] = dict_bb_obs_removed

# Dust; Not sure where this is being used...
fw = open(file_tmp,'w')
Expand Down
66 changes: 34 additions & 32 deletions gsf/plot_sed.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -497,21 +497,27 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=None, f_chind=True, figpdf=Fals
#############
# Main result
#############
conbb_ymax = (xbb>0) & (fybb>0) & (eybb>0) & (fybb/eybb>SNlim)
if len(fybb[conbb_ymax]):
ymax = np.nanmax(fybb[conbb_ymax]*c/np.square(xbb[conbb_ymax])/d_scale) * 1.6
if MB.has_photometry:
conbb_ymax = (xbb>0) & (fybb>0) & (eybb>0) & (fybb/eybb>SNlim)
if len(fybb[conbb_ymax]):
ymax = np.nanmax(fybb[conbb_ymax]*c/np.square(xbb[conbb_ymax])/d_scale) * 1.6
else:
ymax = np.nanmax(fybb*c/np.square(xbb)/d_scale) * 1.6
else:
ymax = np.nanmax(fybb*c/np.square(xbb)/d_scale) * 1.6
ymax = None

ax1.set_xlabel('Observed wavelength [$\mathrm{\mu m}$]', fontsize=11)
ax1.set_ylabel('$f_\lambda$ [$10^{%d}\mathrm{erg}/\mathrm{s}/\mathrm{cm}^{2}/\mathrm{\AA}$]'%(np.log10(scale)),fontsize=11,labelpad=2)

x1max = 100000
if x1max < np.nanmax(xbb):
x1max = np.nanmax(xbb) * 1.5
if len(fybb[conbb_ymax]):
if x1min > np.nanmin(xbb[conbb_ymax]):
x1min = np.nanmin(xbb[conbb_ymax]) / 1.5
if MB.has_photometry:
if x1max < np.nanmax(xbb):
x1max = np.nanmax(xbb) * 1.5
if len(fybb[conbb_ymax]):
if x1min > np.nanmin(xbb[conbb_ymax]):
x1min = np.nanmin(xbb[conbb_ymax]) / 1.5
else:
x1min = 2000

xticks = [2500, 5000, 10000, 20000, 40000, 80000, x1max]
xlabels= ['0.25', '0.5', '1', '2', '4', '8', '']
Expand All @@ -530,27 +536,13 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=None, f_chind=True, figpdf=Fals
scl_yaxis = 0.2
else:
scl_yaxis = 0.1
ax1.set_ylim(-ymax*scl_yaxis,ymax)

if False:
ax1.text(x1min+100,-ymax*0.08,'SNlimit:%.1f'%(SNlim),fontsize=7)
if not ymax == None:
ax1.set_ylim(-ymax*scl_yaxis,ymax)

ax1.set_xticks(xticks)
ax1.set_xticklabels(xlabels)

if False:
dely1 = 0.5
while (ymax-0)/dely1<1:
dely1 /= 2.
while (ymax-0)/dely1>4:
dely1 *= 2.

y1ticks = np.arange(0, ymax, dely1)
ax1.set_yticks(y1ticks)
ax1.set_yticklabels(np.arange(0, ymax, dely1), minor=False)
ax1.yaxis.set_major_formatter(FormatStrFormatter('%.1f'))
ax1.yaxis.labelpad = 1.5

xx = np.arange(100,400000)
yy = xx * 0
ax1.plot(xx, yy, ls='--', lw=0.5, color='k')
Expand Down Expand Up @@ -612,9 +604,17 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=None, f_chind=True, figpdf=Fals
egrism = np.concatenate([eg0,eg1,eg2])
con4000b = (xgrism/zscl>3400) & (xgrism/zscl<3800) & (fgrism>0) & (egrism>0)
con4000r = (xgrism/zscl>4200) & (xgrism/zscl<5000) & (fgrism>0) & (egrism>0)
print('Median SN at 3400-3800 is;', np.median((fgrism/egrism)[con4000b]))
print('Median SN at 4200-5000 is;', np.median((fgrism/egrism)[con4000r]))

MB.logger.info('Median SN at 3400-3800 is; %.1f'%np.median((fgrism/egrism)[con4000b]))
MB.logger.info('Median SN at 4200-5000 is; %.1f'%np.median((fgrism/egrism)[con4000r]))

if MB.has_spectrum and not MB.has_photometry:
con_spec = (eg2 < 1000)
ax1.errorbar(xg2[con_spec], (fg2 * c/np.square(xg2)/d_scale)[con_spec], yerr=(eg2 * c/np.square(xg2)/d_scale)[con_spec], lw=0.5, color='#DF4E00', zorder=10, alpha=1., label='', capsize=0)
con_spec = (eg1 < 1000)
ax1.errorbar(xg1[con_spec], (fg1 * c/np.square(xg1)/d_scale)[con_spec], yerr=(eg1 * c/np.square(xg1)/d_scale)[con_spec], lw=0.5, color='g', zorder=10, alpha=1., label='', capsize=0)
con_spec = (eg0 < 1000)
ax1.errorbar(xg0[con_spec], (fg0 * c/np.square(xg0)/d_scale)[con_spec], yerr=(eg0 * c/np.square(xg0)/d_scale)[con_spec], lw=0.5, color='royalblue', zorder=10, alpha=1., label='', capsize=0)

#
# From MCMC chain
Expand Down Expand Up @@ -1218,11 +1218,9 @@ def plot_sed(MB, flim=0.01, fil_path='./', scale=None, f_chind=True, figpdf=Fals
if MB.f_dust:
label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log M_\mathrm{dust}/M_\odot:%.2f$\n$T_\mathrm{dust}/K:%.1f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$'\
%(ID, zbes, float(fd['Mstel_50']), MD50, TD50, float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']))#, fin_chi2)
ylabel = ymax*0.45
else:
label = 'ID: %s\n$z:%.2f$\n$\log M_\mathrm{*}/M_\odot:%.2f$\n$\log Z_\mathrm{*}/Z_\odot:%.2f$\n$\log T_\mathrm{*}$/Gyr$:%.2f$\n$A_V$/mag$:%.2f$'\
%(ID, zbes, float(fd['Mstel_50']), float(fd['Z_MW_50']), float(fd['T_MW_50']), float(fd['AV_50']))
ylabel = ymax*0.25

if f_grsm:
ax1.text(0.02, 0.68, label,\
Expand Down Expand Up @@ -3121,9 +3119,13 @@ def density_estimation(m1, m2):
else:
scl_yaxis = 0

ymax_bb = np.max(fybb[conbb] * c / np.square(xbb[conbb]) /d_scale) * 1.10
ymax_temp = np.max(ysum * c/ np.square(x0) /d_scale) * 1.10
ymax = np.max([ymax_bb, ymax_temp])
try:
ymax_bb = np.max(fybb[conbb] * c / np.square(xbb[conbb]) /d_scale) * 1.10
ymax_temp = np.max(ysum * c/ np.square(x0) /d_scale) * 1.10
ymax = np.max([ymax_bb, ymax_temp])
except:
ymax_temp = np.max(ysum * c/ np.square(x0) /d_scale) * 1.10
ymax = np.max(ymax_temp)

# Convert into log
Ztmp[kk] /= ACtmp[kk]
Expand Down
2 changes: 1 addition & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ def read(fname):
license = "IPAC",
url = "https://github.com/mtakahiro",
download_url = "https://github.com/",
packages=['gsf'],
packages=['gsf', 'gsf/Logger/'],
package_dir={'gsf': 'gsf'},
requires=['lmfit', 'fsps', 'emcee', 'corner'],
use_scm_version=True,
Expand Down

0 comments on commit b973ecb

Please sign in to comment.