Skip to content

Commit

Permalink
Merge pull request #49 from mtakahiro/nirspec_res
Browse files Browse the repository at this point in the history 
NIRSpec Prism and some bug fixed
  • Loading branch information
@mtakahiro
mtakahiro authored May 11, 2023
2 parents 86b7012 + b99ffee commit 9af5c35
Show file tree
Hide file tree
Showing 6 changed files with 282 additions and 280 deletions.
40 changes: 14 additions & 26 deletions gsf/fitting.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -992,15 +992,17 @@ def fit_redshift(self, xm_tmp, fm_tmp, delzz=0.01, ezmin=0.01, #zliml=0.01, zlim
# Observed data;
# spec;
if self.has_spectrum:
con = (self.dict['ey']<1000) & (self.dict['NR']<self.NRbb_lim)
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;
con = (self.dict['NR']>=self.NRbb_lim)
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)
ax1.scatter(self.dict['x'][con], self.dict['fy'][con], s=100, marker='o',
color='orange', edgecolor='k', label='Observed photometry', zorder=4)
if include_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)
ax1.scatter(self.dict['x'][con], self.dict['fy'][con], s=100, marker='o',
color='orange', edgecolor='k', label='Observed photometry', zorder=4)

# Write prob distribution;
get_chi2(zz_prob, fy_cz, ey_cz, x_cz, fm_tmp, xm_tmp/(1+self.zgal), self.file_zprob)
Expand Down Expand Up @@ -1322,9 +1324,7 @@ def set_param(self):
'''
Set parameters.
'''
print('##################')
print('Setting parameters')
print('##################\n')
self.logger.info('Setting parameters')
agemax = self.cosmo.age(self.zgal).value
fit_params = Parameters()
f_Alog = True
Expand Down Expand Up @@ -1482,9 +1482,7 @@ def check_mainbody(self):
def prepare_class(self, add_fir=None):
'''
'''
print('#################')
print('Preparing library')
print('#################\n')
self.logger.info('Preparing library')
# Load Spectral library;
self.lib = self.fnc.open_spec_fits(fall=0)
self.lib_all = self.fnc.open_spec_fits(fall=1, orig=True)
Expand Down Expand Up @@ -1702,21 +1700,11 @@ def main(self, cornerplot:bool=True, specplot=1, sigz=1.0, ezmin=0.01, ferr=0,
except:
out.params['Z0'].value = out_keep.params['Z0'].value


##############################
print('\n\n')
print('###############################')
print('Input redshift is adopted.')
print('Starting long journey in MCMC.')
print('###############################')
print('\n\n')

################################
self.logger.info('Minimizer Defined')

print('########################')
print('### Starting sampling ##')
print('########################\n')
self.logger.info('Input redshift is adopted')
self.logger.info('Starting MCMC')
# self.logger.info('Minimizer Defined')
# self.logger.info('Starting sampling')
start_mc = timeit.default_timer()

# MCMC;
Expand Down
2 changes: 1 addition & 1 deletion gsf/gsf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ def run_gsf_all(parfile, fplt, cornerplot=True, f_plot_chain=True, f_Alog=True,
f_prior_sfh=False, norder_sfh_prior=3,
f_shuffle=False, amp_shuffle=1e-2, Zini=None, tau_lim=0.001,
skip_sfh=False, f_fancyplot=False, skip_zhist=False, f_sfh_yaxis_force=True, tset_SFR_SED=0.1,
nthin=1, delwave=1, f_plot_resid=False, scale=1e-19, f_plot_filter=True,
nthin=1, delwave=1, f_plot_resid=False, scale=None, f_plot_filter=True,
mmax_param:int=1000
):
'''
Expand Down
153 changes: 91 additions & 62 deletions gsf/maketmp_filt.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,49 +141,51 @@ def check_library(MB, af, nround=3):

flag = True
MB.logger.info('Checking the template library...')
print('Speficied - Template')

# No. of age;
if MB.SFH_FORM==-99:
if len(af['ML']['ms_0']) != len(MB.age):
print('No of age pixels:', len(MB.age), len(af['ML']['ms_0']))
MB.logger.error('No of age pixels:', len(MB.age), len(af['ML']['ms_0']))
flag = False
else:
flag = True

# Matallicity:
for aa in range(len(Zall)):
if Zall[aa] != af['Z%d'%(aa)]:
print('Z:', Zall[aa], af['Z%d'%(aa)])
MB.logger.error('Z:', Zall[aa], af['Z%d'%(aa)])
flag = False

if MB.SFH_FORM==-99:
# Age:
for aa in range(len(MB.age)):
if round(MB.age[aa],nround) != round(af['age%d'%(aa)],nround):
print('age:', MB.age[aa], af['age%d'%(aa)])
MB.logger.error('age:', MB.age[aa], af['age%d'%(aa)])
flag = False
# Tau (e.g. ssp/csp):
for aa in range(len(MB.tau0)):
if round(MB.tau0[aa]) != round(af['tau0%d'%(aa)]):
print('tau0:', MB.tau0[aa], af['tau0%d'%(aa)])
MB.logger.error('tau0:', MB.tau0[aa], af['tau0%d'%(aa)])
flag = False
else:
# Age:
for aa in range(len(MB.ageparam)):
if round(MB.ageparam[aa]) != round(af['age%d'%(aa)]):
print('age:', MB.ageparam[aa], af['age%d'%(aa)])
MB.logger.error('age:', MB.ageparam[aa], af['age%d'%(aa)])
flag = False
for aa in range(len(MB.tau)):
if round(MB.tau[aa]) != round(af['tau%d'%(aa)]):
print('tau:', MB.tau[aa], af['tau%d'%(aa)])
MB.logger.error('tau:', MB.tau[aa], af['tau%d'%(aa)])
flag = False

# IMF:
if MB.nimf != af['nimf']:
print('nimf:', MB.nimf, af['nimf'])
MB.logger.error('nimf:', MB.nimf, af['nimf'])
flag = False

if not flag:
MB.logger.error('# Specified - Template')

return flag


Expand Down Expand Up @@ -512,64 +514,91 @@ def maketemp(MB, ebblim=1e10, lamliml=0., lamlimu=50000., ncolbb=10000,

try:
spec_files = [x.strip() for x in inputs['SPEC_FILE'].split(',')]
ninp0 = np.zeros(len(spec_files), dtype='int')
except:
spec_files = []
MB.logger.info('No spec file is provided.')
pass
ninp0 = np.zeros(len(spec_files), dtype='int')

# THIS PART IS JUST TO GET THE TOTAL ARRAY NUMBER;
for ff, spec_file in enumerate(spec_files):
try:
if spec_file.split('.')[-1] == 'asdf':
id_asdf = int(spec_file.split('_')[2])
fd0 = asdf.open(os.path.join(DIR_EXTR, spec_file))
lm0tmp = fd0[id_asdf]['wavelength'].to(u.angstrom)
ninp0[ff] = len(lm0tmp)
else:
fd0 = np.loadtxt(os.path.join(DIR_EXTR, spec_file), comments='#')
lm0tmp = fd0[:,0]
ninp0[ff] = len(lm0tmp)
except Exception:
MB.logger.warning('File, %s, cannot be open.'%(os.path.join(DIR_EXTR, spec_file)))
pass
# THIS PART IS JUST TO GET THE TOTAL ARRAY NUMBER;
for ff, spec_file in enumerate(spec_files):
try:
if spec_file.split('.')[-1] == 'asdf':
id_asdf = int(spec_file.split('_')[2])
fd0 = asdf.open(os.path.join(DIR_EXTR, spec_file))
lm0tmp = fd0[id_asdf]['wavelength'].to(u.angstrom)
ninp0[ff] = len(lm0tmp)
elif spec_file.split('.')[-1] == 'fits':
fd0 = fits.open(os.path.join(DIR_EXTR, spec_file))[1].data
eobs0 = fd0['full_err']
spec_mask = (eobs0>0)
lm0tmp = fd0['wave'][spec_mask]
ninp0[ff] = len(lm0tmp)
else:
fd0 = np.loadtxt(os.path.join(DIR_EXTR, spec_file), comments='#')
lm0tmp = fd0[:,0]
ninp0[ff] = len(lm0tmp)
except Exception:
MB.logger.error('File, %s, cannot be open.'%(os.path.join(DIR_EXTR, spec_file)))
pass

# Then, Constructing arrays.
lm = np.zeros(np.sum(ninp0[:]),dtype=float)
fobs = np.zeros(np.sum(ninp0[:]),dtype=float)
eobs = np.zeros(np.sum(ninp0[:]),dtype=float)
fgrs = np.zeros(np.sum(ninp0[:]),dtype=int) # FLAG for each grism.
for ff, spec_file in enumerate(spec_files):
if True:#try:
if spec_file.split('.')[-1] == 'asdf':
id_asdf = int(spec_file.split('_')[2])
fd0 = asdf.open(os.path.join(DIR_EXTR, spec_file))
lm0tmp = fd0[id_asdf]['wavelength'].to(u.angstrom).value
fobs0 = fd0[id_asdf]['flux'].value
eobs0 = np.sqrt(fd0[id_asdf]['fluxvar']).value
# Then, Constructing arrays.
lm = np.zeros(np.sum(ninp0[:]),dtype=float)
fobs = np.zeros(np.sum(ninp0[:]),dtype=float)
eobs = np.zeros(np.sum(ninp0[:]),dtype=float)
fgrs = np.zeros(np.sum(ninp0[:]),dtype=int) # FLAG for each grism.
for ff, spec_file in enumerate(spec_files):
try:
if spec_file.split('.')[-1] == 'asdf':
id_asdf = int(spec_file.split('_')[2])
fd0 = asdf.open(os.path.join(DIR_EXTR, spec_file))
lm0tmp = fd0[id_asdf]['wavelength'].to(u.angstrom).value
fobs0 = fd0[id_asdf]['flux'].value
eobs0 = np.sqrt(fd0[id_asdf]['fluxvar']).value
elif spec_file.split('.')[-1] == 'fits':
fd0 = fits.open(os.path.join(DIR_EXTR, spec_file))[1].data
eobs0 = fd0['full_err']
if True:
spec_mask = (eobs0>0)
else:
fd0 = np.loadtxt(os.path.join(DIR_EXTR, spec_file), comments='#')
lm0tmp = fd0[:,0]
fobs0 = fd0[:,1]
eobs0 = fd0[:,2]

for ii1 in range(ninp0[ff]):
if ff==0:
ii = ii1
else:
ii = ii1 + np.sum(ninp0[:ff])
fgrs[ii] = ff
lm[ii] = lm0tmp[ii1]
fobs[ii] = fobs0[ii1]
eobs[ii] = eobs0[ii1]

MB.f_spec = True
data_meta['data_len'][ff] = len(lm0tmp)
data_meta['data_origin'] = np.append(data_meta['data_origin'], '%s'%spec_file)
data_meta['data_index'] = np.append(data_meta['data_index'], '%d'%ff)
spec_mask = ()
lm0tmp = fd0['wave'][spec_mask]
if lm0tmp.max() < 10:
MB.logger.warning('Wave column in the input spec file seems to be um. Scaling to AA.')
lm0tmp *= 1e4
try:
magzp_spec = float(inputs['MAGZP_SPEC'])
magzp = float(inputs['MAGZP'])
C_spec = 10**((magzp-magzp_spec)/(2.5))
except:
MB.logger.info('`MAGZP_SPEC` not found. Assuming same as `MAGZP`')
C_spec = 1
fobs0 = fd0['flux'][spec_mask] * C_spec
eobs0 = fd0['full_err'][spec_mask] * C_spec
else:
fd0 = np.loadtxt(os.path.join(DIR_EXTR, spec_file), comments='#')
lm0tmp = fd0[:,0]
fobs0 = fd0[:,1]
eobs0 = fd0[:,2]

else:#except Exception:
print('No spec data is registered.')
pass
except:
MB.logger.info('No spec file is provided.')
pass
for ii1 in range(ninp0[ff]):
if ff==0:
ii = ii1
else:
ii = ii1 + np.sum(ninp0[:ff])
fgrs[ii] = ff
lm[ii] = lm0tmp[ii1]
fobs[ii] = fobs0[ii1]
eobs[ii] = eobs0[ii1]

data_meta['data_len'][ff] = len(lm0tmp)
data_meta['data_origin'] = np.append(data_meta['data_origin'], '%s'%spec_file)
data_meta['data_index'] = np.append(data_meta['data_index'], '%d'%ff)
MB.f_spec = True

except Exception:
print('No spec data is registered.')
pass

if ncolbb < np.sum(data_meta['data_len']):
MB.logger.info('ncolbb is updated')
Expand Down
Loading

0 comments on commit 9af5c35

Please sign in to comment.