;+ ; NAME: ; fiberc_zfail ; ; PURPOSE: ; Start with the matched file , and deal with fiber collisions and redshift failures. ; ; CALLING SEQUENCE: ; fiberc_zfail,'catalog-file-name','output-file-name' ; ; OPTIONAL INPUTS: ; ; OUTPUTS: ; A fits file containing the matched sample with fiber collision corrections and redshift failures corrections. ; ; COMMENTS: ; Fills in the "imatch", "z", and "zpipe" and "zcom" entries of the structure. ; imatch==0: Not matched to a redshift. ; imatch==1: Matched to a new BOSS redshift. ; imatch==2: Matched to an old SDSS redshift. ; imatch==3: Fiber collision corrected -- the bit we add here! We also assign nearest neighbor redshifts to these objects, for future reference. ; imatch==4: Star. We want to remove it from the target list when ; computing completeness. ; imatch==5: Redshift failure corrected -- we do this here too! We assign nearest neighbor redshifts to these objects, for future reference. ; imatch==6: Fiber collision,initially flagged, and then we cannot ; find a galaxy* redshift for it ; (even when we use nearest neighbor up to ; radius_ratio*62 arcsec). ; imatch==7: Redshift failure,initially flagged and then we cannot ; find a galaxy* redshift for it ; (even when we use nearest neighbor up to ; radius_ratio*62 arcsec). ; imatch==8: had a good redshift, but disgarded a close pair containing a ; legacy galaxy to enforce close pair uniformity. ; ; EXAMPLES: ; ; REQUIREMENTS: ; A bunch of stuff we should document here. ; ; REVISION HISTORY: ; 24-Sep-2011 Written by Shirley Ho (LBNL) ; 27-Mar-2012 Some tidying up. ; 20-Apr-2012 Some more tidying up. ; 11-May-2012 BR: implemented sector-by-sector algorithm, added imatch = 8 ; to indicate pairs thrown out due to fiber collisions with ; legacy galaxies. ; 4-June-2012 BR: changed zfail corrections to be applied in the same sector. ; 14-June-2012 BR: added comments to make code more readable. Hope that helps! ; 5-Sept-2012 BR: fixed bug that was assigning fiber collision weight to imatch=8 galaxies ; near line 220 (where legacy close pairs are being downsampled). ; 13-Sept-2012 BR: filling in nearest neighbor redshifts for fiber-collided and redshift ; failure galaxies. BEWARE! ; pro fiberc_zfail,match_file,output_file print,' ' print,'Starting fiberc_zfail with input ',match_file ; need a random number generator to randomly throw out resolved fiber collisions with legacy galaxies. ;; removing the next line will get a random seed from systime or something. ;; instead, we fix it so the output of this code is deterministic for fixed ;; input files. blahseed = 5678L xx = randomu(blahseed) print,'first random',xx print,'the seed',blahseed data_dir = filepath(root=getenv("BOSSLSS_DIR"),"/data/") file_open = data_dir + match_file print,'Opening matched file:',file_open ee = mrdfits(file_open,1,hdr) hh = histogram(ee.imatch,min=0,max=7,binsize=1) print,'Before correcting for anything:' print,hh ; current correction (almost) precisely follows the method in Anderson et al. catalogs ; fiber collision corrections are handled on a sector-by-sector basis. cnt_legacy_removed = 0L ;; loop through all sectors. for si=min(ee.isect), max(ee.isect) do begin secgalsi = where(ee.isect EQ si, countgal) ;; get all the galaxies in this sector. ;; note we do this slightly differently than Will ;; he uses the known target flag, and we just use imatch = 2. ;; ;; get a list of the galaxies in the sector of interest that aren't imatch = 2, an old SDSS redshift. secgalsiUNKNOWN = where((ee.isect EQ si) AND $ (ee.imatch NE 2) , countunknown) nplates = 0L ;; number of plates covering the sector. platelist = intarr(10) ;; holds plate id's so we can identify unique plates covering the sector. platelist[*] = -10L ;; there are galaxies in this sector, do something. if(countgal GT 0) then begin ;; this for loop figures out how many plates cover the sector of interest. for itmp=0L, n_elements(secgalsi)-1L do begin gindx = secgalsi[itmp] found = 0L plateid = ee[gindx].plate if(plateid GT -10L) then begin ;-10L signifies no plate. this object got no BOSS fiber. for ip=0,nplates-1L do begin if(platelist[ip] EQ plateid) then begin found = 1L endif endfor ; loop over plates found in this sector already. if(found EQ 0L) then begin platelist[nplates] = plateid nplates = nplates + 1 endif ;found EQ 0 endif ;end if plate GT -10L endfor ;end loop over galaxies in sector to get nplate. ; new fiber collision code, sector by sector. ; for now we just do pair-by-pair calculation as in Will's catalog ;; maybe later want to consider larger collision groups than pairs. ; no longer using ingroup flag. ; nothing to do unless nplate > 0! if(nplates GT 0L) then begin ;;; first find nearest neighbor within 62'', then copy Will's code. ;;; NOTE: NOT CURRENTLY USED FOR ANYTHING! ingroup = spheregroup(ee[secgalsiUNKNOWN].ra, ee[secgalsiUNKNOWN].dec,$ (62./3600.),multgroup=mult,firstgroup=first,nextgroup=next) for qw = 0L, n_elements(mult)-1L do begin if(mult[qw] GT 1) then begin indx = first[qw] for xxx = 0L, (mult[qw] - 1L) do begin jndx = next[indx] for yyy = xxx+1L, (mult[qw] - 1L) do begin doubled = sphdist(ee[secgalsiUNKNOWN[indx]].ra, ee[secgalsiUNKNOWN[indx]].dec, ee[secgalsiUNKNOWN[jndx]].ra, ee[secgalsiUNKNOWN[jndx]].dec,/degrees)/62.*3600. jndx = next[jndx] endfor ; yyy indx = next[indx] endfor ; xxx endif ; mult endfor ;; now fill in counts of pairs and good and bad galaxies by looping over pairs. ;; loop over all pairs < 62'' ONCE. only look ahead of you. ;; when you find a fiber collided galaxy, assign it's weight to its closest angular neighbor. ;; also keep track of various pair statistics for use later. ;; make this faster using groups we already calculated. pair_in_sector_good = 0.0 pair_in_sector_tot = 0.0 gal_in_sector_good = 0.0 gal_in_sector_bad = 0.0 for groupid=0L, (n_elements(mult)-1L) do begin if(mult[groupid] LE 1) then continue ; no pairs in this group. indx = first[groupid] for ii = 0L, (mult[groupid]-1L) do begin ra_main = ee[secgalsiUNKNOWN[indx]].ra dec_main = ee[secgalsiUNKNOWN[indx]].dec imatchi = ee[secgalsiUNKNOWN[indx]].imatch jndx = next[indx] for jj = (ii+1L), (mult[groupid]-1L) do begin ra_cur = ee[secgalsiUNKNOWN[jndx]].ra dec_cur = ee[secgalsiUNKNOWN[jndx]].dec imatchj = ee[secgalsiUNKNOWN[jndx]].imatch dist_cur = sphdist(ra_cur,dec_cur,ra_main,dec_main,/degrees) if((dist_cur gt 0.0) and (dist_cur le (62./3600.))) then begin ;; at least one object had fiber or is close-pair corrected if((imatchi ne 0) or (imatchj ne 0)) then begin pair_in_sector_tot += 1.0 endif if((imatchi ne 0) and (imatchi ne 3) and $ (imatchj ne 0) and (imatchj ne 3)) then begin gal_in_sector_good += 1.0 pair_in_sector_good += 1.0 endif ;; one fiber and one with no fiber. if((imatchi eq 0) and (imatchj ne 0) and (imatchj ne 3)) then begin ee[secgalsiUNKNOWN[indx]].imatch = 3 ee[secgalsiUNKNOWN[indx]].z = ee[secgalsiUNKNOWN[jndx]].z ;redshift of NN imatchi = 3 ; this needs to update too so the same galaxy doesn't get updated twice. gal_in_sector_bad += 1.0 ee[secgalsiUNKNOWN[jndx]].weight_cp += 1L endif ;; one fiber and one with no fiber. if((imatchj eq 0) and (imatchi ne 0) and (imatchi ne 3)) then begin ee[secgalsiUNKNOWN[jndx]].imatch = 3 ee[secgalsiUNKNOWN[jndx]].z = ee[secgalsiUNKNOWN[indx]].z ;redshift of NN imatchj = 3 ; this needs to update too so the same galaxy doesn't get updated twice. gal_in_sector_bad += 1.0 ee[secgalsiUNKNOWN[indx]].weight_cp += 1L endif endif ;; within 62'' jndx = next[jndx] endfor ;; loop over jj indx = next[indx] endfor ;;loop over ii endfor ;; loop over groups in the sector. ;; default, when nplate > 1 and there are no pairs in current sector. cp_pair_over_poss = 1.0 cp_gal_over_poss = 1.0 if(pair_in_sector_tot gt 0.0) then begin cp_pair_over_poss = pair_in_sector_good/pair_in_sector_tot endif if((gal_in_sector_good+gal_in_sector_bad) gt 0.0) then begin cp_gal_over_poss = gal_in_sector_good / (gal_in_sector_good + gal_in_sector_bad) endif if(nplates EQ 1L) then begin cp_pair_over_poss=0.0 cp_gal_over_poss = 0.0 endif ;; now look for LEGACY close pairs in this sector. ;; remove some according to probability cp_gal_over_poss ingroup = spheregroup(ee[secgalsi].ra, ee[secgalsi].dec,$ (62./3600.),multgroup=mult,firstgroup=first,nextgroup=next) ;; careful! RNG is wacky. ;; initialize with an undefined variable to get a random seed from systime. for groupid=0L, (n_elements(mult)-1L) do begin if(mult[groupid] LE 1) then continue ; no pairs in this group. indx = first[groupid] for ii = 0L, (mult[groupid]-1L) do begin ra_main = ee[secgalsi[indx]].ra dec_main = ee[secgalsi[indx]].dec imatchi = ee[secgalsi[indx]].imatch jndx = next[indx] for jj = (ii+1L), (mult[groupid]-1L) do begin ra_cur = ee[secgalsi[jndx]].ra dec_cur = ee[secgalsi[jndx]].dec imatchj = ee[secgalsi[jndx]].imatch dist_cur = sphdist(ra_cur,dec_cur,ra_main,dec_main,/degrees) if((dist_cur gt 0.0) and (dist_cur le (62./3600.))) then begin if(((imatchi eq 2) or (imatchj eq 2)) and (imatchi ne 3) and (imatchj ne 3) and $ (imatchi ne 8) and (imatchj ne 8)) then begin r1 = randomu(blahseed) if(r1 GT cp_gal_over_poss) then begin r2 = randomu(blahseed) if(r2 GT 0.5) then begin ee[secgalsi[indx]].imatch = 8 ee[secgalsi[indx]].z = ee[secgalsi[jndx]].z ee[secgalsi[jndx]].weight_cp += 1 endif else begin ee[secgalsi[jndx]].imatch = 8 ee[secgalsi[jndx]].z = ee[secgalsi[indx]].z ee[secgalsi[indx]].weight_cp += 1 endelse cnt_legacy_removed += 1L endif ;; remove a galaxy! endif ;;at least one legacy, no BOSS fiber collision corrections. endif ;; within 62.'' endfor ;;jj loop endfor ;; ii loop endfor ;; groupid loop for xx=0L,(n_elements(secgalsi)-1L) do begin ee[secgalsi[xx]].cp_comp = cp_pair_over_poss endfor ;; new algorithm (in agreement with Will's DR9 algorithm); redshift failures ;; corrected on a sector-by-sector basis. ;; go serially through objects in sector to find closest relevant object. ;; only need to consider UNKNOWN galaxies (so old SDSS redshifts don't receive upweights from redshift failures. zz = where(ee[secgalsiUNKNOWN].imatch EQ 7, countzz) if(countzz GT 0) then begin for zzi = 0L,n_elements(zz)-1L do begin raz = ee[secgalsiUNKNOWN[zz[zzi]]].ra decz = ee[secgalsiUNKNOWN[zz[zzi]]].dec sepmin = 10000.0d0 nbrindx = -1L for znbr = 0L,n_elements(secgalsiUNKNOWN)-1L do begin ;;this won't allow zz[zzi] = znbr, since znbr has imatch = 7. if(ee[secgalsiUNKNOWN[znbr]].imatch EQ 1 or ee[secgalsiUNKNOWN[znbr]].imatch eq 4) then begin ra2 = ee[secgalsiUNKNOWN[znbr]].ra dec2 = ee[secgalsiUNKNOWN[znbr]].dec mysep = sphdist(raz,decz,ra2,dec2,/degrees) if(mysep LT sepmin) then begin sepmin = mysep nbrindx = znbr endif ;LT sepmin endif ;; good imatch endfor ;znbr loop if(nbrindx GT -1L) then begin ee[secgalsiUNKNOWN[zz[zzi]]].imatch = 5 ee[secgalsiUNKNOWN[zz[zzi]]].z = ee[secgalsiUNKNOWN[nbrindx]].z ; record redshift ee[secgalsiUNKNOWN[nbrindx]].weight_noz = $ ee[secgalsiUNKNOWN[nbrindx]].weight_noz + ee[secgalsiUNKNOWN[zz[zzi]]].weight_cp endif else begin ;; good nbrindx print, 'found no neighbor for this guy.' print, raz, decz ;stop endelse endfor ;zz loop over redshift failuers in the sector. endif ;;if any redshift failures in this sector. endif ; nplates GT 0 in this sector. endif ;countgal GT 0 endfor ;loop over sectors. ;; good! agrees with previous calculations of nplate. hh = histogram(ee.imatch,min=0,max=8) print,'After assigning bad-z galaxies redshifts imatch histogram is:' print,hh ; Write out the data to a useful FITS file. print,'Writing:',output_file mwrfits,ee,data_dir+output_file,/create print,'Done fiberc_zfail.' end