This notebook uses Python2.7. You may need to change some syntax if you use Python3.X
The purpose of this notebook was to perform some testing:
We do not investigate photometric redshift accuracies between different photometric surveys
Apart from basic python modules (eg. numpy, matplotlib etc.), the code requires ThreedHST Module available at https://github.com/gbrammer/threedhst.
The other packages can be easily installed using a package manager. eg. pip
In [1]:
%matplotlib inline
from zfire_utils import *
In [2]:
#load data
UDS, UKIDSS_selected = prepare_uds()
COSMOS, MOSDEF_ZFOURGE, VUDS_ZFOURGE, VUDS_extra = prepare_cosmos()
Let's investigate the data
In [3]:
print "Total number of targeted galaxies in UDS", len(UDS)
print "Q=1 N --> ", len(UDS.query('conf==1'))
print "Q=2 N --> ", len(UDS.query('conf==2'))
print "Q=3 N --> ", len(UDS.query('conf==3'))
In [4]:
print "Total number of targeted galaxies in COSMOS", len(COSMOS)
print "Q=1 N --> ", len(COSMOS.query('conf==1'))
print "Q=2 N --> ", len(COSMOS.query('conf==2'))
print "Q=3 N --> ", len(COSMOS.query('conf==3'))
In [5]:
print "Targeted galaxies K band N = ", len(COSMOS[~pd.isnull(COSMOS.zspec_K)])
print "Galaxies outside H-alpha range in K N = ", len(COSMOS.query('conf>1 and (zspec<1.90 | zspec>2.66)'))
print "H-alpha detected galaxies Q=3 N = ", len(COSMOS.query('conf==3 and fha_snr>=5 and zspec>1.90 and zspec<2.66'))
print "H-alpha detected galaxies Q=2 N = ", len(COSMOS.query('conf==2 and fha_snr>=5 and zspec>1.90 and zspec<2.66'))
print "H-alpha non-detected galaxies N = ", len(COSMOS.query('fha_snr<5'))
In [ ]:
Henceforth the plots of the ZFIRE survey paper are made
In [6]:
#give path to eazy here
eazydir = '/Users/temp/Documents/ZFOURGE/Data/release_v2.2/cosmos/EAZY/cosmos_eazy_all/'
In [7]:
def stack_pz():
for index, values in ZFIRE_observed_K.T.iteritems():
object_ID = int(values.Nameobj)
zgrid, pzi = ep.getEazyPz(object_ID-1, MAIN_OUTPUT_FILE='cosmos.v0.10.7.a', OUTPUT_DIRECTORY=eazydir, CACHE_FILE='Same')
ll, ul, sky_loss = get_limits(values)
within_K = (zgrid>ll) & (zgrid<ul)
in_K = sp.integrate.trapz(pzi[within_K],x=zgrid[within_K]) * (1-sky_loss)
total = sp.integrate.trapz(pzi,x=zgrid)
P = in_K/total
if len(all_zgrid)==0:
all_zgrid.append(zgrid)
all_pzi.append(pzi)
inside_all.append(in_K)
total_all.append(total)
all_pzi_sum = np.sum(all_pzi, axis=0)
return all_pzi, inside_all, total_all , all_zgrid, all_pzi_sum
def open_fits_COSMOS(galaxy, band ):
"""
v2.1.1: commented the emission line fits opening since it is not needed for the
purpose of p(z) stacks
version 2.1:
added to open the emission line fit files
upgrading to open multiple object spectra correctly.
version 1.3 21/11/14
upgrade to take spectra from the master table. The file open order should be:
1. the common folder to priorotize the objects observed in multiple observing runs
2. objects in individual observing runs
For H band option 1 doesn't apply
Multiobject option has been removed.
"""
if pd.isnull(galaxy.doubles) is True:
suffix_string = ''
ID = galaxy.Nameobj
else:
double_string = galaxy.doubles
if s.find(double_string, 'b')!=-1:
suffix_string = '-2'
ID = s.rstrip(galaxy.doubles, 'b')
elif s.find(double_string, 'c')!=-1:
suffix_string = '-3'
ID = s.rstrip(galaxy.doubles, 'c')
if band =='H':
path='../../mosdrp/analysis/spectra/spectra_1d/2014feb_1d/after_scaling/spectra/'
try:
fits = glob.glob(str(path)+'Hbandmask*'+ str(ID) +'_*_1D'+suffix_string+'.fits');# print fits; print path
Name = fits[0]; eps = pf.open(Name)
except IndexError:
print str(Object['Nameobj'])+" Object not found in H band"
return -99, -99, -99
print 'Opened ' + str(Name)
if s.find(Name, 'Hbandmask1')!=-1:
mask = 'Hbandmask1'; ET = Hbandmask1; print 'This is----> ' + str(mask)
elif s.find(Name, 'Hbandmask2')!=-1:
mask = 'Hbandmask2'; ET = Hbandmask2; print 'This is----> ' + str(mask)
else:
mask= 'unknown' ; print '**ERROR** mask not recognized: Please Check'
obs_run= 'feb'
elif band=='K':
try:
path='../../mosdrp/analysis/spectra/spectra_1d/common_1d/DR1/after_scaling_common_1D/'
fits = glob.glob(str(path)+'*_'+str(ID)+'_coadd1D.fits')
Name = fits[0]; eps = pf.open(Name); print "opened "+ str(Name)
assert len(fits)==1, "There are 0/multiple matches"
except IndexError:
path='../../mosdrp/analysis/spectra/spectra_1d/201*_1d/after_scaling/spectra/'
fits = glob.glob(str(path)+'*_K_K_*_'+ str(ID) +'_*_1D'+suffix_string+'.fits')
Name = fits[0]; eps = pf.open(Name); print "opened "+ str(Name)
assert len(fits)==1, "There are 0/multiple matches"
return eps
def set_sky_weights(band, wave, w):
"""version 1.0
sets weights around sky lines to be 0 and everything else to be 1
the range is determined according to the spectral resolution"""
if band=='H': sky_lines= sky_H['wavelength']; spec_res = 4.5
elif band=='K': sky_lines= sky_K['wavelength']; spec_res = 5.5
elif band=='J': sky_lines= sky_J['wavelength']; spec_res = 4.0
elif band=='Y': sky_lines= sky_Y['wavelength']; spec_res = 3.5 #not checked
else: print 'Unknown Band'
for i, v in enumerate(sky_lines):
wave_mask = np.ma.masked_outside(wave, sky_lines[i]-spec_res, sky_lines[i]+spec_res)
masked_array = np.ma.getmaskarray(wave_mask)
np.place(w, masked_array==False, [0])
return w
def get_limits(galaxy):
data = open_fits_COSMOS(galaxy,'K')
scidata, sddata, wavelength, hdr = data[0].data, data[1].data, data[2].data, data[0].header
CRVAL1, CD1_1 , CRPIX1 = hdr['CRVAL1'], hdr['CD1_1'], hdr['CRPIX1']
i_w = np.arange(len(scidata))+1 #i_w should start from 1
wavelength = ((i_w - CRPIX1) * CD1_1 ) + CRVAL1
limits = np.nonzero(sddata)#masking procedure is fine. Checked 1/09/14
photometry_mask = np.ma.masked_inside(wavelength, wavelength[limits[0][0]], wavelength[limits[0][-1]] , copy=True)
photometry_mask = np.ma.getmaskarray(photometry_mask)
Lambda_limits = wavelength[photometry_mask]
flux_limits = scidata[photometry_mask]
error_limits = sddata[photometry_mask]
sky_weights = np.ones_like(Lambda_limits)
sky_weights = set_sky_weights('K',Lambda_limits, sky_weights)
print sky_weights[0:20]
fraction_lost = float(len(Lambda_limits[sky_weights==0]))/len(Lambda_limits)
print "fraction lost due to sky = ",fraction_lost
limit_low = (Lambda_limits[0]/6565)-1
limit_upper = (Lambda_limits[-1]/6565)-1
#make_spectra_plot(Lambda_limits,flux_limits,error_limits,galaxy.zspec, galaxy.Nameobj)
return limit_low, limit_upper, fraction_lost
def plot_pz_from_binary(ID,x,y,data,ll, ul, prob,spec_z_COSMOS=0):
fig = plt.figure(figsize=(4,4))
ax1 = fig.add_subplot(111)
xlow=-0.1; xup=4.9
ax1.plot(x[(x<3) & (x>1.5)], (y/max(y))[(x<3) & (x>1.5)], color='black', lw=2, ls='-')
ax1.axvline(ul, c='k', ls=':')
ax1.axvline(ll, c='k', ls=':')
ax1.set_xlabel('$z\mathrm{_{photo}}$', fontsize=18)
ax1.set_ylabel('$\mathrm{Likelihood}$', fontsize=18)
ax1.tick_params(axis='both', which='major', labelsize=15)
ax1.set_xlim(1.4,3.1)
ax1.set_xticks([1.5,2.0,2.5,3.0])
ax1.set_ylim(0,max(y/max(y))+0.1)
props = dict(boxstyle='round', facecolor='c', alpha=0.5)
if ID == 'ALL':
"""the histogram is assigned to the main panel with shared x and y axis"""
print 'This is for the total population'
ax2 = ax1.twinx()
multiples_removed = spec_z_COSMOS[~((spec_z_COSMOS.index.str.contains("b"))| (spec_z_COSMOS.index.str.contains("c")) |
(spec_z_COSMOS.index.str.contains("A")))]
Qz3 = (multiples_removed['zspec'])
Qz3 = np.asarray(Qz3)
x,bins,p = plt.hist( Qz3 , bins=5, color='limegreen',alpha=0.30,
label = 'zspec', normed=False)
print x, bins, p
for item in p: #normalize to a max of 1
item.set_height(item.get_height()/np.max(x))
ax2.set_ylim(0,1.1)
ax2.set_xlim(1.4,3.1)
ax2.set_xticks([1.5,2.0,2.5,3.0])
ax2.set_ylabel('$\mathrm{Normalized\ number\ of\ galaxies}$', color='green', fontsize=18)
ax2.tick_params(axis='both', which='major', labelsize=15)
ax2 = ax1.twiny()
ax2.set_xlabel('$z\mathrm{_{spec}}$', color='green', fontsize=18)
#ax2.set_xlim(xlow, xup)
ax2.set_xlim(1.4,3.1)
ax2.set_xticks([1.5,2.0,2.5,3.0])
ax2.tick_params(axis='both', which='major', labelsize=15)
elif ID == 'conf = 1':
print 'This is for non-detections'
elif values['conf']>1:
plt.axvline(data['zspec'], c='grey', ls='-.')
textstr= '$\mathrm{zspec}=%.2f$\n$\mathrm{zpeak}=%.2f$\n$\mathrm{P_{inK}}=%.2f$'%(data['zspec'],data['z_peak'], prob)
ax.text(0.65, 0.95, textstr, transform=ax.transAxes, fontsize=14,verticalalignment='top', bbox=props)
else:
textstr= '$\mathrm{zpeak}=%.2f$\n$\mathrm{P_{inK}}=%.2f$'%(data['z_peak'], prob)
ax.text(0.65, 0.95, textstr, transform=ax.transAxes, fontsize=14,verticalalignment='top', bbox=props)
plt.tight_layout()
plt.savefig('outputs/'+str(ID)+'.pdf', dpi=400)
plt.close()
In [8]:
sky_K = ascii.read('data/sky_lines/K.dat')
ZFIRE_observed_K = COSMOS[~pd.isnull(COSMOS.zspec_K)]
#remove 3602 since it has been observed twice
ZFIRE_observed_K= ZFIRE_observed_K[ZFIRE_observed_K.index != '3602']
In [48]:
In [49]:
from astropy.io import fits as pf
import threedhst.eazyPy as ep
inside_all=[]
total_all=[]
all_zgrid=[]; all_pzi=[]
ZFIRE_sample = ZFIRE_observed_K[(ZFIRE_observed_K.conf==1) | (((ZFIRE_observed_K.zspec>1.90) & (ZFIRE_observed_K.zspec<2.66))
& ((ZFIRE_observed_K.conf==2) | (ZFIRE_observed_K.conf==3)))]
ZFIRE_observed_K = ZFIRE_observed_K[~pd.isnull(ZFIRE_observed_K.Ks)]
all_pzi, inside_all, total_all , all_zgrid, all_pzi_sum = stack_pz()
P = np.sum(inside_all)/np.sum(total_all)
print "Total area, area inside K, fraction ",np.sum(total_all),np.sum(inside_all), P
spec_z_COSMOS = ZFIRE_sample[ZFIRE_sample['conf']>1]
spec_z_COSMOS.loc[:,'zphoto'] = spec_z_COSMOS['z_peak']
spec_z_COSMOS = spec_z_COSMOS[['zphoto','zspec', 'conf']]
plot_pz_from_binary('ALL', all_zgrid[0],all_pzi_sum/np.sum(total_all),"leave empty",1.90, 2.66,P , spec_z_COSMOS)
In [ ]:
In [ ]: