;------------------------------------------------------------------------------ ;+ ; NAME: ; make_mask ; ; PURPOSE: ; Updates the boss_geometry file so that the weight in each polygon is ; the sector completeness, computed from the galaxy catalogs. ; ; CALLING SEQUENCE: ; make_mask,'catalog-file-name','mask-output-file-name-base',mincomp ; -- catalog-file-name as created by match_sample. ; -- File name base to which .ply or .fits will be appended. ; -- Minimum completeness sectors to keep. ; ; OPTIONAL INPUTS: ; ; OUTPUTS: ; FITS and polygon files containing the mask. ; ; COMMENTS: ; ; EXAMPLES: ; ; REQUIREMENTS: ; A bunch of stuff we should document here. ; ; REVISION HISTORY: ; 24-May-2011 Written by Martin White (LBNL) ; 20-Apr-2012 Edited by Shirley Ho ; 14-Sep-2012 Added different options for mask type. ; 0 -- Standard Anderson et al. LSS choice. ; 1 -- Completenesses appropriate for (1+w) small scale method. ; Fiber collided galaxies are treated as random losses that contribute ; to sector incompleteness. ; Mask still trimmed according to method 0, so that the footprints are ; the same [for the same mincomp, minzcomp]. ; 2 -- Binary mask -- Footprint determined in the same way as method 0, ; but completeness values are only 1 or 0. Useful for angular clustering measurements ; of the targets. ; 20-Feb-2013 Excluding chunks < 7 if lowzchunkcut = 1 ;- pro make_mask,galaxy_file,output_file_base,mincomp,minzcomp,catalog_type,lowzchunkcut print,' ' print,'Make mask.' print,'Called with galaxy_file =',galaxy_file print,'Called with output_file =',output_file_base print,'Called with mincomp =',mincomp print,'Called with minzcomp =',minzcomp PRINT,'Called with catalog_type=',catalog_type print,'called with lowzchunkcut=',lowzchunkcut print,' ' print,'All catalog_types return same footprint for a given input mincompWITHCP,minzcomp' print,'But different catalog_types will have different completeness values in each sector.' ; minzcomp = 0.8 according to redshift completeness set by Anderson et al. Below equation (13) ; Set up some useful paths for the code below. data_dir = filepath(root=getenv("BOSSLSS_DIR"),"/data/") geom_dir = filepath(root=getenv("BOSSLSS_DIR"),"/geometry/") ; Read the galaxies in. galaxy = mrdfits(data_dir+galaxy_file,1,hdr) print,'Read information on ',n_elements(galaxy),' galaxies from ',$ data_dir+galaxy_file ;; make file used by Manera mocks. close,1 marcfile = output_file_base + 'formarc.dat' openw,1,marcfile print, 'opened this file: ',marcfile ; ; Count all targets that have imatch = 0, 1, 3, 5 ; ; Ignore "old" galaxies and stars. ; wall = where((galaxy.imatch NE 2) AND (galaxy.imatch NE 4), cntall) ; ; Count "new" redshifts; since we're going to assign fiber collision ; ; galaxes redshifts, count them as having one. (??) ; wnew = where((galaxy.imatch EQ 1) OR (galaxy.imatch EQ 3), cntnew) ; print,'Have ',cntall,' non-known stars and galaxies.' ; print,'Have ',cntnew,' new BOSS galaxies with redshifts.' ; changes on april 20,2012 - shirley ho ; Sector completeness: (N_obs + N_cp)/(N_targ-N_known) following equation (10) of Anderson et al. ; -> (N_star+N_gal+N_fail +N_cp) / (N_star+N_gal +N_fail+N_cp+ N_missed + N_known - N_known) ; -> (N_star+N_gal+N_fail +N_cp)/ (N_star+N_gal +N_fail+N_cp+ N_missed) ; -> imatch -> im, when we see N_fail or N_cp, I count everything (even those we cannot correct for) ; -> [ im4 + im1 +(im5+im7) +(im3+im6)] / [ im4 + im1 +(im5+im7) +(im3+im6) + im0] ; so I guess sector completeness is really high (since it is only different by im0 -> N_missed) wnew = where( (galaxy.imatch NE 0) AND (galaxy.imatch ne 2), cntnew) wnewnocp = where( (galaxy.imatch NE 0) AND (galaxy.imatch NE 2) AND $ (galaxy.imatch NE 3) AND (galaxy.imatch NE 8), cntnewnocp) wall = where( (galaxy.imatch NE 2), cntall) print, 'Calculating sector completeness...' print, 'The denominator is N_targ - N_known, has ',cntall,' galaxies' print, 'The nominator is N_obs +N_cp, has ',cntnew,' galaxies' print, 'The nominator is N_obs, has ',cntnewnocp,' galaxies' ; Read in the geometry file we'll be updating. We need to read ; both the .ply and .fits files to get the pixel information. ;; just to be safe! polyhdr can't be overwritten by read_angle_polygons for some reason. polyhdr=0 read_mangle_polygons,geom_dir+'boss_geometry_2012_12_04.ply',dummy,hdr=polyhdr read_fits_polygons,geom_dir+'boss_geometry_2012_12_04.fits',bosspoly ;read_mangle_polygons,geom_dir+'boss_geometry_2011_06_10.ply',dummy,hdr=polyhdr ;read_fits_polygons,geom_dir+'boss_geometry_2011_06_10.fits',bosspoly ; Get the pixelisation scheme pixscheme = (strsplit(polyhdr[0],/extract))[1] ; Will die here unless uses simple pix scheme. pixres = nint((strsplit(pixscheme,'s',/extract))[0]) print,'Using ',pixres,'s pixels from boss_geometry' ; Work out the sector completeness. maxsect = max(bosspoly.sector) histall = histogram(galaxy[wall].isect,min=0,max=maxsect) histnew = histogram(galaxy[wnew].isect,min=0,max=maxsect) histnewnocp = histogram(galaxy[wnewnocp].isect,min=0,max=maxsect) comp = histnew/(histall+1.D-10) ;this one used to define footprint; as in Anderson et al. compnocp = histnewnocp/(histall+1.D-10) xxx = where(histall EQ 0,cntxxx) print, 'this many sectors have no targets, setting default completeness to 1.',cntxxx if(cntxxx GT 0) then begin comp[xxx] = 1.0 compnocp[xxx] = 1.0 endif print,'Std completeness ranges from ',min(comp),' to ',max(comp) print,'1+w completeness ranges from ',min(compnocp),' to ',max(compnocp) ww = where(comp LT mincomp,cnt) if (cnt GT 0) then begin comp[ww]=0. compnocp[ww]=0. endif print,'Clipped ',cnt,' sectors below ',mincomp ; Redshift completeness : following Equation (13) of Anderson et al ; Changes on april 20,2012 - shirley ho ; Redshift completeness: N_gal/ (N_obs - N_star) ; -> N_gal /( N_star+N_gal+N_fail - N_star) ; -> im1 / ( im4 + im1 +(im5+im7) - im4) ; -> im1/ (im1 + (im5 + im7) ) ; We drop sectors which has redshift completeness < 0.8 AND (N_obs - N_cp) > 10 ; (N_obs -N_cp) = N_star+N_gal+N_fail-N_cp = im4 + im1 +(im5+im7) - (im3+im6) > 10 ; Strictly speaking, some imatch = 8 objects could have had good redshifts from BOSS or KNOWN, and this doesn't track that. ; Will's code does the same thing. wznom = where( (galaxy.imatch EQ 1), cntnew) wzdenom = where( (galaxy.imatch EQ 1) OR $ (galaxy.imatch EQ 5) OR $ (galaxy.imatch EQ 7), cntall) ; Ngal is same as wznom Ngal = where( (galaxy.imatch EQ 1), cntngal) Nzfail = where( (galaxy.imatch EQ 5) OR (galaxy.imatch EQ 7), cntnzfail) Nstar = where( (galaxy.imatch EQ 4),cntnstar) Nspectot = where( (galaxy.imatch EQ 1) OR (galaxy.imatch EQ 4) OR $ (galaxy.imatch EQ 5) OR (galaxy.imatch EQ 7), cnttot) ;; needed for Manera mocks Nlegacy = where( (galaxy.imatch EQ 2), cntleg) Ncptot = where( (galaxy.imatch EQ 3) OR (galaxy.imatch EQ 8),cntcp) histwznom = histogram(galaxy[wznom ].isect,min=0,max=maxsect) histwzdenom = histogram(galaxy[wzdenom].isect,min=0,max=maxsect) zcomp = histwznom/(histwzdenom+1.D-10) ; for sector cuts, do same as Will. histNgal = histogram(galaxy[Ngal].isect,min=0,max=maxsect) histNstar = histogram(galaxy[Nstar].isect,min=0,max=maxsect) histNzfail = histogram(galaxy[Nzfail].isect,min=0,max=maxsect) histNspec = histogram(galaxy[Nspectot].isect,min=0,max=maxsect) histNlegacy = histogram(galaxy[Nlegacy].isect,min=0,max=maxsect) histNcp = histogram(galaxy[Ncptot].isect,min=0,max=maxsect) histtot = histogram(galaxy.isect,min=0,max=maxsect) ;; currently Will's code just uses Ngal, so maybe this is right. Nobs = histNgal ; we drop sectors with redshift completeness < 0.8 AND ; if the sector has less than 10 BOSS spectra (N_obs-N_cp > 10) print,'Redshift Completeness ranges from ',min(zcomp),' to ',max(zcomp) wwz = where((zcomp LT minzcomp) and (Nobs GT 10), cnt) if (cnt GT 0) then begin comp[wwz]=0. compnocp[wwz]=0. endif print,'Clipped ',cnt,' sectors below redshift completneess',minzcomp,$ ' and have more than 10 observations.' ;; Finally, do some cuts to remove disconnected regions with no observed galaxies. ;; Do a brute force for loop because I don't think in IDL! ra_m90poly = bosspoly.ramid - 90.0d0 yy = where(ra_m90poly LT 0.0) ra_m90poly[yy] += 360.0 yy = where(ra_m90poly LT 0.0,cnt) if(cnt GT 0) then begin print,'ra_m90poly bug' stop endif areasec = fltarr(maxsect+1L) ra_m90sec = fltarr(maxsect+1L) rasec = fltarr(maxsect+1L) decsec = fltarr(maxsect+1L) for seci=0L,maxsect do begin xx = where(bosspoly.sector EQ seci,cntxx) if(cntxx GT 0) then begin areasec[seci] = total(bosspoly[xx].str) if(areasec[seci] GT 0.0) then begin ra_m90sec[seci] = total(ra_m90poly[xx]*bosspoly[xx].str)/areasec[seci] decsec[seci] = total(bosspoly[xx].decmid*bosspoly[xx].str)/areasec[seci] endif else begin ra_m90sec[seci] = ra_m90poly[xx[0]] decsec[seci] = bosspoly[xx[0]].decmid endelse endif ;if there are any polygons in this sector endfor ;end loop over sectors to get area, ra, dec for each. ;; only keep empty sectors that are surrounded. ;; copied from Will's code. for seci=0L,maxsect do begin if((areasec[seci] GT 0.0) AND (histNspec[seci] EQ 0)) then begin cnt_pos_ra_nbrs = 0L cnt_neg_ra_nbrs = 0L cnt_pos_dec_nbrs = 0L cnt_neg_dec_nbrs = 0L ssid = where((sphdist(ra_m90sec[seci],decsec[seci],ra_m90sec,decsec,/degrees) LT 2.0) AND $ ;; old version; no cut on neighboring completeness. ;; (histNspec GT 0),cntssid) ;; new version; suggested in email on Jun 22 to cut on neighbor's completeness too. ;; this should get rid of isolated region Eyal pointed out: ;;RA-Dec box of [215, 216] (RA) and [54.2, 54.6] (Dec). (histNspec GT 0) AND (comp GE mincomp),cntssid) if(cntssid GT 0) then begin for ii=0L,n_elements(ssid)-1L do begin if((ra_m90sec[ssid[ii]] LT ra_m90sec[seci]) AND (ssid[ii] NE seci)) then begin cnt_neg_ra_nbrs += 1L endif if((ra_m90sec[ssid[ii]] GT ra_m90sec[seci]) AND (ssid[ii] NE seci)) then begin cnt_pos_ra_nbrs += 1L endif if((decsec[ssid[ii]] LT decsec[seci]) AND (ssid[ii] NE seci)) then begin cnt_neg_dec_nbrs += 1L endif if((decsec[ssid[ii]] GT decsec[seci]) AND (ssid[ii] NE seci)) then begin cnt_pos_dec_nbrs += 1L endif endfor ;;ii loop endif ;found a 2 degree neighbor with some spectra. if((cnt_neg_ra_nbrs EQ 0L) OR (cnt_pos_ra_nbrs EQ 0L) OR $ (cnt_neg_dec_nbrs EQ 0L) OR (cnt_pos_dec_nbrs EQ 0L)) then begin comp[seci] = 0. compnocp[seci] = 0. endif endif ;sector with no observed galaxies. endfor ;;original line in Will's code, we'll match that even if it' s a bug. ;; if(sector[sid].area*180.*180./(4.*pi)>0.1 && ;; for DR10, we'll switch to 0.1 deg^2 holeinstr = 0.1*!pi*!pi/180./180. print, 'hole size in str: ',holeinstr ;; remove large holes for seci=0L,maxsect do begin if((areasec[seci] GT holeinstr) AND (histNspec[seci] EQ 0)) then begin comp[seci] = 0. compnocp[seci] = 0. endif endfor ;; remove chunks below 7 for lowz if (lowzchunkcut EQ 1) then begin if(n_elements(bosspoly[0].chunkvector) NE 4) then begin print, 'beth assumed number of chunks is 4! Recode this intelligently!!' stop endif cc = where(((bosspoly.chunkvector[0] GT 0) AND (bosspoly.chunkvector[0] LT 7)) OR $ ((bosspoly.chunkvector[1] GT 0) AND (bosspoly.chunkvector[1] LT 7)) OR $ ((bosspoly.chunkvector[2] GT 0) AND (bosspoly.chunkvector[2] LT 7)) OR $ ((bosspoly.chunkvector[3] GT 0) AND (bosspoly.chunkvector[3] LT 7)),cntcc) if(cntcc GT 0) then begin comp[bosspoly[cc].sector] = 0. compnocp[bosspoly[cc].sector] = 0. print,'lowzchunkcut eliminates ',cntcc,' additional polygons.' endif endif ; lowchunkcut ; Set the polygon weights to the completeness. ; choose completeness to fill in based on catalog_type choice. if(catalog_type EQ 0) then begin bosspoly.weight = comp[bosspoly.sector] endif if(catalog_type EQ 1) then begin bosspoly.weight = compnocp[bosspoly.sector] endif if(catalog_type EQ 2) then begin compbinary = comp ;; just cut on a really small completeness to figure out what's > 0. qq = where(compbinary GT 0.001,cntb) if(cntb GT 0) then begin compbinary[qq] = 1. endif bosspoly.weight = compbinary[bosspoly.sector] endif ; Write out the mask as .ply and .FITS files. write_fits_polygons, data_dir+output_file_base+'.fits', bosspoly write_mangle_polygons, data_dir+output_file_base+'.ply' , bosspoly,$ hdr=polyhdr,/pixels ;; print for Manera mocks. for seci=0L,maxsect do begin xxcp = where(galaxy.isect EQ seci,cntxxcp) if(cntxxcp GT 0) then begin cp_comp = galaxy[xxcp[0]].cp_comp endif else begin cp_comp = -1.0 endelse ;print, seci,comp[seci],zcomp[seci],cp_comp printf,1,format='(I, F, F, F, I, I, I, I, I, I)',seci,comp[seci],zcomp[seci],cp_comp, $ histNgal[seci], histNstar[seci], histNzfail[seci], histNlegacy[seci], histNcp[seci], histtot[seci] endfor close,1 end