In [1]:
%pylab inline
import os
from matplotlib import ticker
from astropy.table import Table,join
from simqso.sqgrids import render_gaussians,generateVdBCompositeEmLines


Populating the interactive namespace from numpy and matplotlib

In [2]:
trendsv5 = Table.read('../simqso/data/emlinetrends_v5.fits')
trendsv5['wavelength'][-6,:,1] = 18735.5 # wrong wavelength
trendsv6 = Table.read('../simqso/data/emlinetrends_v6.fits')

In [3]:
def condense_table(trends,M_i=-25):
    tab = Table()
    tab['name'] = trends['name']
    tab['wave'] = polyval(trends['wavelength'][:,0].T,M_i)
    tab['logew'] = polyval(trends['logEW'][:,0].T,M_i)
    tab['ew'] = 10**tab['logew']
    tab['slope'] = trends['logEW'][:,0,0]
    tab['slope'][tab['slope']<0.01] = 0
    tab['sigma'] = polyval(trends['logWidth'][:,0].T,M_i)
    tab['scatter'] = abs(polyval(trends['logEW'][:,1:].T,M_i) - tab['logew']).mean(axis=0)
    return tab

In [4]:
emlinesv5 = condense_table(trendsv5)
emlinesv6 = condense_table(trendsv6)

In [5]:
for i in range(len(emlinesv5)):
    if emlinesv5['name'][i] not in emlinesv6['name']:
        dwave = np.abs(emlinesv5['wave'][i]-emlinesv6['wave'])
        j = dwave.argmin()
        if abs(dwave[j]) < 2:
            print 'replacing ',emlinesv5['name'][i],' with ',emlinesv6['name'][j],dwave[j]
            emlinesv5['name'][i] = emlinesv6['name'][j]
        else:
            print emlinesv5['name'][i],' not in v6',dwave[j]


replacing  Hdelta  with  Hd 0.27002
replacing  Hgamma  with  Hg 0.419922
FeII_Opt37  not in v6 201.15
replacing  [O III]5008  with  [OIII]5008 0.220215
FeII_Opt48  not in v6 297.78
replacing  Halpha  with  HAn 0.0600586
replacing  HeI_7065  with  HeI7067 0.330078
replacing  OI_8446  with  OI8446 0.5
replacing  FeII_9202  with  FeII9202 0.0
replacing  Pa_e  with  Pae 0.400391
replacing  Pa_a  with  Paalpha 0.0
replacing  Pa_d  with  Pad 0.299805
replacing  HeI_10830  with  HeI10830 0.900391
Pa_g  not in v6 79.2002
replacing  Pa_b  with  Pabeta 0.299805
replacing  HeI_20580  with  HeI20580 0.199219

In [6]:
v6tov5 = join(emlinesv6,emlinesv5,'name',join_type='outer',
                 table_names=['v6','v5'])
v6tov5['wave_v6'][v6tov5['wave_v6'].mask] = v6tov5['wave_v5'][v6tov5['wave_v6'].mask]
v6tov5.rename_column('wave_v6','wave')
del v6tov5['wave_v5']
v6tov5.sort('wave')

In [7]:
v6tov5.show_in_notebook()


Out[7]:
<Table masked=True length=65>
idxnamewavelogew_v6ew_v6slope_v6sigma_v6scatter_v6logew_v5ew_v5slope_v5sigma_v5scatter_v5
0LyB1033.830.9966139.92230.04881390.8108570.06716920.9966139.92230.04881390.8108570.0671692
1ArI1068.970.9269248.45130.02931741.172680.04906640.9269248.45130.02931741.172680.0490664
2FeIII:UV11122.010.6553324.522010.03675950.9682470.1497570.6553324.522010.03675950.9682470.149757
3CIII*1174.760.6055634.032390.01.034890.2404380.6055634.032390.01.034890.240438
4LyAn1216.311.2116516.27970.07561950.6249280.09605471.2116516.27970.07561950.6249280.0960547
5LyAb1216.771.848370.51810.09844441.090650.1423721.848370.51810.09844441.090650.142372
6NV1239.281.1265713.38340.00.9419680.05804011.1265713.38340.00.9419680.0580401
7SiII1261.670.8703267.418670.00.9484250.04452750.8703267.418670.00.9484250.0445275
8OI1304.120.5812913.813210.00.9617980.08982510.5812913.813210.00.9617980.0898251
9CII1337.910.5163373.28350.02522751.04430.1119410.5163373.28350.02522751.04430.111941
10SiIV+OIV]1398.691.1279513.42620.01751861.132860.0331691.1279513.42620.01751861.132860.033169
11L14801487.870.1758081.499020.01.021520.05822950.1758081.499020.01.021520.0582295
12CIVn1546.551.0804912.03620.1140110.7479550.163841.0804912.03620.1140110.7479550.16384
13CIVb1548.251.7565757.09150.09233571.311520.05810261.7565757.09150.09233571.311520.0581026
14HeII1636.230.8362276.858460.07143720.9153790.03315120.8362276.858460.07143720.9153790.0331512
15OIII]1666.390.8630687.295720.04039011.109650.02846280.8630687.295720.04039011.109650.0284628
16L16901691.30.5655093.677130.2006240.9534210.01918260.5655093.677130.2006240.9534210.0191826
17NIII]1746.340.6157494.128080.01584370.9705290.0291740.6157494.128080.01584370.9705290.029174
18SiII_18181813.20.5578463.612820.01.143080.06611610.5578463.612820.01.143080.0661161
19AlIII1861.660.7779815.997650.01.15880.02808310.7779815.997650.01.15880.0280831
20SiIII]1892.7-0.03357930.9255940.1030080.6522250.0535518-0.03357930.9255940.1030080.6522250.0535518
21CIII]b1904.121.3497722.37540.020641.259960.01482311.3497722.37540.020641.259960.0148231
22CIII]n1906.870.2701471.862720.00.7226770.06400380.2701471.862720.00.7226770.0640038
23fe21202120.00.2304491.70.01.431360.4604090.2304491.70.01.431360.460409
24fe22202220.00.4771213.00.01.778150.1839880.4771213.00.01.778150.183988
25MgIIb2797.861.4972431.42240.07111381.398510.05875231.4972431.42240.07111381.398510.0587523
26MgIIn2802.951.1073212.80320.06440691.056760.08093231.1073212.80320.06440691.056760.0809323
27OIII_31333127.7-0.06550150.860.00.9722030.05----------
28[NeV]33463345.39-0.4559320.350.00.7403630.05----------
29[NeV]34263425.660.08635981.220.00.9585640.05----------
30[OII]37283729.660.1931251.560.00.5211380.05----------
31[NeIII]38693869.770.1398791.380.00.7250940.05----------
32HeI38893891.03-1.096910.080.00.3053510.05----------
33[NeIII]39683968.43-0.3467870.450.00.7259120.05----------
34Hd4102.730.7032915.050.01.269980.050.698975.00.01.278750.0880456
35Hg4346.421.1010612.620.01.307920.051.1003712.60.01.301030.0880457
36[OIII]43644363.85-0.3372420.460.00.4913620.05----------
37FeII_Opt374565.0----------1.3010320.00.01.792390.0435751
38Hbeta4862.681.6647446.210.01.606810.051.6646446.20.01.602060.0569717
39[OIII]49604960.360.5440683.50.00.5854610.05----------
40[OIII]50085008.221.1215613.230.00.7810370.051.1139413.00.00.7781510.0673493
41FeII_Opt485306.0----------1.3424222.00.01.875060.0395906
42HeI_58775877.410.6937274.940.01.370140.050.698975.00.01.361730.0880456
43[OI]63026303.050.06069781.150.00.496930.05----------
44[OI]63656370.460.1335391.360.01.007750.05----------
45[NII]65496551.06-0.3665320.430.00.3443920.05----------
46HAb6564.242.30816203.310.02.293670.0500001----------
47HAn6564.942.15996144.530.01.72640.05000012.29003195.00.01.67210.0334735
48[NII]65856585.640.3053512.020.00.408240.05----------
49[SII]67186718.850.2174841.650.00.3201460.05----------
50[SII]67326733.720.1731861.490.00.4048340.05----------
51HeI70677065.670.4857213.060.01.18270.050.4771213.00.01.176090.073064
52[OII]73217321.270.4014012.520.01.154120.05----------
53OI84468457.51.0293810.70.02.017870.051.0413911.00.02.017030.0595932
54[SIII]90699076.8-0.096910.80.01.278750.05----------
55FeII92029214.00.5440683.50.01.910620.050.5440683.50.01.908490.0624694
56Pae9534.40.8450987.00.01.600970.050.8450987.00.01.602060.0624694
57Pad10042.31.3242821.10.02.207370.051.3222221.00.02.206830.0532277
58HeI1083010829.11.1461314.00.01.55630.051.556336.00.02.064460.0
59Pag10861.82.12775134.20.02.371990.05----------
60Pa_g10941.0----------0.8450987.00.02.037430.0624694
61OI1128711296.40.5185143.30.01.896530.05----------
62Pabeta12821.31.2648218.40.02.108230.051.2552718.00.02.107210.048455
63Paalpha18735.51.103812.70.02.292480.051.1139413.00.02.292260.0673493
64HeI2058020506.80.9344988.60.02.477840.050.9542439.00.02.477120.048455

In [8]:
figure(figsize=(12,7))
ax = subplot(211)
ii = where(v6tov5['wave'] < 2900)[0]
plot(ii,v6tov5['logew_v6'][ii],drawstyle='steps-mid')
plot(ii,v6tov5['logew_v5'][ii],'o')
ax.xaxis.set_ticks(ii)
ax.xaxis.set_ticklabels(v6tov5['name'][ii],rotation='vertical',size=11)
ax.grid(axis='y')
xlim(ii[0]-0.5,ii[-1]+0.5)
ax = subplot(212)
ii = where(v6tov5['wave'] > 2900)[0]
plot(ii,v6tov5['logew_v6'][ii].filled(-1),drawstyle='steps-mid')
plot(ii,v6tov5['logew_v5'][ii],'o')
ax.xaxis.set_ticks(ii)
ax.xaxis.set_ticklabels(v6tov5['name'][ii],rotation='vertical',size=11)
xlim(ii[0]-0.5,ii[-1]+0.5)
ax.grid(axis='y')
subplots_adjust(hspace=0.5)



In [9]:
def line_index(trends,line):
    return where(trends['name']==line)[0][0]

In [10]:
v7file = '../simqso/data/emlinetrends_v7.fits'
if not os.path.exists(v7file):
    trendsv7 = trendsv6.copy()
    # replace the bad Pa gamma values in v6 with the v5 values from Glikman et al.
    trendsv7[line_index(trendsv7,'Pag')] = trendsv5[line_index(trendsv5,'Pa_g')]
    trendsv7['name'][line_index(trendsv7,'Pa_g')] = 'Pag' # but keep the name for consistency
    # use a single Halpha component, again falling back to v5 from Glikman
    trendsv7[line_index(trendsv7,'HAb')] = trendsv5[line_index(trendsv5,'Halpha')]
    trendsv7.remove_row(line_index(trendsv7,'HAn'))
    # restore the HeI eq width from v5
    trendsv7[line_index(trendsv7,'HeI10830')] = trendsv5[line_index(trendsv5,'HeI_10830')]
    # change some line EWs based on comparison to Jensen et al.
    trendsv7['logEW'][line_index(trendsv7,'LyAb'),:,1] -= 0.1
    trendsv7['logEW'][line_index(trendsv7,'LyAn'),:,1] -= 0.1
    trendsv7['logEW'][line_index(trendsv7,'NV'),:,1] += 0.2
    trendsv7['logEW'][line_index(trendsv7,'CIVb'),:,1] -= 0.1
    trendsv7['logEW'][line_index(trendsv7,'CIVn'),:,1] -= 0.1
    # add some lya forest lines from Moloney & Shull 2014, default is 20% scatter in EW
    def add_line(name,wave,ew,sigma,ewscatter=0.2):
        trendsv7.insert_row(0,(name,[[0,wave]],
                               [[0,log10(ew)],[0,log10((1-ewscatter)*ew)],[0,log10((1+ewscatter)*ew)]],
                               [[0,log10(sigma)]]))
    add_line('NIII991',991,2.0,5.0)
    add_line('CIII977',985.46,6.55,8.95) # actually from VdB composite
    add_line('Lyepsdel',940.93,2.95,4.73) # ditto
    add_line('OII+OIII833',833,3.14,5)
    add_line('NeVIIIcomplex',773,11.73,10)
    add_line('OIII702',702,2.13,5)
    add_line('NIII686',686,4.16,5)
    add_line('OV629',629,6.55,5)
    # increase LyB+OVI (and rename) to better match MS14
    trendsv7['logEW'][line_index(trendsv7,'LyB'),:,1] += 0.35
    trendsv7['name'][line_index(trendsv7,'LyB')] = 'LyB+OVI'
    trendsv7.write(v7file)
trendsv7 = Table.read(v7file)

In [11]:
emlinesv7 = condense_table(trendsv7)
emlinesv7.show_in_notebook()


Out[11]:
<Table length=69>
idxnamewavelogewewslopesigmascatter
0OV629629.00.8162416.550.00.698970.0880456
1NIII686686.00.6190934.160.00.698970.0880456
2OIII702702.00.328382.130.00.698970.0880456
3NeVIIIcomple773.01.069311.730.01.00.0880456
4OII+OIII833833.00.496933.140.00.698970.0880456
5Lyepsdel940.930.4698222.950.00.6748610.0880456
6CIII977985.460.8162416.550.00.9518230.0880456
7NIII991991.00.301032.00.00.698970.0880456
8LyB+OVI1033.831.3466122.21330.04881390.8108570.0671692
9ArI1068.970.9269248.45130.02931741.172680.0490664
10FeIII:UV11122.010.6553324.522010.03675950.9682470.149757
11CIII*1174.760.6055634.032390.01.034890.240438
12LyAn1216.311.1116512.93140.07561950.6249280.0960547
13LyAb1216.771.748356.01450.09844441.090650.142372
14NV1239.281.3265721.21130.00.9419680.0580401
15SiII1261.670.8703267.418670.00.9484250.0445275
16OI1304.120.5812913.813210.00.9617980.0898251
17CII1337.910.5163373.28350.02522751.04430.111941
18SiIV+OIV]1398.691.1279513.42620.01751861.132860.033169
19L14801487.870.1758081.499020.01.021520.0582295
20CIVn1546.550.9804899.560680.1140110.7479550.16384
21CIVb1548.251.6565745.34940.09233571.311520.0581026
22HeII1636.230.8362276.858460.07143720.9153790.0331512
23OIII]1666.390.8630687.295720.04039011.109650.0284628
24L16901691.30.5655093.677130.2006240.9534210.0191826
25NIII]1746.340.6157494.128080.01584370.9705290.029174
26SiII_18181813.20.5578463.612820.01.143080.0661161
27AlIII1861.660.7779815.997650.01.15880.0280831
28SiIII]1892.7-0.03357930.9255940.1030080.6522250.0535518
29CIII]b1904.121.3497722.37540.020641.259960.0148231
30CIII]n1906.870.2701471.862720.00.7226770.0640038
31fe21202120.00.2304491.70.01.431360.460409
32fe22202220.00.4771213.00.01.778150.183988
33MgIIb2797.861.4972431.42240.07111381.398510.0587523
34MgIIn2802.951.1073212.80320.06440691.056760.0809323
35OIII_31333127.7-0.06550150.860.00.9722030.05
36[NeV]33463345.39-0.4559320.350.00.7403630.05
37[NeV]34263425.660.08635981.220.00.9585640.05
38[OII]37283729.660.1931251.560.00.5211380.05
39[NeIII]38693869.770.1398791.380.00.7250940.05
40HeI38893891.03-1.096910.080.00.3053510.05
41[NeIII]39683968.43-0.3467870.450.00.7259120.05
42Hd4102.730.7032915.050.01.269980.05
43Hg4346.421.1010612.620.01.307920.05
44[OIII]43644363.85-0.3372420.460.00.4913620.05
45Hbeta4862.681.6647446.210.01.606810.05
46[OIII]49604960.360.5440683.50.00.5854610.05
47[OIII]50085008.221.1215613.230.00.7810370.05
48HeI_58775877.410.6937274.940.01.370140.05
49[OI]63026303.050.06069781.150.00.496930.05
50[OI]63656370.460.1335391.360.01.007750.05
51[NII]65496551.06-0.3665320.430.00.3443920.05
52[NII]65856585.640.3053512.020.00.408240.05
53Halpha6565.02.29003195.00.01.67210.0334735
54[SII]67186718.850.2174841.650.00.3201460.05
55[SII]67326733.720.1731861.490.00.4048340.05
56HeI70677065.670.4857213.060.01.18270.05
57[OII]73217321.270.4014012.520.01.154120.05
58OI84468457.51.0293810.70.02.017870.05
59[SIII]90699076.8-0.096910.80.01.278750.05
60FeII92029214.00.5440683.50.01.910620.05
61Pae9534.40.8450987.00.01.600970.05
62Pad10042.31.3242821.10.02.207370.05
63HeI_1083010830.01.556336.00.02.064460.0
64Pag10941.00.8450987.00.02.037430.0624694
65OI1128711296.40.5185143.30.01.896530.05
66Pabeta12821.31.2648218.40.02.108230.05
67Paalpha18735.51.103812.70.02.292480.05
68HeI2058020506.80.9344988.60.02.477840.05

In [12]:
vdblines = generateVdBCompositeEmLines(minEW=0)
vdblines = vdblines(1).squeeze()

In [13]:
def plot_templates():
    figure(figsize=(12,8))
    nsplit = 4
    axes = [ subplot(nsplit,1,n) for n in range(1,nsplit+1) ]
    wave = logspace(log10(700),log10(2.11e4),1200)
    chunks = array_split(arange(len(wave)),nsplit)
    for i in range(len(chunks)-1):
        chunks[i] = np.concatenate([chunks[i],chunks[i+1][:20]])
    for name,emlines,pkwargs in [('VdB',vdblines,dict(c='k',ls='--',lw=2)),
                                 ('v5',emlinesv5,dict(c='C2')),
                                 ('v6',emlinesv6,dict(c='C1')),
                                 ('v7',emlinesv7,dict(c='C0'))]:
        if name == 'VdB':
            linepar = emlines
        else:
            linepar = array([emlines['wave'],emlines['ew'],10**emlines['sigma']]).T
        spec = render_gaussians(wave,0,linepar)
        for i,ii in enumerate(chunks):
            axes[i].plot(wave[ii],1+spec[ii],label=name,**pkwargs)
    for i,ax in enumerate(axes):
        ax.set_xlim(wave[chunks[i][0]],wave[chunks[i][-1]])
        ax.set_yscale('log')
        ax.xaxis.set_major_locator(ticker.LinearLocator(12))
        ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())
    legend()

In [14]:
plot_templates()


Comparison to Jensen et al. composite spectra


In [15]:
# download from here http://www.sdss.org/dr14/data_access/value-added-catalogs/?vac_id=composite-spectra-of-quasars-binned-on-spectroscopic-parameters-from-dr12q/
comps = Table.read('QSO_comps.fits')
# these are parameters for a conversion from Lbol vs. M1350, using BC(1350)=4 as in Jensen et al.
Lbolconv = array([-2.5,91.4662])
Lbolbins = array([45.93,46.25,46.63]) # roughly the average low,mid,hi Lbol values from Jensen+
M1350bins = polyval(Lbolconv,Lbolbins)

In [16]:
from simqso.sqbase import fixed_R_dispersion
from simqso.sqrun import buildQsoSpectrum
from simqso.sqmodels import BossDr9_FeScalings
from simqso.sqgrids import *
wave = fixed_R_dispersion(800,1e4,500)
M = AbsMagVar(FixedSampler(M1350bins),restWave=1350)
z = RedshiftVar(FixedSampler([0.0]*len(M1350bins)))
qsos = QsoSimPoints([M,z],cosmo=None,units='luminosity')
contVar = BrokenPowerLawContinuumVar([ConstSampler(-1.5),ConstSampler(-0.5)],[1190.])
emLineVar = generateBEffEmissionLines(qsos.absMag,NoScatter=True,
                                      EmissionLineTrendFilename='emlinetrends_v7',verbose=10)
feVar = FeTemplateVar(VW01FeTemplateGrid(qsos.z,wave,scales=BossDr9_FeScalings))
qsos.addVars([contVar,emLineVar,feVar])
spec = [ buildQsoSpectrum(wave,None,qsos.getVars(SpectralFeatureVar),dat) for dat in qsos.data ]


loading emission line template emlinetrends_v7

In [17]:
def norm_at(wave,spec,normwave):
    i1,i2 = searchsorted(wave,[normwave-10,normwave+10])
    return spec/median(spec[i1:i2])

In [18]:
def compare_jensen(spec,comps,compname):
    sqspec = norm_at(spec.wave,spec.f_lambda,1450)
    cspec = norm_at(comps['BASEWAVE'],comps[compname],1450)
    figure(figsize=(12,7))
    waves = [(950,1750),(1750,3000)]
    for pnum,(w1,w2) in enumerate(waves,start=1):
        subplot(2,1,pnum)
        ii = where( (spec.wave>w1) & (spec.wave<w2) )[0]
        plot(spec.wave[ii],sqspec[ii],label='simqso')
        ii = where( (comps['BASEWAVE']>w1) & (comps['BASEWAVE']<w2) )[0]
        plot(comps['BASEWAVE'][ii],cspec[ii],label='Jensen+')
        xlim(w1,w2)
        ylim(-0.1,ylim()[1])
    legend()

In [19]:
compare_jensen(spec[0],comps,'COMP11')
gcf().axes[0].set_title('Lbol=45.95,  M1350=%.2f'%M1350bins[0]);



In [20]:
compare_jensen(spec[1],comps,'COMP14')
gcf().axes[0].set_title('Lbol=46.28,  M1350=%.2f'%M1350bins[1]);



In [21]:
compare_jensen(spec[2],comps,'COMP17')
gcf().axes[0].set_title('Lbol=46.66,  M1350=%.2f'%M1350bins[2]);



In [22]:
loglbol = array([45.95,46.28,46.66])
logcivew = array([1.894,1.730,1.574])
m1350 = polyval(Lbolconv,loglbol)
def compare_CIV_EW(trends,emlines):
    li = line_index(trends,'CIVb')
    li2 = line_index(trends,'CIVn')
    m1350sq = linspace(-27,-22,10)
    M1450 = m1350sq
    M_i = M1450 - 1.486 + 0.596
    logcivewsq = emlines['logew'][li] + emlines['slope'][li]*(M_i+25)
    logcivew2sq = emlines['logew'][li2] + emlines['slope'][li2]*(M_i+25)
    logtotcivewsq = log10(np.sum(np.power(10,[logcivewsq,logcivew2sq]),axis=0))
    plot(m1350,logcivew,'s--',label='Jensen+')
    plot(m1350sq,logcivewsq,'o--',label='simqso CIVb')
    plot(m1350sq,logtotcivewsq,'o--',label='simqso CIVtot')
    xlim(-22,-27)
    legend();
figure(figsize=(12,4.5))
subplot(121)
compare_CIV_EW(trendsv6,emlinesv6)
title('version 6')
subplot(122)
compare_CIV_EW(trendsv7,emlinesv7)
title('version 7');



In [ ]: