;------------------------------------------------------------------------------ ;+ ; NAME: ; make_catalog_radec ; ; PURPOSE: ; Starting with an input catalog and mask, trim the known galaxies ; by the completeness, trim the new galaxies from low completeness ; polygons and generate Nran randoms which are passed through the veto mask. ; ; CALLING SEQUENCE: ; make_catalog,'catalog-file-name','mask-file-name-base','output-base', 'nbar-base', RoDfac, catalog_type, ransackseed, zseed ; nbarbase and zseed are not used here, but keep consistent calling sequence for make_catalog_radec and make_catalog_z. ; ; OPTIONAL INPUTS: ; ; OUTPUTS: ; Two fits files for the data and randoms. One txt file with nbar(z). ; ; COMMENTS: ; ; EXAMPLES: ; ; REQUIREMENTS: ; A bunch of stuff we should document here. ; ; REVISION HISTORY: ; 25-May-2011 Written by Martin White (LBNL) ; 6-Sep-2011 Update to is_in_window_pix. ; 6-Sep-2011 Put in vetos piece by piece. ; 20-Apr-2012 Take out the explicit inclusion of fiber collided galaxies, since we are doing this by up-weighting the good neighbor galaxy now, Shirley Ho ; 21-Apr-2012 Output a 'full galaxy catalog' which includes all 'galaxies' within mask (so it does not cut on imatch) ; 14-Sep-2012 BR: Adding catalog_type variable. Only difference here is how random redshifts are determined. ; xx-xx-2013 BR: Changed inputs in the following way. ; RoDfac: Instead of requiring a number of randoms, input the desired ratio of randoms to data. ; Std input is 51 to account for randoms that will be veoted. ; Two random seeds are used in the code; one to send to ransack to generate ra/dec, and one for this code ; to randomly assign redshifts to randoms from the data. Both should be unsigned ints. ; 18-03-2014 BR: Breaking make_catalog.pro into make_catalog_radec.pro and make_catalog_z.pro ;- pro make_catalog_radec,galaxy_file,mask_base,output_base,nbar_base,RoDfac,catalog_type,ransackseed,zseed print,' ' print,'make_catalog called with:' print,' galaxy_file =',galaxy_file print,' mask_base =',mask_base print,' output_base =',output_base print,' nbar_base =',nbar_base print,' Nran/Ndata =',RoDfac print,' catalog_type = ',catalog_type print,' ransackseed = ',ransackseed, long(ransackseed) print,' zseed = ',zseed, long(zseed) print,' note known redshift random seed is always 1234L!' print,' this assures you that you get the same catalog ' print,' every time you run with the same input.' print,' Note make_catalog_radec does not use nbar_base or zseed' print,' Calling sequence is same for make_catalog_radec and make_catalog_z' print,' ' print,'testing',comov_vol(0.43,0.7) ;; hard coding the redshift bins we'll use ;; also hard-coding P0! P0 = 20000.0 h = 0.7 z_s = 0.0 ; start of the redshift range for CMASS is 0.43, but we can assign fkp weights for all z. z_e = 1.0 ; end of the redshift range for CMASS is 0.7, but we can assign fkp weights for all z. dz = 0.005 ; reading off the figure 2 (dndz) plot from Anderson et al. zbin= floor((z_e-z_s)/dz)+1 total_sqdeg = 4*!PI*(180/!PI)*(180/!PI) ; in full sky print, 'Total square degree in full sky:',total_sqdeg output_fullbase = output_base+'-full' ; 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/") ; Start with the target list. print,'Opening galaxy file:',galaxy_file dd = mrdfits(data_dir+galaxy_file,1,hdr) print,'Read ',n_elements(dd),' galaxies from ',data_dir+galaxy_file ; Read the mask file. mask_file = data_dir+mask_base+'.fits' read_fits_polygons,mask_file,geom wt = geom[dd.ipoly].weight ; The completeness for each galaxy. dd.comp = geom[dd.ipoly].weight ; assign the completeness for each galaxy ; Subsample known redshifts and trim all galaxies by the mask. rr = randomu(1234L,n_elements(dd)) ww = where((wt GT 0.0) AND ((dd.imatch NE 2) OR (rr LT wt)),cnt) ddf = dd[ww] print,'Have ',cnt,' galaxies after geometry cuts.' ; Write out the data to a useful FITS file with all galaxies. print,'Writing data to :',output_fullbase+'.dat.fits' mwrfits,ddf,data_dir+output_fullbase+'.dat.fits',/create ; Subsample known redshifts and trim all galaxies by the mask. ww = where( (wt GT 0.0) AND $ ( (dd.imatch EQ 1) OR $ ;( (dd.imatch EQ 1) OR (dd.imatch EQ 3) OR $ ; NOT sure about imatch EQ 3, since we are doing this by weight now, shirley ho April 20,2012 ((dd.imatch EQ 2) AND (rr LT wt)) ), cnt ) dd = dd[ww] print,'Have ',cnt,' galaxies after geometry cuts.' ngal = n_elements(dd) ; subsample weights also for the random ; for 2, we don't need redshifts anyway, but just use 0 case anyway to keep code from crashing. if(catalog_type EQ 0 OR catalog_type EQ 2) then begin dd_weight = double(dd.weight_systot)*double(dd.weight_noz + dd.weight_cp -1.0) total_gals = total(double(dd.weight_noz + dd.weight_cp -1.0)) endif if(catalog_type EQ 1) then begin ; subsample weights also for the random ; only difference here is we're NOT including the fiber collision galaxies' NN redshifts ; in the sample to draw random redshifts from. ; this should leave D and R unclustered in the (1+w) method. dd_weight = double(dd.weight_systot)*double(dd.weight_noz) total_gals = total(double(dd.weight_noz)) endif dd_weight_cum = total(dd_weight,/cumulative) weight_max = max(dd_weight_cum) ;; note these sums won't be the same unless you make the vectors double! ;; print,'these should be the same: ',weight_max,total(dd_weight) ; Now generate randoms. These are written to "temp" for now. ; This could cause race problems if more than one copy of the code is ; running at once. mask_file = data_dir+mask_base+'.ply' ;; want different randoms every time. ; seed = '1234' Nran = RoDfac*total_gals print, 'setting Nran = ',Nran ransacktmpfile = data_dir+output_base+'.ransacktmp' print,'using this tmp file to avoid race problems: ',ransacktmpfile cmd = [filepath('ransack',root_dir=getenv('IDLUTILS_DIR'),subdir='bin'),$ '-c',strcompress(ransackseed,/rem),'-r',strcompress(Nran,/rem),mask_file,ransacktmpfile] print,cmd spawn,cmd,/noshell readcol,ransacktmpfile,ra,dec,ipoly,format='D,D,I' ;; tinker adds systematics fields to the structure myNran = n_elements(ra) ;; hopefully this fixed int/double troubles?? rancat = replicate($ create_struct('ra',0d,'dec',0d,'z',0d, 'weight_fkp',0d, 'nz', 0d,'ipoly',0,'isect',0, 'zindx',0L, $ 'skyflux',fltarr(5), 'image_depth',fltarr(5), 'airmass',0.0, 'eb_minus_v',0.0, $ 'psf_fwhm',fltarr(5)),myNran) print,myNran,Nran,n_elements(rancat),n_elements(ra) rancat.ra = ra rancat.dec = dec ;; moved redshift assignment below to grab nbar and fkp weight at same time. ;; for now juts fill in dummy value. rancat.z = 0.5 ;; these are assigned in veto.pro now. ;rancat.ipoly = ipoly ;rancat.isect = geom[ipoly].sector nrstart = n_elements(rancat) ; Pass the randoms through the veto mask. print,'Vetoing randoms...' veto,rancat nrend = n_elements(rancat) print,'Before and after veto masks this many:',nrstart,nrend vetofraction = float(nrend)/float(nrstart) print,'Fractional correction to area from veto regions: ',vetofraction ;; compute nbar in bins, write to a file. npoly = n_elements(geom) mask_area = dblarr(npoly) mask_area_eff = dblarr(npoly) for i=0L, npoly -1L do begin if (geom[i].weight gt 0.001) then begin mask_area(i) = garea(geom[i]) ; total masked area (in str) from all the polygons (regardless of redshift) mask_area_eff(i) = garea(geom[i])*geom[i].weight ; effective area covered endif else begin mask_area(i) = 0.0 endelse endfor mask_area_t = total(mask_area) mask_area_insqdeg = mask_area_t*(180/!PI)*(180/!PI) print,'Total unweighted area in mask:',mask_area_insqdeg mask_area_t = total(mask_area_eff) mask_area_insqdeg = mask_area_t*(180/!PI)*(180/!PI) print,'Total weighted area in mask:',mask_area_insqdeg ;; final weighted area accounts for veto regions. mask_area_insqdeg = mask_area_insqdeg*vetofraction print,'Final effective area: ',mask_area_insqdeg print,'Correponding number of galaxies in effective area: ',total_gals ; print nbar info to a text file. close,1 nbarfname = data_dir+nbar_base+'.dat' openw,1,nbarfname zcenlist = fltarr(zbin) zlowlist = fltarr(zbin) zhighlist = fltarr(zbin) n_bar = fltarr(zbin) w_fkp = fltarr(zbin) veff = fltarr(zbin) mask_vol = fltarr(zbin) cnt_tot = fltarr(zbin) X = dz*indgen(zbin)+z_s for i = 0, (zbin-1L) do begin zi = z_s + dz*i zn = z_s + dz*(i+1) gal_in_z1 = where ( (dd.z ge zi) and (dd.z lt zn), cnt) ;; original code forgot weights to compute nbar. ;; first term computes the probability of z in the range of interest. cnt1 = 0.0 if(cnt GT 0) then begin cnt1 = (total(dd_weight[gal_in_z1])/weight_max)*total_gals endif cnt_tot(i) = cnt1 tot_vol = comov_vol(zi,zn) ; Mpc^3 tot_vol = tot_vol*(h*h*h) ; (Mpc/h)^3 mask_vol(i) = tot_vol/total_sqdeg*mask_area_insqdeg ; masked volume in this redshift bin zcenlist(i) = 0.5*(zi+zn) zlowlist(i) = zi zhighlist(i) = zn n_bar(i) = cnt_tot(i)/mask_vol(i) w_fkp(i) = 1./(1.+ (n_bar(i)*P0)) veff(i) = ((n_bar(i)*P0)/(1+(n_bar(i)*P0)))^2.0*mask_vol(i) print,' ',zi,' ',zn,' ',n_bar(i),' ',w_fkp(i),' ',mask_vol(i),' ',cnt_tot(i) endfor total_eff_vol = total(veff) print,'total effective volume is',total_eff_vol printf,1,format='("# effective area (deg^2), effective volume (Mpc/h)^3: ", E, E)',mask_area_insqdeg,total_eff_vol printf,1,format='("# zcen,zlow,zhigh,nbar,wfkp,shell_vol,total weighted gals")' for i = 0, (zbin-1L) do begin ;;; original. now we'll compute nbar, wfkp, and cnttot after we do systematic weights, which now comes after this step. ; printf,1,format='(F, F, F, E, E, E, E)',zcenlist(i),zlowlist(i),zhighlist(i),n_bar(i),w_fkp(i),mask_vol(i),cnt_tot(i) printf,1,format='(F, F, F, E, E, E, E)',zcenlist(i),zlowlist(i),zhighlist(i),-1.0,-1.0,mask_vol(i),-1.0 endfor close,1 ;; temporary assignments, fill in later. dd.nz = -1.0 dd.weight_fkp = -1.0 ; Write out the data to a useful FITS file with all galaxies. print,'Writing data to :',output_base+'.dat.fits' mwrfits,dd,data_dir+output_base+'.dat.fits',/create rancat.nz = -1.0 rancat.weight_fkp = -1.0 ; Write out the randoms to a useful FITS file. print,'Writing randoms to:',output_base+'.ran.fits' mwrfits,rancat,data_dir+output_base+'.ran.fits',/create ; Clean up the temporary file we used. spawn,'rm '+ransacktmpfile end