The purpose of this code is to compute the Absolute Magnitude of the candidates (Both KDE and RF) and plot them as a function of redshift. The conversion is as follows:

$M = m - 5log_{10}(\frac{d}{10\mathrm{pc}})$ or, in Mpc

$M = m - (5log_{10}(\frac{d}{1\mathrm{Mpc}})-5log_{10}(10^{-5})) == m - 5log_{10}(\frac{d}{1\mathrm{Mpc}})-25$

where d is the Luminosity distance in parsecs. The distance is estimated using the photoz's and a $\Lambda$CDM Cosmology with parameters:

$H_0 = 70 Mpc^{-1}, \Omega_{\Lambda} = 0.726, \Omega_{M} = 0.274$

Using the Friedmann Equation:

$ $

We can compute distance and, using the i_mag I can determine the absolute magnitude at z=0. I Finally need to convert to $M_i[z=2]$ by using the $\alpha_\nu$ values from Richards 2006/Ross2013


In [55]:
%matplotlib inline
import os
import sys
sys.path.insert(0, '/home/john/densityplot/densityplot')
from densityplot.hex_scatter import hex_contour as hex_contour
import numpy as np
from astropy.io import fits as pf
import camb
from camb import model
import matplotlib.pyplot as plt
from matplotlib import gridspec

In [37]:
#open the candidate data
#path = '/Users/johntimlin/Clustering/Combine_SpIES_Shela/Data_sets/Match_SpSh_Cand_wzcorrected_nooutlier_allinfo.fits'
#path = '/Users/johntimlin/Catalogs/QSO_candidates/201606/All_hzcandidate_correctphotoz_fromgaussian_allinfo.fits'
path = '../Data_Sets/QSO_Candidates_allcuts_with_errors_visualinsp.fits'
data = pf.open(path)[1].data
rz = (data.zphotNW>=2.9) & (data.zphotNW<=5.4) & (data.Good_obj == 0) & (data.dec>=-1.2) & (data.dec<=1.2)
rshift = data.zphotNW[rz]

print rshift.dtype
r = rshift.astype('float64') #camb's luminosity distance calculator only accepts float 64 data types
print r.dtype


>f8
float64

In [46]:
print len(r)


8838
[ 20.21052725  20.38572509  19.46829741 ...,  21.39967629  20.56156924
  20.80079899]

In [27]:
#Open the Shen2007 data
shendat = '/Users/johntimlin/Clustering/Shen_test/Data/Shen2007_Clustering_sample.fits'
sdat = pf.open(shendat)[1].data
#Cut to their objects and get array of redshifts
sdx = (sdat.Sfl == 1) #& (sdat.z>=3.5) & (sdat.z<=5.4)
srz = sdat.z[sdx]
sr = srz.astype('float64')
#print simag

In [28]:
#First define Planck 2015 cosmological parameters
H = 70 #H0. 
oc = 0.229 #physical density of CDM 
ob = 0.046 #physical density of baryons
#Set up parameters in CAMB
pars = camb.CAMBparams()
#Conversion to density param: Omega_Matter = (oc+ob)/(H0/100.)**2
#Hard code the cosmolgy params
pars.H0=H #hubble param (No h!!)
pars.omegab=ob #Baryon density parameter
pars.omegac=oc #CDM density parameter
pars.omegav=0.725 #Vacuum density parameter
pars.set_dark_energy()

#Set up parameters in CAMB
pars = camb.CAMBparams()
#H0 is hubble parameter at z=0, ombh2 is the baryon density (physical), omch2 is the matter density (physical)
#mnu is sum of neutrino masses, omk is curvature parameter (set to 0 for flat), meffsterile is effective mass of sterile neutrinos
pars.set_cosmology(H0=H,ombh2=ob, omch2=oc,omk=0)#,mnu=0,meffsterile=0) 
pars.set_dark_energy()

bkg = camb.get_background(pars) #Background parameters 


Ldist = bkg.luminosity_distance(r) #Luminosity distance for SpIES cand
ShLdist = bkg.luminosity_distance(sr) #Luminosity distance for Shen qso


#Make the i_mag line targeted for shen 2007
sampz = np.linspace(0.5,5.4,100)
line = bkg.luminosity_distance(sampz)
const_mag = np.ones(len(line))*20.2
const_magsp = np.ones(len(line))*23.3

In [29]:
#Compute the absolute magnitude at z=0
M = (22.5-2.5*np.log10(data.iflux[rz])) - 5.0*np.log10(Ldist) - 25.0
M202 = const_mag - 5.0*np.log10(line) - 25.0
M23 = const_magsp - 5.0*np.log10(line) - 25.0
shenM = sdat['imag'][sdx] - 5.0*np.log10(ShLdist) - 25.0

In [30]:
#Compute the corrections to apply to M[z=0] to get M[z=2]
def Kcorr(z,alpha = -0.5):
    #Ross13
    K13 = -2.5*np.log10(1+z) - 2.5*alpha*np.log10(1+z)+2.5*alpha*np.log10(1.0+2.0)
    return K13
#import the K-corrections from Richards 2006
K06 = pf.open('./K_correct_Richards06.fits')[1].data
K13 = pf.open('./K_correct_Ross13.fits')[1].data

#Pull out the redshift information from the data for SpIES and Shen
rshifts = data.zphotNW[rz]
srshift = sdat.z[sdx]

#Round the redshifts to 2 decimal places so that I can match to the correction values in Richards 2006
roundz = np.round(rshifts,decimals = 2)
roundt = np.round(sampz,decimals = 2)
roundsz = np.round(srshift,decimals = 2)

#Find the correction value that corresponds to the redshift in the file
Kcor=[]
Ktest = []
Ktestsp = []
Kshen = []
for i in roundz:
    kc = K06.KCorr[np.where(K06.z == i)]
    Kcor.append(kc[0])
    
for j in roundt:
    kt = K06.KCorr[np.where(K06.z == j)] 
    Ktest.append(kt[0])
    Ktestsp.append(kt[0])
for m in roundsz:
    kt = K06.KCorr[np.where(K06.z == j)] 
    Kshen.append(kt[0])


KC = np.asarray(Kcor)
KT = np.asarray(Ktest)
KS = np.asarray(Kshen)

#Correct the Absolute values using the K-corrections found above
dcorrect = M-KC
lcorrect = M202-Ktest
spcorrect= M23 - Ktestsp
scorrect = shenM - Kshen

In [31]:
#Plotting Parameters (Replace with Group code call!)
params = {'legend.fontsize': 16, 'xtick.labelsize': 20, 'ytick.labelsize': 20, 'xtick.major.width':2, 'xtick.minor.width':2, 'ytick.major.width':2, 'ytick.minor.width':2, 'xtick.major.size':8, 'xtick.minor.size':6, 'ytick.major.size':8, 'ytick.minor.size':6}
plt.rcParams.update(params)
plt.rc("axes", linewidth=3.0)

In [32]:
plt.figure(1,figsize = (8,8))
plt.scatter(rshift,dcorrect,color = '#fd8d3c',edgecolor = None,s=1,alpha = 0.9)#,label = 'Timlin 2016 QSO sample' )
plt.scatter(1,1,s=80,color = '#fd8d3c',label = 'This study')
plt.scatter(srz,scorrect,color='#e31a1c',edgecolor = None,s=1,alpha = 0.9)#,label = 'Shen 2007 QSO sample')
plt.scatter(1,1,s=80,color = '#e31a1c',label = 'Shen 2007')

plt.plot(sampz,lcorrect,color = 'k',linewidth = 2,label = r'$i$=20.2')
#plt.plot(sampz,spcorrect,color = 'k',linestyle = '--', dashes=(10,5,10,5),linewidth = 2,label = r'$i$=23.3')
plt.xlim(2.8,5.3)
plt.ylim(-30.5,-22.5)
plt.xlabel('Redshift',fontsize = 18)
plt.ylabel(r'$M_i$[z=2]',fontsize = 18)
plt.gca().invert_yaxis()
plt.minorticks_on()
plt.legend(loc=4,scatterpoints = 1)

#plt.savefig('Absolute_Mag_SpIES_Shen.pdf')


Out[32]:
<matplotlib.legend.Legend at 0x1139acb90>

In [33]:
testcol = np.linspace(-2,3.5,100)
ilim = 21.3-testcol
ch2lim = 22+testcol

plt.figure(2,figsize = (10,10))
ax1 = plt.subplot(211)
ax1.scatter(data.imag-data.ch2_mag,data.ch2_mag, edgecolor = None,s=1, alpha = 0.2, color = 'g')
ax1.set_ylabel('ch2_mag',size = 14)
ax2 = plt.subplot(212,sharex = ax1)
ax2.scatter(data.i_mag-data.ch2_mag,data.i_mag, edgecolor = None,s=1, alpha = 0.2,color = 'r')
ax2.set_ylabel('i_mag',size = 14)

ax1.plot(testcol,ilim,linestyle = '--',dashes = (8,4,8,4),color = 'k',label = 'i_mag = 21.3')
ax2.plot(testcol,ch2lim,linestyle = '--',dashes = (8,4,8,4),color = 'b',label = 'ch2_mag = 22')

ax1.axhline(22.0,linestyle = '--',dashes = (8,4,8,4),color = 'b',label = 'ch2_mag = 22')
ax2.axhline(21.3,linestyle = '--',dashes = (8,4,8,4),color = 'k',label = 'i_mag = 21.3')
plt.legend(loc = 2)
plt.xlabel('i-ch2',size = 14)
plt.xlim(-1.5,3)


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-33-c73e4f8c401e> in <module>()
      5 plt.figure(2,figsize = (10,10))
      6 ax1 = plt.subplot(211)
----> 7 ax1.scatter(data.imag-data.ch2_mag,data.ch2_mag, edgecolor = None,s=1, alpha = 0.2, color = 'g')
      8 ax1.set_ylabel('ch2_mag',size = 14)
      9 ax2 = plt.subplot(212,sharex = ax1)

/Users/johntimlin/anaconda/lib/python2.7/site-packages/numpy/core/records.pyc in __getattribute__(self, attr)
    444             res = fielddict[attr][:2]
    445         except (TypeError, KeyError):
--> 446             raise AttributeError("recarray has no attribute %s" % attr)
    447         obj = self.getfield(*res)
    448 

AttributeError: recarray has no attribute ch2_mag

In [ ]:


In [53]:
imag = 22.5-2.5*np.log10(data.iflux[rz])

num,bins = np.histogram(imag,bins='fd')

print bins

fig = plt.figure(5,figsize = (8,4))
plt.hist(imag,bins, histtype = 'step',normed = False,color = '#FFA500',linewidth = 2)
plt.hist(sdat['imag'][sdx],bins, histtype = 'step',normed = False,color = 'r',linewidth = 2)
plt.xlabel(r'$i$ Magnitude',fontsize = 14)
plt.ylabel('Number',fontsize = 14)
#plt.savefig('imag_hist.pdf',bbox_inches='tight')

plt.show()


[ 17.18215949  17.25545578  17.32875206  17.40204834  17.47534462
  17.54864091  17.62193719  17.69523347  17.76852976  17.84182604
  17.91512232  17.9884186   18.06171489  18.13501117  18.20830745
  18.28160374  18.35490002  18.4281963   18.50149258  18.57478887
  18.64808515  18.72138143  18.79467771  18.867974    18.94127028
  19.01456656  19.08786285  19.16115913  19.23445541  19.30775169
  19.38104798  19.45434426  19.52764054  19.60093683  19.67423311
  19.74752939  19.82082567  19.89412196  19.96741824  20.04071452
  20.1140108   20.18730709  20.26060337  20.33389965  20.40719594
  20.48049222  20.5537885   20.62708478  20.70038107  20.77367735
  20.84697363  20.92026992  20.9935662   21.06686248  21.14015876
  21.21345505  21.28675133  21.36004761  21.43334389  21.50664018
  21.57993646  21.65323274  21.72652903  21.79982531  21.87312159
  21.94641787  22.01971416  22.09301044  22.16630672  22.23960301
  22.31289929  22.38619557  22.45949185  22.53278814  22.60608442
  22.6793807   22.75267698  22.82597327  22.89926955  22.97256583
  23.04586212  23.1191584   23.19245468  23.26575096  23.33904725
  23.41234353  23.48563981  23.5589361   23.63223238  23.70552866
  23.77882494  23.85212123  23.92541751]

In [75]:
imag = 22.5-2.5*np.log10(data.iflux[rz])

num,bins = np.histogram(imag,bins='fd')

fig = plt.figure(5,figsize = (6,12))
gs = gridspec.GridSpec(2, 1, height_ratios=[0.6,0.4])

ax0 = plt.subplot(gs[0])
ax1 = plt.subplot(gs[1],)

plt.axes(ax0)

plt.scatter(rshift,dcorrect,color = '#fd8d3c',edgecolor = None,s=1,alpha = 0.9)#,label = 'Timlin 2016 QSO sample' )
plt.scatter(1,1,s=80,color = '#fd8d3c',label = 'This study')
plt.scatter(srz,scorrect,color='#e31a1c',edgecolor = None,s=1,alpha = 0.9)#,label = 'Shen 2007 QSO sample')
plt.scatter(1,1,s=80,color = '#e31a1c',label = 'Shen 2007')
plt.plot(sampz,lcorrect,color = 'k',linewidth = 2,label = r'$i$=20.2')
#plt.plot(sampz,spcorrect,color = 'k',linestyle = '--', dashes=(10,5,10,5),linewidth = 2,label = r'$i$=23.3')
plt.xlim(2.8,5.3)
plt.ylim(-30.5,-22.5)
plt.xlabel('Redshift',fontsize = 16)
plt.ylabel(r'$M_i$[$z$=2]',fontsize = 16)
plt.gca().invert_yaxis()
plt.minorticks_on()
leg =plt.legend(loc=4,scatterpoints = 1)
leg.get_frame().set_alpha(0.35)

plt.axes(ax1)
plt.hist(sdat['imag'][sdx],bins, histtype = 'step',normed = False,color = 'r',linewidth = 2,label= 'Shen 2007')
plt.hist(imag,bins, histtype = 'step',normed = False,color = '#FFA500',linewidth = 2,label = 'This study')

plt.xlabel(r'$i$-Magnitude',fontsize = 16)
plt.ylabel('Number',fontsize = 16)
plt.minorticks_on()
leg =plt.legend(loc=2)
leg.get_frame().set_alpha(0.35)

plt.savefig('Absolute_Mag_SpIES_Shen.pdf',bbox_inches='tight',pad_inches=0.5)



In [12]:
#Open the full test set
testset = '/Users/johntimlin/Catalogs/QSO_candidates/201606/GTR-ADM-QSO-ir_good_test_2016n.fits'

testobj = pf.open(testset)[1].data

plt.figure(3,figsize = (10,10))
ax1 = plt.subplot(211)
ax1.scatter(testobj.i-testobj.s2,testobj.s2, edgecolor = None,s=1, alpha = 0.2, color = 'g')
#hex_contour(testobj.i-testobj.s2,testobj.s2,std = True,levels =[0.3,0.5,0.7,0.9],smoothing = 3,hkwargs={'gridsize':20},ckwargs = {'colors':'g'})
ax1.set_ylabel('ch2_mag',size = 14)
ax2 = plt.subplot(212,sharex = ax1)
ax2.scatter(testobj.i-testobj.s2,testobj.i, edgecolor = None,s=1, alpha = 0.2,color = 'r')
#hex_contour(testobj.i-testobj.s2,testobj.i,std = True,levels =[0.3,0.5,0.7,0.9],smoothing = 3,hkwargs={'gridsize':20},ckwargs = {'colors':'r'})
ax2.set_ylabel('i_mag',size = 14)

ax1.plot(testcol,ilim,linestyle = '--',dashes = (8,4,8,4),color = 'k',label = 'i_mag = 21.3')
ax2.plot(testcol,ch2lim,linestyle = '--',dashes = (8,4,8,4),color = 'b',label = 'ch2_mag = 22')

ax1.axhline(22.0,linestyle = '--',dashes = (8,4,8,4),color = 'b',label = 'ch2_mag = 22')
ax2.axhline(21.3,linestyle = '--',dashes = (8,4,8,4),color = 'k',label = 'i_mag = 21.3')
plt.legend(loc = 2)
plt.xlabel('i-ch2',size = 14)
plt.xlim(-5,5)


Out[12]:
(-5, 5)

In [ ]:


In [ ]: