diff --git a/Tables/Observations/SFR_Mstar/Jiang16.txt b/Tables/Observations/SFR_Mstar/Jiang16.txt new file mode 100644 index 0000000..0380b42 --- /dev/null +++ b/Tables/Observations/SFR_Mstar/Jiang16.txt @@ -0,0 +1,28 @@ +# ID age(M) mass(G) SFR_Lya SFR_UV E(B-V) + 03 4 0.717 9.5 20.9 0.04 + 04 4 0.336 8.0 9.9 0.00 + 15 8 2.113 19.0 33.9 0.08 + 20 492 3.282 15.1 10.5 0.02 + 23 444 5.794 4.6 16.0 0.00 + 24 4 0.377 12.4 15.2 0.00 + 25 956 6.653 13.3 15.9 0.02 + 27 324 1.735 6.3 6.9 0.00 + 28 4 0.268 16.6 10.4 0.00 + 30 912 2.684 4.2 7.8 0.00 + 31 4 0.238 17.1 8.4 0.04 + 34 992 25.000 3.7 35.0 0.12 + 35 4 8.106 16.5 19.5 0.34 + 36 988 4.639 2.5 16.1 0.00 + 43 972 3.076 4.7 8.6 0.02 + 44 872 17.291 18.1 32.6 0.08 + 47 184 39.111 4.9 48.4 0.18 + 49 212 3.556 5.4 9.7 0.08 + 50 4 0.897 8.7 7.7 0.22 + 54 460 5.154 11.8 14.1 0.00 + 58 960 6.162 4.9 17.5 0.04 + 61 32 6.634 11.4 11.0 0.20 + 62 472 11.501 13.7 23.9 0.04 + 63 8 0.945 16.5 22.4 0.02 + 64 200 4.132 11.8 20.2 0.00 + 66 4 0.135 7.6 7.1 0.00 + 67 4 3.201 20.8 43.3 0.14 diff --git a/Tables/observations_z6 b/Tables/observations_z6 new file mode 100644 index 0000000..757a9c9 Binary files /dev/null and b/Tables/observations_z6 differ diff --git a/sigame/analysis.py b/sigame/analysis.py index ca37d82..516d6d7 100644 --- a/sigame/analysis.py +++ b/sigame/analysis.py @@ -27,8 +27,19 @@ """ Global line luminosities """ #------------------------------------------------------------------------------- -def ISM_line_contributions(split=True,line='CII'): - print('Split L_[CII] into contributions from different ISM phases!') +def ISM_line_contributions(line='CII'): + ''' + Purpose + --------- + Splits line luminosity into contributions from different ISM phases + + Arguments + --------- + line: which line to plot - str + default = 'CII' + ''' + + print('') plt.close('all') # close all windows # Load model results @@ -69,32 +80,18 @@ def ISM_line_contributions(split=True,line='CII'): plt.savefig('plots/line_emission/ISM_contributions/['+line+']_lum_fractions_'+z1+ext_DENSE+ext_DIFFUSE+'.pdf', format='pdf') # .eps for paper! - print('Now plotting L_CII fractions vs mass fractions') - - # Starting figure - fontsize = 13 - mpl.rcParams['xtick.labelsize'] = fontsize - mpl.rcParams['ytick.labelsize'] = fontsize - - xr = [0.1,100] - yr = [0.1,100] - # Set up figure - fig = plt.figure(1,figsize = (8,7)) - ax1 = fig.add_subplot(1,1,1) - # xlog='y',ylog='y',\ - # colored by ISM phase version - plot.simple_plot(add='y',xr=xr, yr=yr,\ - xlab='Mass fraction [%]',ylab='Fraction of total L$_{[\mathrm{CII}]}$ [%]',\ - x1=f_M[0],y1=f_L[0],ma1='x',fill1='y',ms1=6,mew1=2,col1='b',\ - x2=f_M[1],y2=f_L[1],ma2='x',fill2='y',ms2=6,mew2=2,col2='orange',\ - x3=f_M[2],y3=f_L[2],ma3='x',fill3='y',ms3=6,mew3=2,col3='r') - +def ISM_line_efficiency(line='CII'): + ''' + Purpose + --------- + Plots the line efficiency (line luminosity / gas mass) in different ISM phases - plt.show(block=False) - plt.savefig('plots/line_emission/ISM_contributions/['+line+']_lum_mass_fractions_'+z1+ext_DENSE+ext_DIFFUSE+'.pdf', format='pdf') # .eps for paper! + Arguments + --------- + line: which line to plot - str + default = 'CII' + ''' -def ISM_line_efficiency(line='CII'): - print('Split L_[CII] into contributions from different ISM phases!') plt.close('all') # close all windows # Load model results @@ -164,7 +161,13 @@ def ISM_line_efficiency(line='CII'): #------------------------------------------------------------------------------- def ISM_mass_contributions(): - print('Investigate mass from different ISM phases!') + ''' + Purpose + --------- + Plots the mass contribution from different ISM phases + + ''' + plt.close('all') # close all windows # Load model results diff --git a/sigame/plot.py b/sigame/plot.py index 667c67a..f5ce7d6 100644 --- a/sigame/plot.py +++ b/sigame/plot.py @@ -859,7 +859,6 @@ def SFR_Mstar(): Plot of SFR vs stellar mass ''' - print('** Plot SFR vs stellar mass for sample **') plt.close('all') models = aux.find_model_results() @@ -883,25 +882,18 @@ def SFR_Mstar(): age = 0.984 # age of universe, Omega_m = 0.3, Omega_lambda = 0.7, h = 0.65 MS = 10.**((0.84-0.026*age)*np.log10(xr)-(6.51-0.11*age)) # z ~ 6 LBGs and LAEs - # columns = ['No','M_1500','UV_slope','L_Lyalpha','ew_Lyalpha','SFR_UV','SFR_Lyalpha'] - # J16 = pd.read_table('/Users/Karen/Desktop/temp/Observations/Jiang16_paper',names=columns,skiprows=5,sep=r'\t*',engine='python') - # columns = ['ID','age(M)','mass(G)','SFR_Lya','SFR_UV','E(B-V)'] - # J16 = pd.read_table('sigame/temp/Observations/Jiang16.txt',names=columns,skiprows=1,sep=r'\s*',engine='python') - # SFR_J16 = J16['SFR_UV'].values - # No_J16 = J16['ID'].values - # ebv = J16['E(B-V)'].values - # age_J16 = J16['age(M)'].values - # M_star_J13 = np.array([7.2,3.4,21.1,32.8,57.9,3.8,66.5,17.4,2.7,26.8,2.4,250,81.1,46.4,30.8,172.9,391.1,35.6,9.0,51.5,61.6,66.3,115,9.4,41.3,1.4,32])*1e8 - # No_J13 = np.array([3,4,15,20,23,24,25,27,28,30,31,34,35,36,43,44,47,49,50,54,58,61,62,63,64,66,67]) - # SFR_J16 = [float(SFR_J16[No_J16 == No][0]) for No in No_J13] - # # From Linhua Jiang - # klambda = 8.5 #at 2200A - # SFR_J16 = SFR_J16 * 10.**(0.4*klambda*ebv) - - if z1 == 'z2': - # eq. 28 - age = 3.223 # age of universe [Gyr], Omega_m = 0.3, Omega_lambda = 0.7, h = 0.70 - SFR_MS = 10.**((0.84-0.026*age)*np.log10(xr)-(6.51-0.11*age)) + columns = ['ID','age(M)','mass(G)','SFR_Lya','SFR_UV','E(B-V)'] + J16 = pd.read_table(d_t+'Observations/SFR_Mstar/Jiang16.txt',names=columns,skiprows=1,sep=r'\s*',engine='python') + SFR_J16 = J16['SFR_UV'].values + No_J16 = J16['ID'].values + ebv = J16['E(B-V)'].values + age_J16 = J16['age(M)'].values + M_star_J13 = np.array([7.2,3.4,21.1,32.8,57.9,3.8,66.5,17.4,2.7,26.8,2.4,250,81.1,46.4,30.8,172.9,391.1,35.6,9.0,51.5,61.6,66.3,115,9.4,41.3,1.4,32])*1e8 + No_J13 = np.array([3,4,15,20,23,24,25,27,28,30,31,34,35,36,43,44,47,49,50,54,58,61,62,63,64,66,67]) + SFR_J16 = [float(SFR_J16[No_J16 == No][0]) for No in No_J13] + # From Linhua Jiang + klambda = 8.5 #at 2200A + SFR_J16 = SFR_J16 * 10.**(0.4*klambda*ebv) simple_plot(fignum=0,xlog='y',ylog='y',xr=xr,yr=yr,fontsize=16,\ xlab='M$_{\mathrm{*}}$ [M$_{\odot}$]',ylab='SFR [M$_{\odot}$ yr$^{-1}$]',legloc=[0.01,0.75],frameon=False,\ @@ -913,14 +905,14 @@ def SFR_Mstar(): x6=xr,y6=10.**(np.log10(MS)+3.*0.2),col6='k',lw6=1,ls6=':',legend=False) ax1 = plt.gca() - # if z1 == 'z6': - # simple_plot(add='y',x1=M_star_J13[age_J16 > 30],y1=SFR_J16[age_J16 > 30],fill1='y',ma1='x',col1='r',ms1=10,mew1=2,lab1='Old $z\sim6$ LBGs/LAEs [Jiang+16]',\ - # x2=M_star_J13[age_J16 < 30],y2=SFR_J16[age_J16 < 30],fill2='y',ma2='+',col2='b',ms2=11,mew2=2,lab2='Young $z\sim6$ LBGs/LAEs [Jiang+16]') + if z1 == 'z6': + simple_plot(add='y',x1=M_star_J13[age_J16 > 30],y1=SFR_J16[age_J16 > 30],fill1='y',ma1='x',col1='r',ms1=10,mew1=2,lab1='Old $z\sim6$ LBGs/LAEs [Jiang+16]',\ + x2=M_star_J13[age_J16 < 30],y2=SFR_J16[age_J16 < 30],fill2='y',ma2='+',col2='b',ms2=11,mew2=2,lab2='Young $z\sim6$ LBGs/LAEs [Jiang+16]',legloc='upper left') # Plot models SC = ax1.scatter(M_star,SFR,marker='o',lw=2,s=64,c=ages,edgecolor='black',cmap='viridis',alpha=1,label='',zorder=10) cbar = plt.colorbar(SC,pad=0) - cbar.set_label(label='Mass-weighted stellar age [Myr]')#,size=12) # colorbar in it's own axis + cbar.set_label(label='Mass-weighted stellar age [Myr]') # colorbar in it's own axis plt.tight_layout() plt.show(block=False) @@ -942,35 +934,31 @@ def CII_SFR_z6(plot_models=True,plot_fits=True,plot_obs=True,twopanel=True,obs_g models = aux.find_model_results() for key in models.keys(): exec(key + '=models[key].values') + L_CII_tot = L_CII_GMC+L_CII_DNG+L_CII_DIG z1 = 'z6' + # Plotting parameters mpl.rcParams['xtick.minor.size'] = 4 mpl.rcParams['xtick.major.size'] = 6 mpl.rcParams['xtick.minor.width'] = 1.5 mpl.rcParams['xtick.major.width'] = 1.5 fs_labels = 20 if twopanel: fs_labels = 17 - - # Plot parameters, quickly change y or x range here + # Change y or x range here xr = 10.**np.array([0,3.05]) # SFR range yr = 10.**np.array([6.5,9.8]) # [CII] range - if not twopanel: fig = plt.figure(figsize = (12,9)) - if twopanel: fig = plt.figure(figsize = (18,8)) - - # Load [CII] luminosities - sims = cl.sim(sim_paths) - L_CII_tot = L_CII_GMC+L_CII_DNG+L_CII_DIG - - if not twopanel: - ax1 = fig.add_subplot(1,1,1) if twopanel: + fig = plt.figure(figsize = (18,8)) gs = gridspec.GridSpec(1, 2, width_ratios=[6, 5]) ax1 = plt.subplot(gs[1]) + if not twopanel: + fig = plt.figure(figsize = (12,9)) + ax1 = fig.add_subplot(1,1,1) + ax1.set_xlabel('SFR [M$_{\odot}$ yr$^{-1}$]',fontsize=fs_labels) ax1.set_ylabel('L$_{[\mathrm{CII}]}$ [L$_{\odot}$]',fontsize=fs_labels) ax1.set_xlim(xr) ax1.set_ylim(yr) - ax1.set_xlabel('SFR [M$_{\odot}$ yr$^{-1}$]',fontsize=fs_labels) xlab = ax1.get_xticks() ylab = ax1.get_yticks() ax1.set_xticklabels(xlab,fontsize=fs_labels) @@ -984,46 +972,37 @@ def CII_SFR_z6(plot_models=True,plot_fits=True,plot_obs=True,twopanel=True,obs_g # Models! # ====================================================================== - # Powerlaw fit to high-z observations [de Looze et al. 2014] - # logL_CII = np.array([4,11]) - # logSFR = -8.52+1.18*logL_CII - # fit0 = ax1.plot(10.**logSFR,10.**logL_CII,c='purple',ls='--',lw=lw) - # Powerlaw fit to z=0 metal-poor dwarf galaxies + # Powerlaw fit to z=0 metal-poor dwarf galaxies [de Looze et al. 2014] logL_CII = np.array([4,11]) logSFR = -5.73+0.80*logL_CII ax1.fill_between(10.**logSFR, 10.**(logL_CII-0.37), 10.**(logL_CII+0.37), facecolor='lightgrey', linewidth=0, alpha=0.5) fit0 = ax1.plot(10.**logSFR,10.**logL_CII,c='grey',ls='--',lw=lw) - # Powerlaw fit to z=0 starburst galaxies + # Powerlaw fit to z=0 starburst galaxies [de Looze et al. 2014] logL_CII = np.array([4,11]) logSFR = -7.06+1.00*logL_CII ax1.fill_between(10.**logSFR, 10.**(logL_CII-0.27), 10.**(logL_CII+0.27), facecolor='lightgrey', linewidth=0, alpha=0.5) fit1 = ax1.plot(10.**logSFR,10.**logL_CII,c='grey',ls='--',dashes=(1,4),lw=lw) - print('\nAverage factor below de Looze+14 metal-poor dwarfs') - - # Powerlaw fit to z=0 observations [Malhotra et al. 2001] - # L_FIR = xr*3.4e9 - # L_CII = 10.**(-2.51+np.log10(L_FIR)) - # fit1 = ax1.plot(xr,L_CII,'--',dashes=(2,5),c='teal',lw=lw) # Powerlaw fit to model galaxies at z=6 slope_models,intercept_models,slope_dev,inter_dev = aux.lin_reg_boot(np.log10(SFR),np.log10(L_CII_tot)) + # Making different line styles for the SFR-range covered by the simulations and + # the SFR-range where an extrapolation is necessary L_CII_model = 10.**(slope_models*np.log10(xr)+intercept_models) xr2 = [min(SFR),max(SFR)] L_CII_model2 = 10.**(slope_models*np.log10(xr2)+intercept_models) # dex spread around relation: dex = np.std(np.abs(np.log10(L_CII_tot)-(slope_models*np.log10(SFR)+intercept_models))) print('Dex: %s' % dex) - if not twopanel: if plot_models: fit = ax1.plot(xr,L_CII_model,ls='--',color='purple',lw=1.8) fit = ax1.plot(xr2,L_CII_model2,ls='-',color='purple',lw=1.8) fit = ax1.plot(xr,L_CII_model,ls='--',color='purple',lw=1.8) fit = ax1.plot(xr2,L_CII_model2,ls='-',color='purple',lw=1.8) - if plot_fits: # Powerlaw fit to Vallini+15, z~7 eq. 8 + sims = cl.sim(sim_paths) Zfit = np.mean(sims.Zmw()) L_CII_V15 = 10.**(7.0 + 1.2*np.log10(xr) + 0.021*np.log10(Zfit) + \ 0.012*np.log10(xr)*np.log10(Zfit) - 0.74*np.log10(Zfit)**2.) @@ -1033,7 +1012,7 @@ def CII_SFR_z6(plot_models=True,plot_fits=True,plot_obs=True,twopanel=True,obs_g # PCA fit to L_[CII](SFR,Z) Zfit2 = np.mean(Zsfr) - L_CII_PCA = 10.**(7.26+0.35*np.log10(xr)+0.38*np.log10(Zfit2)) + L_CII_PCA = 10.**(7.17+0.55*np.log10(xr)+0.23*np.log10(Zfit2)) ax1.plot(xr,L_CII_PCA,c='purple',ls='--',dashes=(6,4,1,4,1,4),lw=lw) if len(galnames) > 1: @@ -1050,19 +1029,20 @@ def CII_SFR_z6(plot_models=True,plot_fits=True,plot_obs=True,twopanel=True,obs_g rms_dex = np.sqrt(np.sum(dev**2./len(dev))) print('rms error of: '+str(rms_dex)+' dex') - print('On average this far below observed relation:') - L_CII_obs = 10.**((np.log10(SFR)+7.06)/1.00) - print(str(np.mean((L_CII_obs-L_CII_tot)/L_CII_obs)*100.)+' %') - print(str(np.mean(np.log10(L_CII_obs)-np.log10(L_CII_tot)))+' dex') - print('Power law fit this far above Vallini+15:') - L_CII_V15 = 10.**(7.0 + 1.2*np.log10(min(SFR)) + 0.021*np.log10(Zfit) + \ - 0.012*np.log10(min(SFR))*np.log10(Zfit) - 0.74*np.log10(Zfit)**2.) - L_CII = 10.**(slope_models*np.log10(min(SFR))+intercept_models) - print('at min SFR : '+str(np.log10(L_CII)-np.log10(L_CII_V15))+' dex') - L_CII_V15 = 10.**(7.0 + 1.2*np.log10(max(SFR)) + 0.021*np.log10(Zfit) + \ - 0.012*np.log10(max(SFR))*np.log10(Zfit) - 0.74*np.log10(Zfit)**2.) - L_CII = 10.**(slope_models*np.log10(max(SFR))+intercept_models) - print('at max SFR : '+str(np.log10(L_CII)-np.log10(L_CII_V15))+' dex') + if len(galnames) > 1: + print('On average this far below observed relation:') + L_CII_obs = 10.**((np.log10(SFR)+7.06)/1.00) + print(str(np.mean((L_CII_obs-L_CII_tot)/L_CII_obs)*100.)+' %') + print(str(np.mean(np.log10(L_CII_obs)-np.log10(L_CII_tot)))+' dex') + print('Power law fit this far above Vallini+15:') + L_CII_V15 = 10.**(7.0 + 1.2*np.log10(min(SFR)) + 0.021*np.log10(Zfit) + \ + 0.012*np.log10(min(SFR))*np.log10(Zfit) - 0.74*np.log10(Zfit)**2.) + L_CII = 10.**(slope_models*np.log10(min(SFR))+intercept_models) + print('at min SFR : '+str(np.log10(L_CII)-np.log10(L_CII_V15))+' dex') + L_CII_V15 = 10.**(7.0 + 1.2*np.log10(max(SFR)) + 0.021*np.log10(Zfit) + \ + 0.012*np.log10(max(SFR))*np.log10(Zfit) - 0.74*np.log10(Zfit)**2.) + L_CII = 10.**(slope_models*np.log10(max(SFR))+intercept_models) + print('at max SFR : '+str(np.log10(L_CII)-np.log10(L_CII_V15))+' dex') # Powerlaw fit to model galaxies at z=2 [Olsen+15] L_CII_model_z2 = 0.78e7*xr**1.27 @@ -1113,11 +1093,9 @@ def CII_SFR_z6(plot_models=True,plot_fits=True,plot_obs=True,twopanel=True,obs_g ax1.fill_between(10.**logSFR, 10.**(logL_CII-0.27), 10.**(logL_CII+0.27), facecolor='lightgrey', linewidth=0, alpha=0.5,zorder=0) fit1 = ax1.plot(10.**logSFR,10.**logL_CII,c='grey',ls='--',dashes=(1,4),lw=lw,zorder=0) - # pdb.set_trace() if plot_obs: add_observations_to_plot(slope_models,intercept_models,mark_reasons,mark_uplim,mark_det,z1=z1) # Legends for observations! - # get handles if plot_models: ax1.scatter([1e20,1e30],[1e20,1e30],marker='o',s=64,color='lightseagreen',label=r'S$\mathrm{\'I}$GAME (this work)$_{}$',lw=2,edgecolor='black') handles2, labels2 = ax1.get_legend_handles_labels() # remove the errorbars @@ -1130,14 +1108,12 @@ def CII_SFR_z6(plot_models=True,plot_fits=True,plot_obs=True,twopanel=True,obs_g ax1.legend(handles2, labels2, loc='upper left',fontsize=11,frameon=False,numpoints=1,borderpad=1) # Model galaxies - # L_CII_tot = L_CII_tot*0.+10.**6.6 ax1.xaxis.set_major_formatter(ScalarFormatter()) if plot_models: plot = ax1.scatter(SFR,L_CII_tot,marker='o',edgecolor='black',c=np.log10(SFRsd),s=64,lw=2,cmap='viridis',zorder=200) cbar = fig.colorbar(plot,ax=ax1,pad=0)#,cax = cbar_ax) cbar.set_label(label=getlabel('lSFRsd'),size=13)#,weight='bold') - for axis in ['top','bottom','left','right']: ax1.spines[axis].set_linewidth(1.) @@ -1599,10 +1575,8 @@ def grid_parameters(histo_color='teal',FUV=5,ISM_phase='GMC',figsize=(10,7)): if histo_color == 'colsel': ax1.legend(loc='upper left') panel += 1 - # plt.suptitle('Grid points (dashed lines) on histograms (solid lines) in diffuse gas') plt.show(block=False) plt.savefig('plots/diffuse_gas/grid_parameters/difgrid_points_on_histos'+ext_DIFFUSE+'_'+z1+'.pdf',format='pdf') # .eps for paper! - # fig.savefig('plots/dif/grids/grid_points_on_histos_paper.png', dpi=fig.dpi) #=============================================================================== """ Gas distribution plots """