Simple tests for QSO templates

A simple notebook that plays around with the templates.QSO Class.


In [1]:
import numpy as np
import warnings

import matplotlib.pyplot as plt
from matplotlib.patches import Polygon

from desisim.templates import QSO

%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [2]:
seed = 123
nmodel = 10

In [3]:
qso = QSO(minwave=3000, maxwave=5e4)

Make templates with and without the fast Lyman-alpha forest.


In [4]:
flux, wave, meta = qso.make_templates(nmodel=nmodel, zrange=(2.0, 4.0), seed=seed, 
                                      nocolorcuts=True, lyaforest=False)

In [5]:
flux_forest, _, meta_forest = qso.make_templates(nmodel=nmodel, zrange=(2.0, 4.0), seed=seed, 
                                                 nocolorcuts=True, lyaforest=True)

In [6]:
meta


Out[6]:
<Table length=10>
OBJTYPESUBTYPETEMPLATEIDSEEDREDSHIFTMAGDECAM_FLUX [6]WISE_FLUX [2]OIIFLUXHBETAFLUXEWOIIEWHBETAD4000VDISPOIIDOUBLETOIIIHBETAOIIHBETANIIHBETASIIHBETAZMETALAGETEFFLOGGFEH
erg / (cm2 s)erg / (cm2 s)AngstromAngstromkm / sdexdexdexdexGyrKm / s2
str10str10int64int64float64float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32
QSO029913123822.8462129202521.8451.13045 .. 1.934743.5364 .. 25.9223-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO130621197893.9615283967720.45620.0979406 .. 8.0199932.9949 .. 54.5157-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO212289591023.3696594771720.43862.50919 .. 6.380445.58704 .. 33.0771-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO318402686102.9618638029721.32891.82292 .. 3.175543.74577 .. 40.9519-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO49743195802.7842350363921.32961.84122 .. 2.867044.0844 .. 40.5941-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO529673278422.686356032321.5861.62835 .. 3.141947.27288 .. 37.9606-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO623678788863.4580994147722.12360.122572 .. 1.645967.81299 .. 15.1341-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO730887270572.8771444893621.81110.928251 .. 2.7343614.4657 .. 32.5649-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO830900956992.1193557932221.52763.99881 .. 2.3398952.4078 .. 131.512-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSO921093397542.7960885106621.80611.15508 .. 2.36017.46708 .. 35.2731-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0

In [7]:
meta_forest


Out[7]:
<Table length=10>
OBJTYPESUBTYPETEMPLATEIDSEEDREDSHIFTMAGDECAM_FLUX [6]WISE_FLUX [2]OIIFLUXHBETAFLUXEWOIIEWHBETAD4000VDISPOIIDOUBLETOIIIHBETAOIIHBETANIIHBETASIIHBETAZMETALAGETEFFLOGGFEH
erg / (cm2 s)erg / (cm2 s)AngstromAngstromkm / sdexdexdexdexGyrKm / s2
str10str10int64int64float64float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32
QSOLYA029913123822.8462129202521.8451.02201 .. 1.934653.53665 .. 25.8862-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA130621197893.9615283967720.45620.0970256 .. 8.7027535.8276 .. 59.1432-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA212289591023.3696594771720.43862.22059 .. 6.378215.58765 .. 32.5178-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA318402686102.9618638029721.32891.62206 .. 3.175283.7461 .. 40.7308-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA49743195802.7842350363921.32961.68465 .. 2.86694.08463 .. 40.5309-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA529673278422.686356032321.5861.4656 .. 3.141877.2734 .. 37.9339-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA623678788863.4580994147722.12360.104318 .. 1.645667.81413 .. 14.92-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA730887270572.8771444893621.81110.816304 .. 2.7343314.4665 .. 31.7727-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA830900956992.1193557932221.52763.83194 .. 2.3398452.3716 .. 131.417-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0
QSOLYA921093397542.7960885106621.80611.01241 .. 2.360047.46744 .. 28.6596-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0-1.0

Show the forest


In [8]:
for ii in range(nmodel):
    plt.plot(wave, flux_forest[ii, :])
    plt.plot(wave, flux[ii, :])
    plt.xlim(3000, 6300)
    plt.show()


Show the effect of extrapolation


In [18]:
for ii in range(nmodel):
    plt.plot(wave, flux[ii, :])
    #plt.xlim(3000, 200)
    plt.xscale('log')
    plt.show()


Look at the color-cuts.


In [10]:
flux1, _, meta1 = qso.make_templates(nmodel=100, seed=1, lyaforest=True, nocolorcuts=True)
flux2, _, meta2 = qso.make_templates(nmodel=100, seed=1, lyaforest=True, nocolorcuts=False)


WARNING:templates.py:1914:make_templates: 2 spectra could not be computed given the input priors!
WARNING:templates.py:1914:make_templates: 2 spectra could not be computed given the input priors!

In [11]:
fail = np.where(np.sum(flux2, axis=1) == 0)[0]
fail


Out[11]:
array([], dtype=int64)

In [12]:
def qso_colorbox(ax, plottype='grz'):
    """Draw the QSO selection boxes."""
    rmaglim = 22.7
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    if plottype == 'grz-r':
        verts = [(xlim[0]-0.05, 17.0),
                 (22.7, 17.0),
                 (22.7, ylim[1]+0.05),
                 (xlim[0]-0.05, ylim[1]+0.05)
                ]
    if plottype == 'rW1-rz':
        verts = None
        ax.axvline(x=-0.3, ls='--', color='k')
        ax.axvline(x=1.3, ls='--', color='k')

    if plottype == 'gr-rz':
        verts = [(-0.3, 1.3),
                 (1.1, 1.3),
                 (1.1, ylim[0]-0.05),
                 (-0.3, ylim[0]-0.05)
                ]
    if verts:
        ax.add_patch(Polygon(verts, fill=False, ls='--', color='k'))

In [13]:
def flux2colors(cat):
    """Convert DECam/WISE fluxes to magnitudes and colors."""
    colors = dict()
    with warnings.catch_warnings(): # ignore missing fluxes (e.g., for QSOs)
        warnings.simplefilter('ignore')
        for ii, band in zip((1, 2, 4), ('g', 'r', 'z')):
            colors[band] = 22.5 - 2.5 * np.log10(cat['DECAM_FLUX'][..., ii].data)
        colors['grz'] = 22.5-2.5*np.log10((cat['DECAM_FLUX'][..., 1] + 
                                           0.8 * cat['DECAM_FLUX'][..., 2] +
                                           0.5 * cat['DECAM_FLUX'][..., 4]).data / 2.3)

        colors['gr'] = colors['g'] - colors['r']
        colors['rz'] = colors['r'] - colors['z']
    
    return colors

In [14]:
nocuts = flux2colors(meta1)
cuts = flux2colors(meta2)

In [15]:
fig, ax = plt.subplots()
ax.scatter(nocuts['rz'], nocuts['gr'], s=14, label='No Color-cuts')
ax.scatter(cuts['rz'], cuts['gr'], s=14, marker='s', alpha=0.7, label='With Color-cuts')
ax.set_xlabel('$r - z$')
ax.set_ylabel('$g - r$')
ax.set_xlim(-1, 2.2)
ax.set_ylim(-1, 2.0)
ax.legend(loc='upper right')
qso_colorbox(ax, 'gr-rz')



In [ ]: