In [2]:

    

import absorberspec, ebossspec, sdssspec, datapath, fitsio, starburstspec
import cookb_signalsmooth as cbs
from scipy.stats import nanmean, nanmedian
import cosmology as cosmo

# Read in
masterwave, allflux, allivar = ebossspec.rest_allspec_readin()
objs_ori = ebossspec.elg_readin()
nobj = objs_ori.size
median_sn = zeros(objs_ori.size)
galaxytype = objs_ori['CLASS']
zgood = objs_ori['zGOOD']
z = objs_ori['Z']

# Select wavelength range
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]]
tmploglam = log10(tmpwave)

# Calculate median S/N
for i in arange(nobj):
    iuse = (where(tmpivar[:,i]>0))[0]
    median_sn[i] = median(tmpflux[iuse,i]*sqrt(tmpivar[iuse,i]))
    
# Median Composite of All sources with zgood==1
tmpmedian = zeros(tmpwave.size)
for i in arange((tmpflux.shape)[0]):
    iuse = (where(np.logical_and(np.logical_and(tmpivar[i,:]>0, zgood==1), median_sn>0)))[0]
    tmpmedian[i] = median(tmpflux[i,iuse])


    

    




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






    




/Users/Benjamin/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/_methods.py:59: RuntimeWarning: Mean of empty slice.
  warnings.warn("Mean of empty slice.", RuntimeWarning)

In [6]:

    

inuv = (where(logical_and(logical_and(logical_and(zgood, 
                                                   z>0.6), z<1.2), galaxytype=='GALAXY')))[0]
usemedian = zeros(tmpwave.size)
for i in arange((tmpflux.shape)[0]):
    iuse = (where(np.logical_and(tmpivar[i,inuv]>0, median_sn[inuv]>0)))[0]
    usemedian[i] = median(tmpflux[i,inuv[iuse]])


    

In [16]:

    

fig = figure(figsize=(20,7))
ax = fig.add_subplot(111)
ax.plot(tmpwave, usemedian)
ax.set_xlim(2200, 4700)
ax.set_ylim(0.,1)
ax.plot([2796.53, 2796.53], [0,1])


    

    
Out[16]:





[<matplotlib.lines.Line2D at 0x22edc0b10>]






    

In [188]:

    

# Calculate oii_ew and oii_luminosity, only for INUV objects
# 1Mpc = 3.08568025x10^24 cm
Mpc_cm = 3.08568025E24
oiilum = zeros(nobj)
oii_ew = zeros(nobj)
index_oii = searchsorted(tmpwave, 3728.)
dnoiiwave = 10
dwave = median(tmpwave[index_oii-dnoiiwave:index_oii+dnoiiwave]-tmpwave[index_oii-dnoiiwave-1:index_oii+dnoiiwave-1])
print dwave, tmpwave[index_oii-dnoiiwave], tmpwave[index_oii+dnoiiwave]
oiisum = sum(tmpflux[index_oii-dnoiiwave:index_oii+dnoiiwave, inuv]*
             (tmpivar[index_oii-dnoiiwave:index_oii+dnoiiwave, inuv]>0), axis=0)*dwave # Need to subtract the negligible continuum
oiilum[inuv] = oiisum*power(cosmo.luminosity_distance(z[inuv]), 2)*4.*pi*power(Mpc_cm,2)*1E-17
oii_left = sum(tmpflux[index_oii-25:index_oii-15, inuv]*(tmpivar[index_oii-25:index_oii-15, inuv]>0), axis=0)
oii_right = sum(tmpflux[index_oii+15:index_oii+25, inuv]*(tmpivar[index_oii+15:index_oii+25, inuv]>0), axis=0)
oii_ew[inuv] = oiisum/(oii_left+oii_right)*2. #This is just wrong!

# Calculate MgII continum for normalization
index_mgii = (where(logical_or(logical_and(tmpwave>2770., tmpwave<2788.), logical_and(tmpwave>2812., tmpwave<2830.))))[0]
mgii_cont = median(tmpflux[index_mgii,:], axis=0)
tmpflux_mgii = tmpflux/mgii_cont.reshape(1, nobj)

index_2600 = (where(logical_or(logical_and(tmpwave>2554., tmpwave<2572.), logical_and(tmpwave>2630., tmpwave<2648.))))[0]
feii2600_cont = median(tmpflux[index_2600,:], axis=0)
tmpflux_2600 = tmpflux/feii2600_cont.reshape(1, nobj)

index_2400 = (where(logical_or(logical_and(tmpwave>2302., tmpwave<2320.), logical_and(tmpwave>2400., tmpwave<2420.))))[0]
feii2400_cont = median(tmpflux[index_2400,:], axis=0)
tmpflux_2400 = tmpflux/feii2400_cont.reshape(1, nobj)


    

    




0.858309079237 3719.87215282 3737.04230254

In [131]:

    

fig = figure(figsize=(9,7))
ax = fig.add_subplot(111)
n, bins, patches = ax.hist(oii_ew[inuv], 
                           bins=arange(0, 35, 1), edgecolor='black')
this_ylim = [0, 1500]
ax.set_ylim(this_ylim)
#ax.set_xscale('log')
#ax.set_xlim(3680, 3780)
ax.set_title(r'[OII] EW Distribution', fontsize=15)
ax.set_xlabel(r'$W_{\rm O\,II}^{\lambda3727}$ [$\AA$]', fontsize=20)
ax.set_ylabel(r'$\Delta$ N ($\Delta W_{\rm O\,II}^{\lambda3727}$ = 1 $\AA$)', fontsize=20)

# Luminosity
fig = figure(figsize=(9,7))
ax = fig.add_subplot(111)
n, bins, patches = ax.hist(log10(oiilum[inuv]), 
                           bins=arange(40., 43.5, 0.1), edgecolor='black')
this_ylim = [0, 1100]
ax.set_ylim(this_ylim)
#ax.set_xscale('log')
#ax.set_xlim(3680, 3780)
ax.set_title(r'[OII] Lum Distribution', fontsize=15)
ax.set_xlabel(r'$W_{\rm O\,II}^{\lambda3727}$ [$\AA$]', fontsize=20)
ax.set_ylabel(r'$\Delta$ N ($\Delta W_{\rm O\,II}^{\lambda3727}$ = 1 $\AA$)', fontsize=20)


    

    
Out[131]:





<matplotlib.text.Text at 0x113fbe5d0>






    













    

In [227]:

    

#oiiew_logbins = np.arange(log10(4.), log10(100.), log10(1.5))
#oiiew_logbins[0] = log10(2.)
#oiiew_logbins[-1] = 4.
#oiilum_logbins = np.arange(41., 42.5, log10(1.5))
#oiilum_logbins[0] = 40.5
#oiilum_logbins[-1] = 43.
#oiiew_logbins = array([log10(1.), log10(2.8), log10(5.2), log10(10.2), log10(200.)])
#oiilum_logbins = array([40., 41.2, 41.6, 42.08, 44.])
oiiew_logbins = array([log10(1.), log10(2.), log10(4.), log10(8.), log10(16.), log10(200.)])
oiilum_logbins = array([40., 40.8, 41.4, 41.8, 42.4, 44.])
median_oii_ew = zeros(oiiew_logbins.size-1)
for i in np.arange(oiiew_logbins.size-1):
    oiitmp = oii_ew[(oii_ew>power(10., oiiew_logbins[i])) & (oii_ew<power(10., oiiew_logbins[i+1]))]
    median_oii_ew[i] = median(oiitmp)
    print oiitmp.shape, median(oiitmp)
median_oii_lum = zeros(oiilum_logbins.size-1)
for i in np.arange(oiilum_logbins.size-1):
    oiitmp = oiilum[(oiilum>power(10., oiilum_logbins[i])) & (oiilum<power(10., oiilum_logbins[i+1]))]
    median_oii_lum[i] = median(oiitmp)
    print oiitmp.shape, median(oiitmp)


    

    




(150,) 1.77546598331
(2475,) 3.14377383494
(3915,) 5.53537009896
(1623,) 9.94432399986
(352,) 21.4187182271
(179,) 4.33465663019e+40
(2064,) 1.7144831496e+41
(3622,) 3.97922045352e+41
(2576,) 9.67100290622e+41
(171,) 3.08728489513e+42

In [118]:

    




    

    




(2165,) 3.04571087249
(2227,) 4.58434676408
(2254,) 6.71715501795
(1719,) 11.117764339
(2238,) 1.98247943113e+41
(2314,) 3.75295735842e+41
(2267,) 6.80719617252e+41
(1365,) 1.44976813375e+42

In [228]:

    

mgii_oiiew_median = zeros((tmpwave.size, oiiew_logbins.size-1))
mgii_oiilum_median = zeros((tmpwave.size, oiilum_logbins.size-1))
feii2600_oiiew_median = zeros((tmpwave.size, oiiew_logbins.size-1))
feii2600_oiilum_median = zeros((tmpwave.size, oiilum_logbins.size-1))
feii2400_oiiew_median = zeros((tmpwave.size, oiiew_logbins.size-1))
feii2400_oiilum_median = zeros((tmpwave.size, oiilum_logbins.size-1))


    

In [229]:

    

# Note tmpflux includes all objects but for those inapproriate oiiew=0, oiilum=0
for i in arange(oiiew_logbins.size-1):
    iuse = (where(logical_and(oii_ew>power(10., oiiew_logbins[i]), oii_ew<power(10., oiiew_logbins[i+1]))))[0]
    mgii_oiiew_median[:,i] = median(tmpflux_mgii[:,iuse], axis=1)
    feii2600_oiiew_median[:,i] = median(tmpflux_2600[:,iuse], axis=1)
    feii2400_oiiew_median[:,i] = median(tmpflux_2400[:,iuse], axis=1)

for i in arange(oiilum_logbins.size-1):
    iuse = (where(logical_and(oiilum>power(10., oiilum_logbins[i]), oiilum<power(10., oiilum_logbins[i+1]))))[0]
    mgii_oiilum_median[:,i] = median(tmpflux_mgii[:,iuse], axis=1)
    feii2600_oiilum_median[:,i] = median(tmpflux_2600[:,iuse], axis=1)
    feii2400_oiilum_median[:,i] = median(tmpflux_2400[:,iuse], axis=1)


    

In [260]:

    

tmpvel = zeros(tmpwave.size)
tmpvel[1:] = (tmpwave[1:]-2803.53)/tmpwave[1:]*3E5
tmpvel[0] = tmpvel[1]
fig = figure(figsize=(20,8))
#this_ylim = [0.35, 1.45]
this_ylim = [0.1, 2.0]
#this_ylim = [0.2, 12.5]
ax = fig.add_subplot(111)
dwave = median(tmpwave[1:]-tmpwave[:-1])
colors = ['red', 'brown', 'green','blue', 'black']
#colors = ['black', 'violet', 'blue', 'green', 'yellow', 'brown', 'magenta', 'red']
#for i in arange(oiiew_logbins.size-1):
#    ax.plot(tmpwave+dwave/2., oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
for i in arange(oiiew_logbins.size-2)+1:
    #ax.plot(tmpwave+dwave/2., mgii_oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
    ax.plot(tmpvel+69./2., mgii_oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
ax.set_ylim(this_ylim)
#ax.set_xlim(3680, 3780)
ax.set_xlim(-2000., 1000.)
#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'Normalized $f(\lambda)$', fontsize=20)
vel2796 = -(2803.53-2796.35)/2800.*3E5
ax.plot([0,0], this_ylim, '--k')
ax.plot([vel2796, vel2796], this_ylim, '--k')


    

    
Out[260]:





[<matplotlib.lines.Line2D at 0x12be9f610>]






    

In [250]:

    

tmpvel = zeros(tmpwave.size)
tmpvel[1:] = (tmpwave[1:]-2600.17)/tmpwave[1:]*3E5
tmpvel[0] = tmpvel[1]
fig = figure(figsize=(20,8))
#this_ylim = [0.35, 1.45]
this_ylim = [0.2, 1.5]
ax = fig.add_subplot(111)
dwave = median(tmpwave[1:]-tmpwave[:-1])
colors = ['red', 'brown', 'green','blue']
#colors = ['black', 'violet', 'blue', 'green', 'yellow', 'brown', 'magenta', 'red']
#for i in arange(oiiew_logbins.size-1):
#    ax.plot(tmpwave+dwave/2., oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
#for i in arange(oiilum_logbins.size-1):
#    ax.plot(tmpwave+dwave/2., feii2600_oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
i = 0
#ax.plot(tmpwave+dwave/2., feii2600_oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
#ax.plot(tmpwave+dwave/2., feii2600_oiiew_median[:,i], color=colors[i], linewidth=3)
ax.plot(tmpvel+69./2., feii2600_oiiew_median[:,i], color=colors[i], linewidth=3, drawstyle='steps')
i = 3
#ax.plot(tmpwave+dwave/2., feii2600_oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
#ax.plot(tmpwave+dwave/2., feii2600_oiiew_median[:,i], color=colors[i], linewidth=3)
ax.plot(tmpvel+69./2., feii2600_oiiew_median[:,i], color=colors[i], linewidth=3, drawstyle='steps')
ax.set_ylim(this_ylim)
#ax.set_xlim(2550, 2650)
ax.set_xlim(-4000, 4000)
#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'Normalized $f(\lambda)$', fontsize=20)
ax.plot([0,0], this_ylim, '--k')


    

    
Out[250]:





[<matplotlib.lines.Line2D at 0x12133d050>]






    

In [251]:

    

fig = figure(figsize=(20,8))
#this_ylim = [0.35, 1.45]
this_ylim = [0.2, 1.5]
ax = fig.add_subplot(111)
dwave = median(tmpwave[1:]-tmpwave[:-1])
colors = ['red', 'brown', 'green','blue']
#colors = ['black', 'violet', 'blue', 'green', 'yellow', 'brown', 'magenta', 'red']
#for i in arange(oiiew_logbins.size-1):
#    ax.plot(tmpwave+dwave/2., oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
#for i in arange(oiilum_logbins.size-1):
#    ax.plot(tmpwave+dwave/2., feii2600_oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
i = 0
ax.plot(tmpwave+dwave/2., feii2400_oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
i = 3
ax.plot(tmpwave+dwave/2., feii2400_oiiew_median[:,i], color=colors[i], drawstyle='steps', linewidth=3)
ax.set_ylim(this_ylim)
#ax.set_xlim(2300, 2420)
ax.set_xlim(2335, 2390)
#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'Normalized $f(\lambda)$', fontsize=20)


    

    
Out[251]:





<matplotlib.text.Text at 0x123fca190>






    

In [216]:

    

(2803.53-2796.35)/2800.*3E5


    

    
Out[216]:





769.2857142857455

In [214]:

    

iwave = (where(logical_and(tmpwave>2780., tmpwave<2810.)))[0]
dwave = median(tmpwave[iwave[1]:iwave[-1]]-tmpwave[iwave[0]:iwave[-2]])


    

In [215]:

    

69.*(2803.53-2796.35)/dwave


    

    
Out[215]:





769.8920382429967

In [222]:

    

feii2600_oiiew_median.shape


    

    
Out[222]:





(5786, 5)

In [261]:

    

inuv.shape


    

    
Out[261]:





(8620,)

In [ ]: