In [1]:
%matplotlib inline
from __future__ import division
import om10,os
import numpy as np
import matplotlib.pyplot as plt
import triangle
In [2]:
# redshift sigma g/Reff r i z g mag r mag i mag z w1 w2 w3 w4
lrgs = np.loadtxt(os.path.expandvars('$OM10_DIR/data/LRGo.txt'))
#kill photometric outliers in r and i bands
lrgs = lrgs[np.logical_and(lrgs[:,8]<24,lrgs[:,9]<22)]
#color outliers:
lrgs = lrgs[(lrgs[:,6] - lrgs[:,7]) < 3]
lrgs = lrgs[(lrgs[:,7] - lrgs[:,8]) > 0]
lrgs = lrgs[(lrgs[:,8] - lrgs[:,9]) > 0]
#redshift g r i z w1 w2 w3 w4
qsos = np.loadtxt(os.path.expandvars('$OM10_DIR/data/QSOo.txt'))
qsos = qsos[np.logical_and(qsos[:,3]>0,qsos[:,1]>0)] #require g and i band > 0 to kill -10,000 outliers
In [3]:
_ = plt.hist(lrgs[:,0].flatten(),bins=100)
plt.xlabel('z')
Out[3]:
Well, we don't have that great a handle on redshifts past .8, but hopefully ok...
In [3]:
#load and decorate OM10 catalog
db = om10.DB(catalog='$OM10_DIR/data/qso_mock.fits')
#PS1-like sample!
db.select_random(maglim=21.4,area=100000.0,IQ=1.0)
print 'Total number of sims: ', len(db.sample)
db.get_sky_positions(input_cat='$OM10_DIR/data/CFHTLS_LRGs.txt')
#db.assign_sky_positions()
db.paint(lrg_input_cat='$OM10_DIR/data/LRGo.txt',qso_input_cat='$OM10_DIR/data/QSOo.txt')
#which objects got correctly colored in our desired parameter range?
matched = db.sample[db.sample['ZLENS'] < 3]
#matched = matched[matched['ZLENS'] < 3]
matched = matched[np.logical_and(matched['SDSS_FLAG_SRC'] == 0,matched['SDSS_FLAG_LENS'] == 0)]
matched = matched[matched['ZSRC'] < 4.5]
matched = matched[np.logical_and(matched['VELDISP'] > 200,matched['VELDISP'] < 300)]
print 'Fraction remaining after bad color matching: ', len(matched)/len(db.sample)
#more restrictive range that improves color matching
matched = matched[matched['ZLENS'] < .98]
print 'Fraction remaining after good color matching: ', len(matched)/len(db.sample)
#which objects got assigned a position?
#matched = db.sample[np.logical_and(db.sample['RA']> 0,db.sample['DEC'] > -20)]
#matched = matched[np.logical_and(matched['RA']> 0,matched['DEC'] > -20)]
print 'Fraction remaining after position-matching: ', len(matched)/len(db.sample)
Going from ZLENS < 3 to ZLENS < 1 costs us 9% of our training set, but I think will greatly improve our ability to apply XD to this problem, because otherwise, there will be massive color peaks in the lens color distributions. The fixes Adri suggested are good, but I propose we save them for pass 2...
So the final prediction is that there would be 379 lenses that we could hope to find in our PS1QLS sample!
In [40]:
len(matched) #total number of systems in input cat!
Out[40]:
In [34]:
db.sample = matched
out_cat = db.make_sim_input_catalog()
In [35]:
out_cat
Out[35]:
In [38]:
out_cat.write('input_cat_for_eric.fits',format='fits',overwrite=True)
In [4]:
to_plot = np.array([lrgs[:,0], lrgs[:,1], lrgs[:,8], lrgs[:,6]-lrgs[:,7], lrgs[:,7]-lrgs[:,8],\
lrgs[:,8]-lrgs[:,9], lrgs[:,9]-lrgs[:,10], lrgs[:,10]-lrgs[:,11]]).transpose()
painted_lenses = np.array([matched['ZLENS'],matched['VELDISP'], matched['MAGI_LENS'], matched['MAGG_LENS']-matched['MAGR_LENS'], \
matched['MAGR_LENS']-matched['MAGI_LENS'],matched['MAGI_LENS']-matched['MAGZ_LENS'], \
matched['MAGZ_LENS']-matched['MAGW1_LENS'],matched['MAGW1_LENS']-matched['MAGW2_LENS']]).transpose()
fig = triangle.corner(to_plot,labels=['z','sigma','i mag','g-r','r-i','i-z','z-w1','w1-w2'],color='Blue')
_ = triangle.corner(painted_lenses,labels=['z','sigma','i mag','g-r','r-i','i-z','z-w1','w1-w2'],color='Green',fig=fig)
plt.suptitle('SDSS Lenses')
Out[4]:
In [10]:
to_plot = np.array([qsos[:,0], qsos[:,3], qsos[:,1]-qsos[:,2], qsos[:,2]-qsos[:,3], qsos[:,3]-qsos[:,4],\
qsos[:,4]-qsos[:,5],qsos[:,5]-qsos[:,6]]).transpose()
painted_qsos = np.array([matched['ZSRC'],matched['MAGI'],matched['MAGG_SRC']-matched['MAGR_SRC'], \
matched['MAGR_SRC']-matched['MAGI_SRC'],matched['MAGI_SRC']-matched['MAGZ_SRC'], \
matched['MAGZ_SRC']-matched['MAGW1_SRC'],matched['MAGW1_SRC']-matched['MAGW2_SRC']]).transpose()
fig = triangle.corner(to_plot,labels=['z','magi','g-r','r-i','i-z','z-w1','w1-w2'],color='Blue')
_ = triangle.corner(painted_qsos,labels=['z','magi','g-r','r-i','i-z','z-w1','w1-w2'],color='Green',fig=fig)
plt.suptitle('SDSS Quasars')
Out[10]:
In [ ]: