diff --git a/gsf/fitting.py b/gsf/fitting.py index 88773e7..eeee666 100644 --- a/gsf/fitting.py +++ b/gsf/fitting.py @@ -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)): @@ -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']=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) @@ -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: @@ -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' diff --git a/gsf/maketmp_filt.py b/gsf/maketmp_filt.py index 317440b..7db2e8e 100644 --- a/gsf/maketmp_filt.py +++ b/gsf/maketmp_filt.py @@ -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])) @@ -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') diff --git a/gsf/plot_sed.py b/gsf/plot_sed.py index cdb55c5..cde4331 100644 --- a/gsf/plot_sed.py +++ b/gsf/plot_sed.py @@ -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', ''] @@ -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') @@ -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 @@ -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,\ @@ -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] diff --git a/setup.py b/setup.py index 4620a82..4a403b1 100644 --- a/setup.py +++ b/setup.py @@ -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,