In [66]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [67]:

    

import sys; sys.path.insert(0, '/global/homes/e/elliek/imaginglss/imaginglss')
import sys; sys.path.insert(0, '/global/homes/e/elliek/imaginglss/fitsio')
import imaginglss
CONF = '/global/project/projectdirs/m779/imaginglss/dr2.conf.py'
!cat $CONF
decals = imaginglss.DECALS(CONF)
DR = decals.datarelease
CAT = DR.catalogue

import sys; sys.path.insert(0, '/global/homes/e/elliek/kdcount')
from kdcount import KDTree, KDAttr

import numpy
from numpy.testing import assert_equal, run_module_suite
from kdcount.utils import constant_array

from imaginglss.model.catalogue import C


    

    




decals_root = '/global/project/projectdirs/cosmo/data/legacysurvey/dr2'
decals_cache = '/global/project/projectdirs/m779/imaginglss/cache'
decals_release = "DR2"
tycho_dir = os.path.abspath(os.path.join(__path__, 'tycho.fit'))
dust_dir = os.path.abspath(os.path.join(__path__, 'SFD98')) 

In [68]:

    

def get_completeness(obj_type):
        
    # Get fluxes from catalog
    pathname = '/global/homes/e/elliek/imaginglss/imaginglss/nersc/'
    filename = pathname + obj_type + '.txt'
    gflux, rflux, zflux, comp = numpy.loadtxt(filename, usecols=(4,5,7,15), unpack=True)
    colors = numpy.concatenate([gflux.reshape(-1, 1), rflux.reshape(-1,1), zflux.reshape(-1, 1)], axis=-1)
    # We just want the 100% complete part to build our estimator
    colors = colors[comp == 1]
    
    ## Later, the (ra,dec,depths) will be fed directly in but for now...
    ra, dec = numpy.loadtxt(filename, usecols=(0,1), unpack=True)
    gdepth, rdepth, zdepth = numpy.loadtxt(filename, usecols=(10,11,13), unpack=True)
    ra = ra[comp == 0]
    dec = dec[comp == 0]
    gdepth = gdepth[comp == 0]
    rdepth = rdepth[comp == 0]
    zdepth = zdepth[comp == 0]
    gdepth = 5.0 * gdepth
    rdepth = 5.0 * rdepth
    zdepth = 3.0 * zdepth
    band_depths = numpy.concatenate([gdepth.reshape(-1, 1), rdepth.reshape(-1,1), zdepth.reshape(-1, 1)], axis=-1)

    # Select bands depending on object type
    # 0,1,2 = g,r,z
    if obj_type == 'LRG':
        mybands = [1, 2]
    elif obj_type == 'ELG':
        mybands = [0, 1, 2]
    elif obj_type == 'QSO':
        mybands = [0, 1]
    else:
        print 'Unknown object. Options: LRG, ELG, QSO. --CHANGE THIS TO ACTUALLY THROW EXCEPTION LATER.'
    # Calculate min/max in each band
    npts = len(ra)
    band_min = numpy.zeros((npts,3)) - numpy.inf
    band_max = numpy.zeros((npts,3)) + numpy.inf
    for band in mybands:
        band_min[:,band] = numpy.min(colors[:,band])
    # Generate tree
    tree = KDTree(colors[:,mybands])
    root = tree.root
    # Ensure fcomplete does not exceed 1 where depth is more than sufficient
    band_depths[band_depths < band_min] = band_min[band_depths < band_min]
    # Search KDTree and calculate fcomplete
    seen = root.integrate((band_depths[:,mybands]),(band_max[:,mybands]))
    total = root.integrate((band_min[:,mybands]),(band_max[:,mybands]))
    # Calculate completeness fraction
    fcomplete = seen*1.0 / total*1.0
    
    # Write to output.
    myoutput = numpy.concatenate([ra.reshape(-1, 1), dec.reshape(-1, 1), fcomplete.reshape(-1, 1)], axis=-1)
    output_name = obj_type+'_completeness.txt'
    numpy.savetxt(output_name, myoutput)


    

In [ ]:

    

get_completeness('LRG')


    

In [69]:

    

get_completeness('ELG')


    

In [65]:

    

get_completeness('QSO')