;+ ; NAME: ; desi_quicksim ; ; PURPOSE: ; Simple S/N calculator for DESI spectra ; ; CALLING SEQUENCE: ; desi_quicksim, model= [, wave, flux, infile=, airmass=, exptime=, $ ; nread=, outfile=, simdat= ] ; ; INPUTS: ; model - Throughput model for light entering the fiber; ; 'lrg' for a galaxy with r_deV=1.0 arcsec ; 'elg' for a galaxy with r_exp=0.45 arcsec ; 'star' for point source ; 'perfect' for no losses ; ; OPTIONAL INPUTS: ; wave - Input spectrum wavelength vector [NPIX] ; flux - Input spectrum flux vector [NPIX] ; infile - Input ASCII file ; airmass - Airmass for atmospheric extinction; default to 1 ; exptime - Exposure time; default to value in parameter file ; nread - Number of reads for computing read noise; default to 1 ; (the only effect is to increase read noise by sqrt(NREAD)) ; outfile - Optional output file name; 'sn-'+INFILE if set as /OUTFILE ; ; OUTPUTS: ; ; OPTIONAL OUTPUTS: ; simdat - Output structure with simulated data ; ; COMMENTS: ; Either WAVE,FLUX must be set or INFILE for the input spectrum, ; where the wavelength is in units of Ang and the flux in units ; of 1e-17 erg/cm^2/s/Ang. This spectrum can have any wavelengh sampling, ; and will be linearly interpolated to the sky spectrum wavelengths. ; ; The sky spectrum file is required to have a linear wavelength scale. ; ; The following parameters are drawn from the data file 'desi.yaml': ; exptime: default exposure time [seconds] ; fiberdiam: fiber diameter on the sky [arcsec] ; M1_diameter: primary mirror diameter [meters] ; obscuration_diameter: central obscuration diameter [meters] ; M2_support_width: spider support widths [meters] ; The following parameters should be defined in each of the 3 cameras ; in order b,r,z: ; readnoise: read noise [electrons] ; ang_per_row: average Ang per row in the dispersion direction ; psf_fwhm_wave: average PSF FWHM in wavelength direction [pix] ; psf_fwhm_spatial: average PSF FWHM in spatial direction [pix] ; ; Full documentation is on the DESI project wiki here: ; https://desi.lbl.gov/trac/wiki/Pipeline/QuickSim ; ; EXAMPLES: ; ; BUGS: ; Descritization noise from the amplifier has been ignored, which ; would only be important if the read noise approaches 1 and the gain ; is larger than 1. ; ; The final spectrum output is hardwired as 5X binning of the input ; sky spectrum. ; ; DATA FILES: ; $DESIMODEL/data/desi.yaml ; $DESIMODEL/data/specpsf/psf-quicksim.dat ; $DESIMODEL/data/spectra/spec-sky.dat ; $DESIMODEL/data/spectra/ZenithExtinction-KPNO.dat ; $DESIMODEL/data/throughput/fiberloss-perfect.fits ; $DESIMODEL/data/throughput/fiberloss-star.fits ; $DESIMODEL/data/throughput/fiberloss-elg.fits ; $DESIMODEL/data/throughput/fiberloss-lrg.fits ; $DESIMODEL/data/throughput/thru-b.fits ; $DESIMODEL/data/throughput/thru-r.fits ; $DESIMODEL/data/throughput/thru-z.fits ; ; REVISION HISTORY: ; 10-Jun-2014 Written by D. Schlegel, LBL ;- ;------------------------------------------------------------------------------ pro desi_quicksim, wave_in, flux_in, infile=infile, airmass=airmass1, $ model=model, exptime=exptime1, nread=nread1, outfile=outfile1, simdat=simdat if (keyword_set(airmass1)) then airmass = airmass1 $ else airmass = 1.0 outfile = 0 if (size(outfile1,/tname) EQ 'STRING') then outfile = outfile1 $ else if (keyword_set(outfile1) AND keyword_set(infile)) $ then outfile = 'sn-'+file_basename(infile) if (keyword_set(nread1)) then nread = nread1 $ else nread = 1 topdir = getenv('DESIMODEL') if (NOT keyword_set(topdir)) then $ message, 'Enviroment DESIMODEL must be set' ; Read the instrument parameter file ; Compute the effective area in m^2 instfile = filepath('desi.yaml', root_dir=topdir, subdir='data') readcol, instfile, pname, pval, format='(a,f)', /silent fiberdiam = pval[(where(pname EQ 'diameter_arcsec:'))[0]] fibarea = !pi * (fiberdiam/2)^2 readnoise = pval[where(pname EQ 'readnoise:')] * sqrt(nread) geomarea = pval[(where(pname EQ 'geometric_area:'))[0]] ; Read the PSF parameters bpsf = mrdfits(topdir+'/data/specpsf/psf-quicksim.fits', 'QUICKSIM-B') rpsf = mrdfits(topdir+'/data/specpsf/psf-quicksim.fits', 'QUICKSIM-R') zpsf = mrdfits(topdir+'/data/specpsf/psf-quicksim.fits', 'QUICKSIM-Z') if (keyword_set(exptime1)) then exptime = exptime1 $ else exptime = pval[(where(pname EQ 'exptime:'))[0]] ; Read the sky spectrum and the atmospheric extinction model ; Everything will be reampled to the sky spectrum pixels readcol, filepath('spec-sky.dat', root_dir=topdir, $ subdir=['data','spectra']), wave, sky, format='(d,d)', /silent nwave = n_elements(wave) readcol, filepath('ZenithExtinction-KPNO.dat', root_dir=topdir, $ subdir=['data','spectra']), wave_tmp, extinct_tmp, format='(d,d)', /silent extinct = interpol(extinct_tmp, wave_tmp, wave) ; Read the input spectrum if (keyword_set(wave_in)*keyword_set(flux_in)) then begin flux = interpol(flux_in, wave_in, wave) endif else begin if (keyword_set(infile) EQ 0) then message, 'INFILE must be set' readcol, infile, wave_tmp, flux_tmp, format='(d,d)', /silent flux = interpol(flux_tmp, wave_tmp, wave) wave_tmp = 0 ; clear memory flux_tmp = 0 ; clear memory endelse ; Read the telescope+instrument throughput for the sky and object bthru = mrdfits(topdir+'/data/throughput/thru-b.fits', 'THROUGHPUT', /silent) rthru = mrdfits(topdir+'/data/throughput/thru-r.fits', 'THROUGHPUT', /silent) zthru = mrdfits(topdir+'/data/throughput/thru-z.fits', 'THROUGHPUT', /silent) thru_b = interpol(bthru.throughput, bthru.wavelength, wave) $ * (wave GT min(bthru.wavelength)) * (wave LT max(bthru.wavelength)) thru_r = interpol(rthru.throughput, rthru.wavelength, wave) $ * (wave GT min(rthru.wavelength)) * (wave LT max(rthru.wavelength)) thru_z = interpol(zthru.throughput, zthru.wavelength, wave) $ * (wave GT min(zthru.wavelength)) * (wave LT max(zthru.wavelength)) ; Account for fiber aperture losses lossfile = filepath('fiberloss-'+model+'.dat', root_dir=topdir, $ subdir=['data','throughput']) readcol, lossfile, wave_tmp, fibfac_tmp, format='(d,d)', /silent if (NOT keyword_set(wave_tmp)) then $ message, 'Error reading file '+lossfile fibfac = interpol(fibfac_tmp, wave_tmp, wave) ; Convert the object and sky counts to photons detected erg_per_photon = 6.62606957e-27 * 2.99792458e10 / (wave * 1e-8) dwave = shift(wave,-1) - wave dwave[nwave-1] = dwave[nwave-2] photon_per_flux_b = 1e-17 * exptime * (geomarea*100.^2) * thru_b $ * dwave / erg_per_photon photon_per_flux_r = 1e-17 * exptime * (geomarea*100.^2) * thru_r $ * dwave / erg_per_photon photon_per_flux_z = 1e-17 * exptime * (geomarea*100.^2) * thru_z $ * dwave / erg_per_photon nphoton_obj_b = flux * photon_per_flux_b * fibfac nphoton_obj_r = flux * photon_per_flux_r * fibfac nphoton_obj_z = flux * photon_per_flux_z * fibfac nphoton_sky_b = sky * fibarea * photon_per_flux_b nphoton_sky_r = sky * fibarea * photon_per_flux_r nphoton_sky_z = sky * fibarea * photon_per_flux_z ; Apply atmospheric absorption to the object but not the sky nphoton_obj_b *= 10^(-extinct/2.5) nphoton_obj_r *= 10^(-extinct/2.5) nphoton_obj_z *= 10^(-extinct/2.5) ; Smooth by the PSF of the instrument nhalf = 50L ; half-width for PSF convolution flux_smooth_b = fltarr(nwave) flux_smooth_r = fltarr(nwave) flux_smooth_z = fltarr(nwave) nphoton_obj_smooth_b = fltarr(nwave) nphoton_obj_smooth_r = fltarr(nwave) nphoton_obj_smooth_z = fltarr(nwave) nphoton_sky_smooth_b = fltarr(nwave) nphoton_sky_smooth_r = fltarr(nwave) nphoton_sky_smooth_z = fltarr(nwave) psf_sigma_b = interpol(bpsf.fwhm_wave, bpsf.wavelength, wave) / 2.35482 psf_sigma_r = interpol(rpsf.fwhm_wave, rpsf.wavelength, wave) / 2.35482 psf_sigma_z = interpol(zpsf.fwhm_wave, zpsf.wavelength, wave) / 2.35482 for i=nhalf+0L, nwave-nhalf-1L do begin thispsf_b = $ exp(-0.5*((wave[i-nhalf:i+nhalf]-wave[i])/psf_sigma_b[i])^2) thispsf_b /= total(thispsf_b) ; normalize the smoothing kernel flux_smooth_b[i] = total(flux[i-nhalf:i+nhalf]*thispsf_b) nphoton_obj_smooth_b[i] = total(nphoton_obj_b[i-nhalf:i+nhalf]*thispsf_b) nphoton_sky_smooth_b[i] = total(nphoton_sky_b[i-nhalf:i+nhalf]*thispsf_b) thispsf_r = $ exp(-0.5*((wave[i-nhalf:i+nhalf]-wave[i])/psf_sigma_r[i])^2) thispsf_r /= total(thispsf_r) ; normalize the smoothing kernel flux_smooth_r[i] = total(flux[i-nhalf:i+nhalf]*thispsf_r) nphoton_obj_smooth_r[i] = total(nphoton_obj_r[i-nhalf:i+nhalf]*thispsf_r) nphoton_sky_smooth_r[i] = total(nphoton_sky_r[i-nhalf:i+nhalf]*thispsf_r) thispsf_z = $ exp(-0.5*((wave[i-nhalf:i+nhalf]-wave[i])/psf_sigma_z[i])^2) thispsf_z /= total(thispsf_z) ; normalize the smoothing kernel flux_smooth_z[i] = total(flux[i-nhalf:i+nhalf]*thispsf_z) nphoton_obj_smooth_z[i] = total(nphoton_obj_z[i-nhalf:i+nhalf]*thispsf_z) nphoton_sky_smooth_z[i] = total(nphoton_sky_z[i-nhalf:i+nhalf]*thispsf_z) endfor ; Bin to the spectra to 0.5 Ang bins, with the centering at ; the same center as the high-resolution sampling binfac = 5 nnew = nwave / binfac wave_bin = dblarr(nnew) flux_bin_b = fltarr(nnew) flux_bin_r = fltarr(nnew) flux_bin_z = fltarr(nnew) nphoton_obj_bin_b = fltarr(nnew) nphoton_obj_bin_r = fltarr(nnew) nphoton_obj_bin_z = fltarr(nnew) nphoton_sky_bin_b = fltarr(nnew) nphoton_sky_bin_r = fltarr(nnew) nphoton_sky_bin_z = fltarr(nnew) psf_fwhm_bin_b = fltarr(nnew) psf_fwhm_bin_r = fltarr(nnew) psf_fwhm_bin_z = fltarr(nnew) for i=0L, nnew-1L do begin i1 = (i*binfac - (binfac-1)/2) > 0L i2 = (i*binfac + (binfac-1)/2) < (nwave-1L) wave_bin[i] = wave[i*binfac] flux_bin_b[i] = mean(flux_smooth_b[i1:i2]) flux_bin_r[i] = mean(flux_smooth_r[i1:i2]) flux_bin_z[i] = mean(flux_smooth_z[i1:i2]) nphoton_obj_bin_b[i] = total(nphoton_obj_smooth_b[i1:i2]) nphoton_obj_bin_r[i] = total(nphoton_obj_smooth_r[i1:i2]) nphoton_obj_bin_z[i] = total(nphoton_obj_smooth_z[i1:i2]) nphoton_sky_bin_b[i] = total(nphoton_sky_smooth_b[i1:i2]) nphoton_sky_bin_r[i] = total(nphoton_sky_smooth_r[i1:i2]) nphoton_sky_bin_z[i] = total(nphoton_sky_smooth_z[i1:i2]) psf_fwhm_bin_b[i] = mean(psf_sigma_b[i1:i2]) * 2.35482 psf_fwhm_bin_r[i] = mean(psf_sigma_r[i1:i2]) * 2.35482 psf_fwhm_bin_z[i] = mean(psf_sigma_z[i1:i2]) * 2.35482 endfor dwave_bin = shift(wave_bin,-1) - wave_bin dwave_bin[nnew-1] = dwave_bin[nnew-2] ; Read noise has an effective number of pixels for optimal extraction, ; which is equal to sqrt(4*!pi)*sigma for a Gaussian in the limit that ; it is well-sampled. This is the best-case where the trace location ; and PSF model is known. ; The terms with throughput optimally weight the read noise contribution ; from each camera, such that it increases by sqrt(2) at the dichroic ; crossover between cameras. rdnoise_bin_b = readnoise[0] * sqrt(dwave_bin $ * interpol(bpsf.neff_spatial / bpsf.angstroms_per_row, $ bpsf.wavelength, wave_bin)) rdnoise_bin_r = readnoise[1] * sqrt(dwave_bin $ * interpol(rpsf.neff_spatial / rpsf.angstroms_per_row, $ rpsf.wavelength, wave_bin)) rdnoise_bin_z = readnoise[2] * sqrt(dwave_bin $ * interpol(zpsf.neff_spatial / zpsf.angstroms_per_row, $ zpsf.wavelength, wave_bin)) ; The noise in the S/N has components from object Poisson noise, ; sky noise, and read noise (combining two cameras in the dichroic overlap) qq_b = nphoton_obj_bin_b GT 0 $ AND wave_bin GT min(bpsf.wavelength) AND wave_bin LT max(bpsf.wavelength) qq_r = nphoton_obj_bin_r GT 0 $ AND wave_bin GT min(rpsf.wavelength) AND wave_bin LT max(rpsf.wavelength) qq_z = nphoton_obj_bin_z GT 0 $ AND wave_bin GT min(zpsf.wavelength) AND wave_bin LT max(zpsf.wavelength) ; Add the S/N in quadrature between each of the 3 cameras sn2_b = qq_b * nphoton_obj_bin_b^2 / $ ((nphoton_obj_bin_b + nphoton_sky_bin_b + rdnoise_bin_b^2 + (qq_b EQ 0))) sn2_r = qq_r * nphoton_obj_bin_r^2 / $ ((nphoton_obj_bin_r + nphoton_sky_bin_r + rdnoise_bin_r^2 + (qq_r EQ 0))) sn2_z = qq_z * nphoton_obj_bin_z^2 / $ ((nphoton_obj_bin_z + nphoton_sky_bin_z + rdnoise_bin_z^2 + (qq_z EQ 0))) snvec_bin = sqrt(sn2_b + sn2_r + sn2_z) ; In dichroic overlaps, Average the PSF from the cameras weighting by the ; throughput rather than the (S/N)^2; as was done for the object spectrum; ; set to 1 when out-of-bounds qgood = (thru_b + thru_r + thru_z) GT 0 psf_fwhm_bin = (thru_b * psf_fwhm_bin_b + thru_r * psf_fwhm_bin_r $ + thru_z * psf_fwhm_bin_z) / (thru_b + thru_r + thru_z + (qgood EQ 0)) $ + (qgood EQ 0) ; Approximate the output flux spectrum as throughput-weighted between ; the channels. flux_bin = (thru_b * flux_bin_b + thru_r * flux_bin_r $ + thru_z * flux_bin_z) / (thru_b + thru_r + thru_z + (qgood EQ 0)) ; The below should be defined even in the case of no object flux!??? invvar_bin = snvec_bin^2 / (flux_bin^2 + (flux_bin LE 0)) ;err1=1./sqrt(invvar_bin) ;ibin = lindgen(nnew)*binfac ;err2= sqrt(nphoton_obj_bin_b + nphoton_sky_bin_b + rdnoise_bin_b^2) $ ; /(photon_per_flux_b[ibin]*fibfac[ibin]) ; Write output file if (keyword_set(outfile)) then begin openw, olun, outfile, /get_lun printf, olun, '# INFILE= '+infile printf, olun, '# AIRMASS= '+strtrim(airmass,2) printf, olun, '# MODEL= '+model printf, olun, '# EXPTIME= '+strtrim(exptime,2) printf, olun, '#' printf, olun, '# Median (S/N)= '+strtrim(median(snvec_bin),2) printf, olun, '# Total (S/N)^2= '+strtrim(total(snvec_bin^2),2) printf, olun, '#' printf, olun, '# Wave Flux Invvar S/N Counts_obj Counts_sky Counts_read FWHM' printf, olun, '# [Ang] [e-17 erg/s/cm^2/Ang] [1/(e-17 erg/s/cm^2/Ang)^2] [] [electrons] [electrons] [electrons] [Ang]' for i=0L, n_elements(wave_bin)-1L do $ printf, olun, wave_bin[i], flux_bin[i], invvar_bin[i], snvec_bin[i], $ nphoton_obj_bin_b[i]+nphoton_obj_bin_r[i]+nphoton_obj_bin_z[i], $ nphoton_sky_bin_b[i]+nphoton_sky_bin_r[i]+nphoton_sky_bin_z[i], $ rdnoise_bin_b[i]+rdnoise_bin_r[i]+rdnoise_bin_z[i], $ psf_fwhm_bin[i], format='(f9.2,6e12.4,f8.3)' close, olun free_lun, olun endif ; Construct the output data structure if (arg_present(simdat)) then begin simdat = create_struct( $ 'wavelength' , wave_bin, $ 'flux' , flux_bin, $ 'invvar' , invvar_bin, $ 'snvec' , snvec_bin, $ 'counts_obj' , nphoton_obj_bin_b+nphoton_obj_bin_r+nphoton_obj_bin_z, $ 'counts_sky' , nphoton_sky_bin_b+nphoton_sky_bin_r+nphoton_sky_bin_z, $ 'counts_read', rdnoise_bin_b+rdnoise_bin_r+rdnoise_bin_z, $ 'psf_fwhm' , psf_fwhm_bin) endif return end ;------------------------------------------------------------------------------