; 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:
No comments:
Post a Comment