In [1]:
import absorberspec, ebossspec, sdssspec, datapath, fitsio
import cookb_signalsmooth as cbs
from scipy.stats import nanmean, nanmedian
import cosmology as cosmo

In [2]:
masterwave, allflux, allivar = ebossspec.rest_allspec_readin()
objs_ori = ebossspec.elg_readin()
nobj = objs_ori.size
index_wave = searchsorted(masterwave, [2300., 4001.])


Reading /Users/Benjamin/AstroData/AllInOne/AIO_ELG_eBOSS_SDSSRestFrame_Wave01800_03600A.fits.
Reading /Users/Benjamin/AstroData/AllInOne/AIO_ELG_eBOSS_SDSSRestFrame_Wave03600_07200A.fits.

In [5]:
index_oii = searchsorted(masterwave, 3727.)
ii = (where(logical_and(logical_and(allivar[index_wave[0], :]>0, allivar[index_wave[1], :]>0), objs_ori['zGOOD']==1)))[0]
oiisum = sum(allflux[index_oii-10:index_oii+10, ii], axis=0)
oiilum = oiisum*power(cosmo.luminosity_distance(objs_ori[ii]['Z']), 2)
oii_left = sum(allflux[index_oii-25:index_oii-15, ii], axis=0)
oii_right = sum(allflux[index_oii+15:index_oii+25, ii], axis=0)
oii_ew = oiisum/(oii_left+oii_right)*2.
xx = allflux[index_wave[0]:index_wave[1], ii]

In [74]:
sort_index = argsort(oii_ew)
xmean1 = nanmean(xx[:,sort_index[100:2000]], 1)
xmean2 = nanmean(xx[:,sort_index[1000:3000]], 1)
xmean3 = nanmean(xx[:,sort_index[2000:4000]], 1)
xmean4 = nanmean(xx[:,sort_index[3000:4000]], 1)
xmean5 = nanmean(xx[:,sort_index[4000:-100]], 1)
xmean = nanmean(xx[:,sort_index[2000:-100]], 1)
xmeanall = nanmean(xx,1)

In [1]:
figure(figsize=(10,6))
plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean5, 'b', drawstyle='steps')
xlim(2750, 2850); ylim(0.15,0.5)
plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean1*0.1, 'r', drawstyle='steps')
plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean*0.6, 'g', drawstyle='steps')
plot([2796.75,  2796.75], [0.0, 3.5], 'k')
plot([2803.27,  2803.27], [0.0, 3.5], 'k')


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-b9ee5227f899> in <module>()
      1 figure(figsize=(10,6))
----> 2 plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean5, 'b', drawstyle='steps')
      3 xlim(2750, 2850); ylim(0.15,0.5)
      4 plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean1*0.1, 'r', drawstyle='steps')
      5 plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean*0.6, 'g', drawstyle='steps')

NameError: name 'masterwave' is not defined
<matplotlib.figure.Figure at 0x108100cd0>

In [38]:
figure(figsize=(10,6))
plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean2*0.75, 'r', drawstyle='steps')
#plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean4*1.1, 'g', drawstyle='steps')
plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean5*1.15, 'b', lw=2, drawstyle='steps')
xlim(2780, 2820); ylim(0.25,0.6)
plot([2796.75,  2796.75], [0.0, 3.5], 'k')
plot([2803.27,  2803.27], [0.0, 3.5], 'k')


Out[38]:
[<matplotlib.lines.Line2D at 0x2554b1390>]

In [89]:
figure(figsize=(10,6))
plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean2*0.85, 'r', drawstyle='steps')
#plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean4*1.1, 'g', drawstyle='steps')
plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean5*1.15, 'b', lw=2, drawstyle='steps')
xlim(2580, 2635); ylim(0.15,0.8)
#plot(masterwave[index_wave[0]:index_wave[1]]+0.5, xmean*1.8, 'g', drawstyle='steps')
plot([2586.65, 2586.65], [0.0, 1.25], 'g')
plot([2600.17, 2600.17], [0.0, 1.25], 'g')
plot([2626.45, 2626.45], [0.0, 1.25], 'r')
plot([2612.65, 2612.65], [0.0, 1.25], 'r')
#text(2586.65+0.5, 1.15, 'Fe II', color='green')
#text(2600.17+0.5, 1.15, 'Fe II', color='green')
#text(2626.45+0.5, 0.07, 'Fe II*', color='red')
#text(2612.65+0.5, 0.07, 'Fe II*', color='red')


Out[89]:
[<matplotlib.lines.Line2D at 0x24a5d2dd0>]

In [75]:
figure(figsize=(20,6))
plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean, 'r', drawstyle='steps')
#plot(masterwave[index_wave[0]:index_wave[1]]+0.6, xmean5, 'b', drawstyle='steps')
xlim(2300, 2870); ylim(0.30,0.55)
plot([2586.65, 2586.65], [0.0, 1.25], 'g')
plot([2600.17, 2600.17], [0.0, 1.25], 'g')
plot([2626.45, 2626.45], [0.0, 1.25], 'b')
plot([2612.65, 2612.65], [0.0, 1.25], 'b')
plot([2344.21, 2344.21], [0.05, 1.25], 'g')
plot([2374.46, 2374.46], [0.05, 1.25], 'g')
plot([2382.77, 2382.77], [0.05, 1.25], 'g')
plot([2365.55, 2365.55], [0.05, 1.25], 'b')
plot([2396.35, 2396.35], [0.05, 1.25], 'b')
plot([2796.75,  2796.75], [0.0, 3.5], 'k')
plot([2803.27,  2803.27], [0.0, 3.5], 'k')


Out[75]:
[<matplotlib.lines.Line2D at 0x23d8c0450>]

In [88]:
figure(figsize=(20,6))
plot(masterwave[index_wave[0]:index_wave[1]]+0.55, xmean2*0.82, 'r', drawstyle='steps')
#plot(masterwave[index_wave[0]:index_wave[1]]+0.55, xmean4*1.1, 'g', drawstyle='steps')
plot(masterwave[index_wave[0]:index_wave[1]]+0.55, xmean5*1.0, 'b', lw=2, drawstyle='steps')
xlim(2310, 2410)
ylim(0.25,0.75)
plot([2344.21, 2344.21], [0.25, 1.25], 'g')
plot([2374.46, 2374.46], [0.25, 1.25], 'g')
plot([2382.77, 2382.77], [0.25, 1.25], 'g')
plot([2365.55, 2365.55], [0.25, 1.25], 'r')
plot([2396.35, 2396.35], [0.25, 1.25], 'r')


Out[88]:
[<matplotlib.lines.Line2D at 0x24a594b90>]

In [3]:
index_wave_all = searchsorted(masterwave, [1900., 7200.])
tmpflux = allflux[index_wave_all[0]:index_wave_all[1],:]
tmpivar = allivar[index_wave_all[0]:index_wave_all[1],:]
tmpwave = masterwave[index_wave_all[0]:index_wave_all[1]]

In [4]:
tmpmedian = zeros((tmpflux.shape)[0])
for i in arange((tmpflux.shape)[0]):
    iuse = (where(np.logical_and(tmpivar[i,:]>0, objs_ori['zGOOD']==1)))[0]
    tmpmedian[i] = median(tmpflux[i,iuse])

In [228]:
fig = figure(figsize=(20,8))
ax = fig.add_subplot(111)
ax.plot(tmpwave, tmpmedian)
ax.set_xlim(2000, 7250)
ax.set_ylim(-0.05, 2.0)
ax.set_title('Composite Spectrum of Emission Line Galaxies (ELGs) from eBOSS', fontsize=15)
ax.set_xlabel(r'$\lambda$ [$\AA$]', fontsize=20)
ax.set_ylabel(r'$f(\lambda)$ [arbitrary unit]', fontsize=20)
for axis in ['top','bottom','left','right']:
  ax.spines[axis].set_linewidth(2.0)
ax.tick_params(axis='both', which='both', length=7, width=2, labelsize=15)
#emission_name = array(['FeII*', 'FeII*', '[OII]', '[NeIII]','[OIII]', '[OIII]', 'HeI', '[OI]',
#                       'NII', 'NII', 'SII', 'SII', '[ArIII]'])
#emission_xpos = array([2300.,   2550.,   3727.,   3868.,     4363.,   4959.,    5876., 6300.,
#                       6548., 6583., 6716., 6730., 7137.])
for i in arange(nemlines):
    etmp = emission[i]
    ax.text(etmp['XPOS'], etmp['YPOS'], etmp['NAME'], fontsize=14)
for i in arange(nabslines):
    atmp = absorption[i]
    ax.text(atmp['XPOS'], atmp['YPOS'], atmp['NAME'], fontsize=14)
for i in arange(nbalmer):
    btmp = balmer[i]
    ax.text(btmp['XPOS'], btmp['YPOS'], btmp['NAME'], color='green', fontsize=20)
plot([3934.78, 3934.78], [0.11, 0.17], 'k')
plot([3969.59, 3969.59], [0.11, 0.17], 'k')
plot([5891.58, 5891.58], [0.23, 0.29], 'k')
plot([5897.56, 5897.56], [0.23, 0.29], 'k')


Out[228]:
[<matplotlib.lines.Line2D at 0x241f341d0>]

In [227]:
nemlines = 15
emission = zeros(nemlines, dtype=[('NAME','S10'),('XPOS','f4'),('YPOS','f4'), ('WAVE','f4')])
emission[0] = ('FeII*', 2344.-0., 0.52, 2344.)
emission[1] = ('FeII*', 2600.-25., 0.52, 2600.)
emission[2] = ('[OII]', 3727.+20., 1.75, 3727.)
emission[3] = ('[NeIII]', 3868.-100., 0.63, 3868.)
#emission[4] = ('[OIII]', 4363., 0.8, 4363.)
emission[5] = ('[OIII]', 4959.-100., 0.96, 4959.)
emission[6] = ('[OIII]', 5007.-50., 1.75, 5007.)
emission[7] = ('HeI', 5876.-55., 0.56, 5876.)
emission[8] = ('[OI]', 6302.-70., 0.51, 6302.)
emission[9] = ('[NII]', 6548.-120., 0.57, 6548.)
emission[10] = ('[NII]', 6583.-5., 0.89, 6583.)
emission[11] = ('[SII]', 6716.+20., 0.93, 6716.)
#emission[12] = ('SII', 6730., 1.2, 6730.)
emission[13] = ('[ArIII]', 7137.-100., 0.45, 7137.)
emission[14] = ('CII]', 2327.34-60., 0.60, 2327.34)
 
nabslines = 10
absorption = zeros(nabslines, dtype=[('NAME','S20'),('XPOS','f4'),('YPOS','f4'),('WAVE', 'f4')])
absorption[0] = ('FeII', 2300.+10., 0.20, 2300.)
absorption[1] = ('FeII', 2600.-60., 0.16, 2600.)
absorption[2] = ('MgII', 2800.-60., 0.17, 2800.)
absorption[3] = ('MgI', 2853.-30., 0.24, 2853.)
absorption[4] = ('CaII', 3934.78-180., 0.05, 3934.78)
absorption[5] = ('K', 3934.78-25., 0.05, 3934.78)
absorption[6] = ('H', 3969.59-10., 0.05, 3969.59)
#absorption[7] = ('MgI b', 5100.-35., 0.26, 5100.)
absorption[8] = (r'NaI D$_{2,1}$', 5900.-165., 0.15, 5900.)

nbalmer = 7
balmer = zeros(nbalmer, dtype=[('NAME','S10'),('XPOS','f4'),('YPOS','f4'),('WAVE', 'f4')])
balmer[0] = (r'$\alpha$', 6563.-35., 0.25, 6563.)
balmer[1] = (r'$\beta$', 4861.-35., 0.25, 4861.)
balmer[2] = (r'$\gamma$', 4341.-35., 0.265, 4341.)
balmer[3] = (r'$\delta$', 4102.-35., 0.24, 4102.)
balmer[4] = (r'$\epsilon$', 3970.-30., 0.22, 3907.)
balmer[5] = (r'$\zeta$', 3889.-30., 0.23, 3889.)
balmer[6] = (r'$\eta$', 3835.-30., 0.22, 3835.)

In [229]:
fig = figure(figsize=(20,8))
ax = fig.add_subplot(111)
ax.plot(tmpwave+0.55, tmpmedian, drawstyle='steps')
ax.set_xlim(2200, 2880)
ax.set_ylim(0.15, 0.7)
ax.set_title('Composite Spectrum of Emission Line Galaxies (ELGs) from eBOSS', fontsize=15)
ax.set_xlabel(r'$\lambda$ [$\AA$]', fontsize=20)
ax.set_ylabel(r'$f(\lambda)$ [arbitrary unit]', fontsize=20)
for axis in ['top','bottom','left','right']:
  ax.spines[axis].set_linewidth(2.0)
ax.tick_params(axis='both', which='both', length=7, width=2, labelsize=15)
#plot([2423.46, 2423.46], [0.3, 0.6], 'k')
#plot([2467.95, 2467.95], [0.3, 0.6], 'k')
#plot([2683.88, 2683.88], [0.3, 0.6], 'k')
#plot([2719.83, 2719.83], [0.3, 0.6], 'k')
#plot([2740.52, 2740.52], [0.3, 0.6], 'k')



In [276]:
reload(ebossspec)


Out[276]:
<module 'ebossspec' from 'ebossspec.py'>

In [277]:
(outwave, fluxmean, fluxmedian, norm_fluxmean, norm_fluxmedian) = ebossspec.feiimgii_composite()


Reading /Users/Benjamin/AstroData/AllInOne/AIO_ELG_eBOSS_SDSSRestFrame_Wave01800_03600A.fits.
Reading /Users/Benjamin/AstroData/AllInOne/AIO_ELG_eBOSS_SDSSRestFrame_Wave03600_07200A.fits.

In [282]:
figure(figsize=(20,8))
plot(outwave, norm_fluxmedian)
ylim(0.4,1.6)
xlim(2200, 2900)
plot([2200, 2900], [1,1])


Out[282]:
[<matplotlib.lines.Line2D at 0x2521e4650>]

In [275]:
outloglam = log10(outwave)
if True:
    fluxflag = np.ones(fluxmedian.size)
    wave_pos = np.array([2200.])
    rest_loc = np.searchsorted(outwave, wave_pos)
    fluxflag[0:rest_loc[0]] = 0
    # Fe II 2350
    wave_pos = np.array([2330., 2420])
    rest_loc = np.searchsorted(outwave, wave_pos)
    fluxflag[rest_loc[0]:rest_loc[1]] = 0.
    # Fe II 2600 
    wave_pos = np.array([2570., 2640])
    rest_loc = np.searchsorted(outwave, wave_pos)
    fluxflag[rest_loc[0]:rest_loc[1]] = 0.
    # Mg II 2800
    wave_pos = np.array([2770., 2820])
    rest_loc = np.searchsorted(outwave, wave_pos)
    fluxflag[rest_loc[0]:rest_loc[1]] = 0.
    wave_pos = np.array([2843., 2863])
    rest_loc = np.searchsorted(outwave, wave_pos)
    fluxflag[rest_loc[0]:rest_loc[1]] = 0.
    # right 2900
    wave_pos = np.array([2900.])
    rest_loc = np.searchsorted(outwave, wave_pos)
    fluxflag[rest_loc[0]:] = 0.

    imask = (np.where(fluxflag>0.))[0]
    if imask.size>10:
       x = outloglam[imask]
       # Mean
       y = fluxmedian[imask]
       z = np.polyfit(x, y, 3)
       p = np.poly1d(z)
       continuum = p(outloglam)
       norm_fluxmedian = fluxmedian/continuum
        
plot(outwave[imask], fluxmedian[imask], outwave[imask], continuum[imask])


Out[275]:
[<matplotlib.lines.Line2D at 0x2522ba550>,
 <matplotlib.lines.Line2D at 0x2522ba7d0>]

In [246]:
hist(objs_ori['Z'], bins=arange(0.4, 1, 0.02))


Out[246]:
(array([ 247.,  204.,  206.,  232.,  254.,  307.,  215.,  250.,  307.,
         330.,  458.,  426.,  447.,  470.,  555.,  542.,  600.,  765.,
         719.,  656.,  561.,  594.,  596.,  520.,  483.,  449.,  306.,
         320.,  353.]),
 array([ 0.4 ,  0.42,  0.44,  0.46,  0.48,  0.5 ,  0.52,  0.54,  0.56,
         0.58,  0.6 ,  0.62,  0.64,  0.66,  0.68,  0.7 ,  0.72,  0.74,
         0.76,  0.78,  0.8 ,  0.82,  0.84,  0.86,  0.88,  0.9 ,  0.92,
         0.94,  0.96,  0.98]),
 <a list of 29 Patch objects>)

In [248]:
ivar = np.where(allivar==0.)

In [252]:
ivar[1].shape


Out[252]:
(204714002,)

In [ ]: