Wednesday, April 21, 2010

Known Quasars and Set Theory

It's been kind of a nightmare to put a file of all known quasars act as a truth table for my likelihood testing. Adam Myers has a list of all known quasars from SDSS DR 7 and other data sets, and then Nic Ross has a truth table with BOSS Quasars.

So here is code to combine them without repeats:

;Read in BOSS Truth table
collatefile = '/home/jessica/boss/BOSS_Quasars_3PCv3_cat.fits.gz'
collateinfo = mrdfits(collatefile, 1)

;Trim to only be QSOs (redshift > 2 with a high human confidence)
confcut = where(collateinfo.z_conf_person GT 2 and collateinfo.Z_PERSON GE 2)
confTT = collateinfo[confcut]

;Read in Myers known quasars
knownqsofile = '/home/jessica/boss/knownquasarstar.041310.fits'
knownqsos = mrdfits(knownqsofile, 1)

;Constrain to only be QSOS in the BOSS redshift range
knownqso = knownqsos[where(strmatch(knownqsos.source, "*QSO*"))]
knownbossqsos = knownqso[where((knownqsos.zem GE 2) AND (knownqsos.zem LE 5.0))]

;Spherematch these two sets to find overlapping qsos (qsooverlap)
spherematch, confTT.ra, confTT.dec, knownbossqsos.ra, knownbossqsos.dec, 2./3600, i1, qsooverlap, d12

;Subtract the overlapping QSOs from the Myers known QSOs
knownqsoindex = indgen(n_elements(knownbossqsos.dec))
knownqsos = setdifference(knownqsoindex, qsooverlap)

(the functions setdifference, setunion and setintersection can be found here)

;Combine the two sets to get the ra, dec, and redshift of all known QSOs
quasarsRa = [knownbossqsos[knownqsos].ra, confTT.ra]
quasarsDec = [knownbossqsos[knownqsos].dec, confTT.dec]
quasarsZ = [knownbossqsos[knownqsos].zem, confTT.z_person]

1 comment:

  1. An elegant way of removing these sorts of duplicates, I think:

    a=mrdfits('somefilewithdups',1)
    ;ADM 3.6 is an example, change depending on coordinate accuracy
    ;ADM maxmatch=0 records all instances of pairs
    spherematch, a.ra, a.dec, a.ra, a.dec, 3.6/3600, m1, m2, maxmatch=0
    ;ADM dups is unique indices of all duplicates
    dups = m1[where(m1 gt m2)]
    ;ADM an array set to 1 for everything...everything's good
    good = [indgen(n_elements(a))*0]+1
    ;ADM zero the array where there are duplicates
    good[dups]=0
    ;ADM discard the duplicates
    dupsremoved = a[where(good)]
    mwrfits, dupsremoved, 'somefilewithoutdups',/create

    So you could make one large file of all the known quasars and then remove all duplicates. This has the added bonus that it will remove duplicates in the 3PC catalog (There usually are some....). Plus, you can retrieve the duplicates:

    justthedups = a[where(good eq 0)]

    ReplyDelete