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
In [46]:
print len(r)
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]:
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)
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()
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]:
In [ ]:
In [ ]: