I hate IDL

I've been working a little bit in IDL for this new project because all the data is in fits files, Alexie has a lot of code that is in IDL which we are using to handle the CFHT data.

Anyway, I spent all morning struggling trying to do very simple input and output from this stupid stupid language.

I learned a few things, and since I hardly ever code in IDL, I am sure I'll forget unless I write them down right now. So here we go:

1) Tag names

You can get the tag names from a structure like this:

tname = tag_names(struct)

2) Merging structures

You can merge two structures as follows:

str3 = create_struct(str1,str2,NAME='str3')

3) writecol bug

For some reason when I try to use writecol to write out data, if there are more that 6 columns, and I don't specify the format, it puts a return carriage in the middle of the data. This goes away if you specify the format:

thisfile = './qsosMergedInfo.dat'
writecol,thisfile, qsos.ra, qsos.dec, qsos.z, qsos.zwarning, $
qsos.flux_clip_mean[0],qsos.flux_clip_mean[1], qsos.flux_clip_mean[2], $
qsos.flux_clip_mean[3], qsos.flux_clip_mean[4], fmt='(f,f,f,f,f,f,f,f,f)'

Testing QSO-Galaxy Cross CF

I spent the last couple days making sure that I trusted the correlation function code I had developed for this project. I compared it to the results from a (slow) sphere-match code on a small data set.

The following are the testing codes (also here ../logs/cross_corr_check.idl). In IDL using spherematch:

; Read in data/random files

readcol,'cfht_data_tiny.dat',rad1,decd1, format = ('D,D')

readcol, 'qso_data_tiny.dat',rad2,decd2,cd, format = ('D,D,D') ;cd is comoving distance of qsos

readcol,'cfht_randoms_tiny.dat',rar1,decr1, format = ('D,D')

; Find number of objects in files

nr1 = n_elements(rar1)

nd1 = n_elements(rad1)

nd2 = n_elements(rad2)

; Correlate out to a large theta to get all pairs

thetamax = 50.0

print,'Starting cross-correlation:'

print,'Estimating DD...'

spherematch,rad1,decd1,rad2,decd2,thetamax,ind1a,ind1b,dist_dd,maxmatch=0

; Convert from angular separation to comoving distance separation

this_dd = 2*sin(dist_dd*!pi/360)*cd[ind1b]

;Bins go from 0.1 to 10 with 15 bins.

corrmin = 0.1D

corrmax = 10.0D

nbins = 15.0D

; Find bins lower, upper and centers

bins_lower = (corrmax-corrmin)/(nbins)*findgen(nbins)+corrmin

bins_upper = (corrmax-corrmin)/(nbins)*(findgen(nbins)+1)+corrmin

rmean = fltarr(nbins)

for i = 0,(nbins-1) do rmean[i] = (bins_lower[i]+bins_upper[i])/2.

; Bin the DD separation distances

dd = fltarr(nbins)

for i = 0,(nbins-1) do dd[i] = n_elements(where(this_dd gt bins_lower[i] AND this_dd le bins_upper[i]))

print,'Estimating DR...'

spherematch,rar1,decr1,rad2,decd2,thetamax,ind1,ind2,dist_dr1,maxmatch=0

this_dr = 2*sin(dist_dr1*!pi/360)*cd[ind2]

dr = fltarr(nbins)

for i = 0,(nbins-1) do dr[i] = n_elements(where(this_dr ge bins_lower[i] AND this_dr le bins_upper[i]))

corr1 = 1L*dd/dr*1L*(nd2*nr1)/(1L*nd1*nd2)-1L

for i = 0,(nbins-1) do print, rmean[i], corr1[i]

Separation omega

0.430000 -0.115686

1.09000 -0.104478

1.75000 -0.120804

2.41000 -0.0914845

3.07000 -0.0393971

3.73000 -0.0268416

4.39000 0.0134841

5.05000 0.0596094

5.71000 0.0227162

6.37000 0.102554

7.03000 0.0929233

7.69000 0.0900670

8.35000 0.0591398

9.01000 0.0284724

9.67000 0.0598689

(Note that these "tiny" files only have 6 qsos and 9000 galaxies, so the correlation function values are very noisy, this was just to test and I used small files to)

By comparison I also have the python/C code which runs much faster (../Jessica/qsobias/Correlate/runCorrelation.py):

import numpy as N

from pylab import *

from correlationFunctions import *

#------------------------------------------------------------------------

# Create file names (tiny catalogs)

#------------------------------------------------------------------------

workingDir = 'tinyrun'

makeworkingdir(workingDir)

galaxyDataFile, qsoDataFile, randomDataFile, corr2dCodefile, argumentFile, runConstantsFile = makeFileNamesTiny(workingDir)

oversample = 5. # Amount that randoms should be oversampled

corrBins = 25.0 # Number of correlation bins (+1)

mincorr = 0.1 # (Mpc/h comoving distance separation) Must be great than zero if log-binning

maxcorr = 10.0 # (Mphc/h comoving distance separation)

convo = 180./pi # conversion from degrees to radians

tlogbin = 1 # = 0 for uniform spacing, = 1 for log spacing in theta

#------------------------------------------------------------------------

# Write run constants to a file

#------------------------------------------------------------------------

writeRunConstantsToFile(runConstantsFile, galaxyDataFile, qsoDataFile, \

randomDataFile, corr2dCodefile, argumentFile, oversample, corrBins, \

mincorr, maxcorr, tlogbin)

#------------------------------------------------------------------------

# Compute the Angular Correlation Function

#------------------------------------------------------------------------

runcrossCorrelation(workingDir, argumentFile, corr2dCodefile, galaxyDataFile,\

qsoDataFile, randomDataFile, mincorr, maxcorr, corrBins, tlogbin)

# separation (Mpc/h) crossw (Mpc/h)

0.4300000000 -0.1156862745

1.0900000000 -0.1044776119

1.7500000000 -0.1208039566

2.4100000000 -0.0914845135

3.0700000000 -0.0393970538

3.7300000000 -0.0268417043

4.3900000000 0.0134841235

5.0500000000 0.0596093513

5.7100000000 0.0227161938

6.3700000000 0.1025539385

7.0300000000 0.0929232804

7.6900000000 0.0900670231

8.3500000000 0.0591397849

9.0100000000 0.0284723490

9.6700000000 0.0598689436

As you can see the correlation functions match!

Final Numbers for Likelihood Paper Table 1

Ok, I finally finished doing the calculations for efficiency/completeness table. I've tried to comment the below code pretty well, so hopefully it makes sense to me when I try to reproduce this in the future (also in the following log file: ../logs/110215log.pro):

; Read in the targeting file, made by Adam Myers and is 219.9 deg^2 in area
adamfile = '/home/jessica/boss/chunk1_wcomp.fits'
targets = mrdfits(adamfile, 1)
stripe82targets1 = targets
newratio = stripe82targets1.LIKE_RATIO_CORE

;Read in coadded targeting file, given to me by Adam, and has the
;coadded photometry and the coadded likelihood ratios (like_ratio_bonus)
coaddfile = '/home/jessica/boss/star82-varcat-bound-ts-wgalex-wnnvarbonus-galexbitset.fits.gz'
coadd = mrdfits(coaddfile, 1)

;Match the coadded targeting file to the 219.9 file above so that they
;are the same area
spherematch, stripe82targets1.ra, stripe82targets1.dec, coadd.ra, coadd.dec, 2./3600, i1, i2, d12
coadd = coadd[i2]

; Read in the truthtable file
truthfile = '/home/jessica/boss/BOSS_Quasars_3PCplus_v1.1_wcomp.fits'
collateinfo1 = mrdfits(truthfile, 1)
stripe82data1 = collateinfo1
;Constrain to High human confidence quasars (with z > 2.15)
confcut = where(stripe82data1.z_conf_person GE 2 and stripe82data1.Z_PERSON GE 2.15)
stripe82quasars1 = stripe82data1[confcut] ; with redshift > 2.15
lowqsocut = where(stripe82data1.z_conf_person GE 2 and stripe82data1.Z_PERSON GE 0.5 and stripe82data1.Z_PERSON LT 2.15)
stripe82lowquasars1 = stripe82data1[lowqsocut] ; with 0.5 < redshift < 2.15

;Cut to objects which are in areas of high spectroscopic completeness
compthresh = 0.5 ;adjust to make different rows of Adam's table
wphot = where(stripe82targets1.completeness ge compthresh)
wspec = where(stripe82quasars1.completeness ge compthresh)
wdata = where(stripe82data1.completeness ge compthresh)
wlow = where(stripe82lowquasars1.completeness ge compthresh)

;Find Area of the area with above spectroscopic completeness
areas = mrdfits('/home/jessica/boss/area_wcomp.fits',1)
areacol = where(areas.comp_thresh ge compthresh)
thisArea = areas[areacol[0]].area

;Limit to objects with above spectroscopic completeness
stripe82targets = stripe82targets1[wphot]
stripe82quasars = stripe82quasars1[wspec]
stripe82data = stripe82data1[wdata]
stripe82lowqso = stripe82lowquasars1[wlow]

;Find likelihood threshold at (targetdensity) tpsd
area = 219.9 ; we know the exact area of the chunk1 file
targetdensity = 40.0 ;change this number to get the rows in table 1
tpsd = round(targetdensity*area)
;Use the full targeting file (no cuts)
sortRatio = targets[reverse(sort(targets.LIKE_RATIO_CORE))].LIKE_RATIO_CORE
;Threshold is the threshold such that we target at the above targetdensity
threshold = sortRatio[tpsd-1]

;Similarly find the threshold using the coadded photometry (will be
;different then the single-epoch photometry threshold
sortCoadd = coadd[reverse(sort(coadd.LIKE_RATIO_BONUS))].LIKE_RATIO_BONUS
thresholdCoadd = sortCoadd[tpsd-1]

;Spectroscopic Completeness threshold
print, compthresh

;Area
print, thisArea

;Likelihood Threshold
print, threshold

;Likelihood Targets per square degree
thistpsd = n_elements(where(stripe82targets.like_ratio_core gt threshold))/thisArea
print, thistpsd

spherematch, stripe82targets.ra, stripe82targets.dec, stripe82quasars.ra, stripe82quasars.dec, 2./3600, i1, i2, d12
; Likelihood Quasars targeted per square degree (not corrected)
qpsd = n_elements(where(stripe82targets[i1].like_ratio_core gt threshold))/thisArea
print, qpsd

spherematch, stripe82targets.ra, stripe82targets.dec, stripe82lowqso.ra, stripe82lowqso.dec, 2./3600, i1, i2, d12
; Likelihood Low Redshift Quasars targeted per square degree
lqpsd = n_elements(where(stripe82targets[i1].like_ratio_core gt threshold))/thisArea
print, lqpsd

spherematch, stripe82targets.ra, stripe82targets.dec, stripe82data.ra, stripe82data.dec, 2./3600, i1, i2, d12
; Others targeted per square degree
otpsd = n_elements(where(stripe82targets[i1].like_ratio_core gt threshold))/thisArea
print, otpsd - (qpsd + lqpsd)

;Likelihood fibers not targeted
print, thistpsd-otpsd

;Total z > 2.15 qsos
totalqsos = n_elements(stripe82quasars)/thisArea
print, totalqsos

;Total 0.5 < z < 2.15 qsos
totallowqsos = n_elements(stripe82lowqso)/thisArea
print, totallowqsos

;Efficiency (not corrected)
print, qpsd/thistpsd

;Completeness (not corrected)
print, qpsd/totalqsos

;Dealing with the non-targeted quasars that should have been targeted
;in the coadded photometry, correcting the efficiency and completeness

;Get likelihood targets (pass like threshold)
liketargets = stripe82targets[where(stripe82targets.like_ratio_core gt threshold)]
;Find targets that were targeted
spherematch, liketargets.ra, liketargets.dec, stripe82data.ra, stripe82data.dec, 2./3600, i1, i2, d12

;All likelihood targets (above threhsold)
all = where(liketargets.ra GT -99999)
;Targets that were targeted
lobs = liketargets[i1]
;Targets that were missed
lnobs = setdifference(all, i1)

;Of the ones that were missed, get the coadded photometry
likemissed = liketargets[lnobs]
spherematch, likemissed.ra, likemissed.dec, coadd.ra, coadd.dec, 2./3600, i1, i2, d12
likemissedcoadd = coadd[i2]

;Number that were missed that should have been targeted, because their
;coadded like_ratio_bonus passes the coadded threshold
nmissedcoadd = n_elements(where(likemissedcoadd.like_ratio_bonus gt thresholdCoadd))

; Number per square degree that would have been quasars if they were
; targeted. Multiply the number missed by the (not corrected) efficiency
; We will add these into our numbers to make them correct
nmpsd = round(nmissedcoadd*qpsd/thistpsd)/thisArea

;Line from table 1 for the table
; targets/deg^2, threshold, total like targets, total QSOs found, QSOs missed, completeness, efficiency
print, targetdensity, threshold, thistpsd*thisarea, (qpsd+nmpsd)*thisarea, (totalqsos-qpsd-nmpsd)*thisarea, (qpsd+nmpsd)/(totalqsos+nmpsd), (qpsd+nmpsd)/thistpsd

And here is the pretty LaTeX table:

Using the Whole Footprint with Masks

For the past few days I've been working on using masks to run correlation functions on the whole BOSS footprint. It's required learning how to use polygon masks and also to create randoms in these masks. Here is the code/testing:

To mask BOSS data (spectroscopic set)

;Read in BOSS data
bossgalaxies = '/home/jessica/boss/spAll-v5_4_14.fits'
spall = mrdfits(bossgalaxies, 1)

;Trim BOSS data to be a galaxy (remove QSOs, starts, etc)
isgal = where(spall.specprimary EQ 1 $
AND (spall.boss_target1 AND 2L^0+2L^1+2L^2+2L^3+2L^7) NE 0 $
AND strmatch(spall.class,'GALAXY*') $
AND spall.zwarning EQ 0)

;Select galaxies
galaxies = spall[isgal]
ra = galaxies.ra
dec = galaxies.dec
z = galaxies.z

;Select objects in the masks
weights = maskdataboth(ra,dec,z)

; Readout ra, dec, redshift of spectroscopic objects inside mask
thisfile = '/home/jessica/boss/spectroRaDecZMasked.dat'
writecol,thisfile,ra,dec,z

The code for maskdataboth is here:
/home/jessica/repository/ccpzcalib/Jessica/pythonJess/maskdataboth.pro

'maskdataboth' filters the data through both the BOSS mask and Shirley's mask to make sure the data sets are in the same ra/dec footprint.


To Mask Shirley's Data (photometric set)

;Read in Shirley's Photometric Data
file1 = "/home/shirleyho/research/SDSS_PK/powerspec_dir/DATA/ra_dec_z_run_rerun_camcol.z_gt_0.01"
readcol, file1, pra, pdec, pz, RUN, rerun, camcol, x, x, x, x, x, x, x, x, format='(F,F,F,D,D,D,D,D,F,F,F,F,F,F)'

;Select objects in the masks
pweights = maskdataboth(pra,pdec,pz)

; Readout ra, dec, redshift of spectroscopic objects inside mask
thisfile = '/home/jessica/boss/photoRaDecZMasked.dat'
writecol,thisfile,pra,pdec,pz

Here are plots of the data sets to show that they do indeed fall in the same region of the sky:




Both these data sets are saved in the following directory: /home/jessica/boss/
The code to run this is here: ../logs/101210log.pro

maskdataboth.pro

Miscellaneous IDL Commands

qselect -u jessica -s Q | xargs qdel

window,xsize=700,ysize=600

colorobject = qsos.z_tot
xobject = ugcoaddcolor
yobject = grcoaddcolor

thiscolorbin = 0.2
thiscolorstart = 2.0

xtitle = 'u-g magnitude'
ytitle = 'g-r magnitude'
mtitle = 'Color-Color Diagram All Input QSOs'

tempplot, colorobject, xobject, yobject, thiscolorbin, thiscolorstart, xtitle, ytitle, mtitle

data = qsos.z_tot
datamin = 0.0
datamax = 5.0
binsize = 0.1

xtitle = 'Redshift Bin'
ytitle = 'QSOs in Bin'
mtitle = 'Redshift Distribution'

.com histplot.pro
histplot, data, datamin, datamax, binsize, xtitle, ytitle, mtitle

Useful IDL Commands

I need to do some work for the likelihood project and thus code in IDL again. Because I am not working in IDL very frequently, I tend to forgot commands and syntax. I thought it would be useful to compile a list of commands here so that I don't have to keep looking these up again and again.

number of elements in an array:
print, N_ELEMENTS(array)

unix command in idl: dollar sign ($) before command
$pwd
$ls

list all variables currently being used:
help

doc_library, 'spherematch' -- man page for the function

which, 'spherematch' --- tells you where in the IDL path the function is.

when recompiling idlutils need to use evilmake

!pi is constant pi

plothist (makes histogram plot)

when looking at a structure you use the /str command:
help, mystructure, /str

reading a fits file using mrdfits:
filedata = mrdfits(file2read, 1)

Combine two structures of the same size:
struct3 =struct_combine(struct1,struct2)
Concatenate two structure of the same type:
struct3 = [struct1, struct2]

Plotting Histogram (from here):

PRO t_histogram
data = [[-5,  4,   2,  -8,  1], $
[ 3,  0,   5,  -5,  1], $
[ 6, -7,   4,  -4, -8], $
[-1, -5, -14,   2,  1]]
hist = HISTOGRAM(data)
bins = FINDGEN(N_ELEMENTS(hist)) + MIN(data)
PRINT, MIN(hist)
PRINT, bins
SPLOT, bins, hist, YRANGE = [MIN(hist)-1, MAX(hist)+1], PSYM = 10, XTITLE = 'Bin Number', YTITLE = 'Density per Bin'
END 

If you add maxmatch=1 to your spherematch it will avoid duplicates (should default to this)


total with higher precision: total(...., /double)


in for loops:
for j = 0L, rsize-1 do begin

the L makes j a long int.


;Show color table
colorName = PickColorName(startColorName)

; Make heat map (warning, not very sophisticated)


colorobject = qsotemplate1.z_sim
xobject = ugqcolor1
yobject = grqcolor1

colorbin = 0.3
colorstart = 1.0
colorstop = colorstart + colorbin

plot, xobject[0:1], yobject[0:1], ps=null, color=fsc_color('white'), xr=[-0.5,4],yr=[-0.5,2], XTITLE = 'u-g magnitude', YTITLE = 'g-r magnitude', TITLE = 'Color-Color Diagram QSO Catalog Jiang Combo', charsize = 1.5, charthick = 1

zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('dark red') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('red') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('orange red') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('orange') ;Fake QSO
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('gold') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('lawn green') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('lime green') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('dark green') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('teal') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('dodger blue') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('royal blue') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('blue') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('navy') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('dark slate blue') ;Fake QSOs
colorstart = colorstop
colorstop = colorstart + colorbin
zrange = where(colorobject GE colorstart AND colorobject LT colorstop)
oplot, xobject[zrange], yobject[zrange], ps=3, color=FSC_COLOR('dark orchid') ;Fake QSOs

The hope is that I will expand on this list every time I come across something I think I might use again (but don't use enough to memorize).

Commissioning Quasars

I am trying to answer these questions about the likelihood targets in the commissioning data. Today I took the collated targets for the QSO commissioning region from the wiki (here) and used spherematch to align it with the truth table that Nic Ross has provided me (based on visually inspected spectra). For truth table see the Nov 20th email from Ross with subject line BOSS Quasars "Truth Table." Both are also on riemann in the directory:

/home/jessica/boss/BOSS_QSOs_First14plates.fits
/home/jessica/boss/bosstarget-qso-comm-collate.fits

Target bitmasks from the wiki:

# restrictive qso selection
maskbits BOSS_TARGET1 10 QSO_CORE
# permissive qso selection
maskbits BOSS_TARGET1 11 QSO_BONUS
# known qso between [2.2,9.99]
maskbits BOSS_TARGET1 12 QSO_KNOWN_MIDZ
# known qso outside of miz range. never target
maskbits BOSS_TARGET1 13 QSO_KNOWN_LOHIZ
# Neural Net that match to sweeps/pass cuts
maskbits BOSS_TARGET1 14 QSO_NN
# UKIDSS stars that match sweeps/pass flag cuts
maskbits BOSS_TARGET1 15 QSO_UKIDSS
# kde targets from the stripe82 coadd
maskbits BOSS_TARGET1 16 QSO_KDE_COADD
# likelihood method
maskbits BOSS_TARGET1 17 QSO_LIKE

These are bitmasks, so in IDL to select items that are only selected by the likelihood method we would do the following:
LIKEONLY = where(truthmatch.BOSS_TARGET1 EQ ishft(long(1),17))

To select items which are selected by the likelihood (and others):
LIKE = where((truthmatch.BOSS_TARGET1 AND ishft(long(1),17)) GT 0)

To select items which are not selected by the likelihood (selected by others):
NOTLIKE = where((truthmatch.BOSS_TARGET1 AND ishft(long(1),17)) EQ 0)

To select objects with high redshift:
HIGHZ = where(truthmatch.Z_PERSON GT 3)

Below is a color-color plot of the above selections on top of the templates used for the likelihood selection. Red is everything. Green is QSOs. Magenta are objects that the likelihood method selected (LIKE). Cyan are objects that only the likelihood method selected (LIKEONLY). White are objects that have a high redshift (HIGHZ):

Color-Color (u-g vs g-r)
of likelihood selected objects

Below are objects that are QSOs but were missed by the likelihood (NOTLIKE). Need to figure out why.

Color-Color (u-g vs g-r)
of not-likelihood QSO objects

Pretty Plots

I got Joe's plotting code to work. Below are color-color plots of quasars where the "temperature" of the data points is the redshift of the quasars. They are plotted on top of a black/gray stellar locus. Now to plot these same objects where the temperature is a likelihood instead of redshift:

Pretty colors!

Plotting with Joe

Joe Hennawi and I worked on trying to figure out where exactly we are going wrong with the likelihood method. He suggested making a two-dimensional histogram of the likelihoods binned in color-color space. This would allow us to see where in color-color space were are getting objects which have high likelihoods and will help us debug if there is a population of objects which are being falsely selected.

It is quite easy to make the 2D histogram, but plotting it has proved to be an issue. I've spent most of the day trying to do this. Here is how I create the histogram:

qsos = testqsos

b_u = 1.4
b_g = 0.9
b_r = 1.2
b_i = 1.8
b_z = 7.4

z_temp = qsos.Z_TOT
fu1 = qsos.PSFFLUX[0]
fg1 = qsos.PSFFLUX[1]
fr1 = qsos.PSFFLUX[2]
fi1 = qsos.PSFFLUX[3]
fz1 = qsos.PSFFLUX[4]

ivar_u1 = qsos.PSFFLUX_IVAR[0]
ivar_g1 = qsos.PSFFLUX_IVAR[1]
ivar_r1 = qsos.PSFFLUX_IVAR[2]
ivar_i1 = qsos.PSFFLUX_IVAR[3]
ivar_z1 = qsos.PSFFLUX_IVAR[4]


u1 = sdss_flux2mags(fu1, b_u)
g1 = sdss_flux2mags(fg1, b_g)
r1 = sdss_flux2mags(fr1, b_r)
i1 = sdss_flux2mags(fi1, b_i)
z1 = sdss_flux2mags(fz1, b_z)

sigu = sdss_ivar2magerr(ivar_u1, fu1, b_u)
sigg = sdss_ivar2magerr(ivar_g1, fg1, b_g)
sigr = sdss_ivar2magerr(ivar_r1, fr1, b_r)
sigi = sdss_ivar2magerr(ivar_i1, fi1, b_i)
sigz = sdss_ivar2magerr(ivar_z1, fz1, b_z)

ug = u1-g1
gr = g1-r1
ri = r1-i1
iz = i1-z1

icut = where(i1 GT 19.5)

nx = 50
ny = 50
imageQSO = fltarr(nx,ny)
imageAll = fltarr(nx,ny)
imageRatio = fltarr(nx,ny)

likeQSO = qsoqlike[icut]
likeAll = qsoslike[icut]
likeRatio = likeQSO/likeAll

qsoug = ug[icut]
qsogr = gr[icut]

populate_image, imageQSO, qsoug, qsogr, weight = likeQSO
populate_image, imageALL, qsoug, qsogr, weight = likeAll
populate_image, imageRatio, qsoug, qsogr, weight = likeRatio

Now I need code to plot this histogram in color-color space with hotter colors in bins with higher likelihoods and cooler colors in bins with smaller likelihoods. Anyone have code that does this in IDL or python? Joe has some, but I'm having a hard time getting it to work.