Spectral analysis of SN spectra

Look through the CC and Ia spectra to:
a) Identify emission lines and features. Annotate the features and make figure.
b) Determine the spectroscopic redshift and uncertainty of this.
c) Measure the equivalent width and fluxes of important spectral lines.
d) Calculate SFR based on available spectral lines, [OII] or H$\alpha$

Start with a good SNR galaxy, SN56.

In [118]:

import pylab as p
import scipy as sp
import pyfits as pyf
import argparse
import matplotlib.cm as cm
from matplotlib.ticker import MultipleLocator
from stsci.convolve import boxcar
from astroconv import z2dL, pc, SFR2LHa
from math import pi

LHa2SFR=1./SFR2LHa #Kennicutt 1998
LOII2SFR=6.58e-42 #Kewley 2004
lam0=p.array([3726.03,4861.33,5006.84,6562.8,6583.41])
linelist=[r'OII-3726',r'H$\beta$-4861',r'OIII-5007',r'H$\alpha$-6563',r'NII-6583']

def plotspec(ax,spectrum,atmlines,lam0,linelist,z=0.):
    """Plot spectrum in an axis"""
    minortick_wl = MultipleLocator(100)
    minortick_sp = MultipleLocator(0.5)
    ax.plot(spectrum[:,0],spectrum[:,1],'k-')
    # set upper and lower limits for error
    llim=spectrum[:,1]-spectrum[:,2]
    ulim=spectrum[:,1]+spectrum[:,2]
    ax.fill_between(spectrum[:,0],ulim,llim,color='r',alpha=.3)
    # plot encircled + signs at the position of possible sky contamination
    ax.plot(atmlines,(spectrum[:,1].mean()-0.8*spectrum[:,1].std())*p.ones(atmlines.size),'g+',ms=10,mfc='None')
    ax.plot(atmlines,(spectrum[:,1].mean()-0.8*spectrum[:,1].std())*p.ones(atmlines.size),'go',ms=10,mfc='None')
    # plot dashed lines at the location of spectral lines, shifted by 1+z
    ylimmax=spectrum[:,1].mean()+(spectrum[:,1].max()-spectrum[:,1].mean())*1.5
    for ii,lam in enumerate(lam0):
        ax.axvline(lam*(1.+z),color='#2AFA05',ls='--',lw=2)
        ypos=ylimmax*0.95
        if ii%2==1: ypos=ylimmax*0.9
        ax.text(lam*(1.+z)-200,ypos,linelist[ii],color='k',size=12)
    # add label with SN info and redshift
    ax.text(6000.,5.,'SN56, type Ibc \n z='+str(z),size=15)
    ax.set_ylim(spectrum[:,1].mean()-1.*spectrum[:,1].std(),ylimmax)
    ax.set_xlim(5000,9999)
    ax.xaxis.set_minor_locator(minortick_wl)
    ax.yaxis.set_minor_locator(minortick_sp)
    ax.set_xlabel(r'$\lambda$ (\AA)',size=15)
    ax.set_ylabel(r'F ($10^{-18}$ erg/cm$^2$/s/\AA)',size=15)
    return

def smoothspec(spectrum,bc):
    smspec=p.array([spectrum[:,0],boxcar(spectrum[:,1],(bc,),
			mode="nearest"),boxcar(spectrum[:,2],(bc,),
			mode="nearest")]).transpose()
    return smspec

def getlum(flux,zz,ext,h=0.7,om_m=0.3):
    """Get luminosity in erg/s (if input is in erg/s/cm2)"""
    #use cosmology calcs from astroconv, convert to cm
    distance=z2dL(zz,h*100.,om_m,(1.-om_m))*pc
    return flux/ext*4*pi*distance**2

def extCCM(lam,Rv,EBV):
    #first get k(x), NB only valid for optical wavelengths
    x=1./(lam/1e-6)
    y=x-1.82
    k=((1+0.17699*y-0.50447*y**2-0.02427*y**3+0.72085*y**4+0.01979*y**5-0.77530*y**6+0.32999*y**7)*Rv+
        1.41338*y+2.28305*y**2+1.07233*y**3-5.38434*y**4-0.62551*y**5+5.30260*y**6-2.09002*y**7)
    return k,10**(0.4*k*EBV)
 

atmlines=p.array([5577.,5896.,6300.,7650.,8400,8850.,9400.,9850.])

In [119]:

sn56sp=p.loadtxt('CC/spec_SNH_56.dat')

In [120]:

#The mean error is used with splot to get an error estimate for the line centers and EQWs
print sn56sp[:,2].mean()

0.00134045841012

In [121]:

# Get line centers, EQWs, and fluxes (+errors) for OII, Hb, OIII, Ha, NII from splot
# OII, Hb, OIII, Ha, NII
linec=p.array([5568.6,7260.29,7479.79,9801.826,9833.95])
dlinec=p.array([0.19678,0.10547,0.37,0.15332,0.25])
eqws=p.array([17.88,8.367,1.665,68.59,40.0])
deqws=p.array([0.5337,0.1628,0.1974,2.257,2.599])
fluxes=p.array([0.352228,0.3411,0.0705263,0.885772,0.424059])*1e-16
dflux=p.array([0.01051,0.00664,0.00836,0.01255,0.01279])*1e-16

In [122]:

zsp=linec/lam0 - 1.
print zsp.mean()
print zsp.std()
#Ha/Hb=2.86 or not? Higher indicates extinction, but stellar absorption in Hb!!!
print fluxes[3]/fluxes[1]

0.493839168037
0.000370392667729
2.59681031955

Now calculate the SFR from Ha and OII. Assume CCM extinction with Rv=3.1 and E(B-V)=0.7 (from photoz fitting). Kewley et al 2004

In [123]:

k_OII,ext_OII=extCCM(3726e-10,3.1,0.7)
k_Ha,ext_Ha=extCCM(6563e-10,3.1,0.7)

In [124]:

SFR56_Ha=getlum(fluxes[3],0.49,1./ext_Ha,0.7,0.3)*LHa2SFR
SFR56_OII=getlum(fluxes[0],0.49,1./ext_OII,0.7,0.3)*LOII2SFR
print SFR56_Ha, SFR56_OII

3.28386031244 4.59406840266

In [125]:

# convert flux units to 10^-18 erg/s/Å
sn56sp[:,1]=100.*sn56sp[:,1]
sn56sp[:,2]=100.*sn56sp[:,2]
# and smooth spectra with a boxcar filter
smspec=smoothspec(sn56sp,10)
spfig=p.figure(figsize=(12,6))
spax=spfig.add_subplot(111)
plotspec(spax,smspec,atmlines,lam0,linelist,0.4938)
p.savefig('SN56_specfinal.png',dpi=200)
p.show()

Spectral analysis of SN spectra

CC SNe

SN56