diff --git a/py/desisim/data/_quickcat.yaml b/py/desisim/data/_quickcat.yaml new file mode 100644 index 000000000..658019ea1 --- /dev/null +++ b/py/desisim/data/_quickcat.yaml @@ -0,0 +1,34 @@ +ELG: + EFFICIENCY: + SIGMA_FUDGE: 2.213938897168332 + SNR_CONTINUUM_SCALE: 7.202754668765929 + SNR_LINES_SCALE: 0.9651499811981232 + FAILURE_RATE: 0.006540594340964026 + UNCERTAINTY: + POWER_LAW_INDEX: 0.962171322 + SIGMA_17: 0.000553701818 +LOWZ_QSO: + EFFICIENCY: + SIGMOID_CUTOFF: 22.957808255020066 + SIGMOID_FUDGE: 0.12355429502045497 + FAILURE_RATE: 0.004928698166524282 + UNCERTAINTY: + POWER_LAW_INDEX: 2.36704293 + SIGMA_17: 0.00306028 +LRG: + EFFICIENCY: + SIGMOID_CUTOFF: 1501.931900023393 + SIGMOID_FUDGE: 58.85348636544879 + FAILURE_RATE: 0.0 + UNCERTAINTY: + POWER_LAW_INDEX: 0.846356643 + SIGMA_17: -0.00000126317571 +LYA_QSO: + EFFICIENCY: + SIGMOID_CUTOFF: 23.038943272889487 + SIGMOID_FUDGE: 0.09444368920722161 + FAILURE_RATE: 0.0019649874955341195 + UNCERTAINTY: + POWER_LAW_INDEX: 7.24851166 + SIGMA_17: 0.00161001792 +QSO_ZSPLIT: 2.1 diff --git a/py/desisim/data/quickcat.yaml b/py/desisim/data/quickcat.yaml index 2e537b0d1..a9231dead 100644 --- a/py/desisim/data/quickcat.yaml +++ b/py/desisim/data/quickcat.yaml @@ -31,4 +31,4 @@ LYA_QSO: UNCERTAINTY: POWER_LAW_INDEX: 0.17072327837058882 SIGMA_17: 0.0004556498419243293 -QSO_ZSPLIT: 2.0 +QSO_ZSPLIT: 2.0 \ No newline at end of file diff --git a/py/desisim/quickcat.py b/py/desisim/quickcat.py index a8550c670..ab468c063 100644 --- a/py/desisim/quickcat.py +++ b/py/desisim/quickcat.py @@ -16,7 +16,7 @@ import numpy as np from astropy.io import fits -from astropy.table import Table, Column, vstack +from astropy.table import Table, Column, vstack, join import sys import scipy.special as sp import desisim @@ -30,7 +30,7 @@ from desiutil.log import get_logger log = get_logger() -#- redshift errors, zwarn, cata fail rate fractions from +#- redshift errors, zwarn, cata. fail rate fractions from #- /project/projectdirs/desi/datachallenge/redwood/spectro/redux/redwood/ #- sigmav = c sigmaz / (1+z) _sigma_v = { @@ -38,6 +38,7 @@ # 'LRG': 67.38, 'BGS': 37.70, # 'QSO': 182.16, + 'MWS': 50.00, 'STAR': 51.51, 'WD':54.35, 'SKY': 9999, #- meaningless @@ -111,9 +112,10 @@ def get_zeff_obs(simtype, obsconditions): x = obsconditions['SEEING'] - np.mean(obsconditions['SEEING']) px = p_x[0] + p_x[1]*x + p_x[2] * (x**2 - np.mean(x**2)) - # transparency - if 'LINTRANS' in obsconditions : - y = obsconditions['LINTRANS'] - np.mean(obsconditions['LINTRANS']) + # transparency; + # https://desidatamodel.readthedocs.io/en/latest/DESISURVEY_OUTPUT/exposures.html + if 'TRANSP' in obsconditions : + y = obsconditions['TRANSP'] - np.mean(obsconditions['TRANSP']) py = p_y[0] + p_y[1]*y + p_y[2] * (y**2 - np.mean(y**2)) else : py = np.ones(ncond) @@ -148,7 +150,7 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi 'TILEID': array of tile IDs 'AIRMASS': array of airmass values on a tile 'EBMV': array of E(B-V) values on a tile - 'LINTRANS': array of atmospheric transparency during spectro obs; floats [0-1] + 'TRANSP': array of atmospheric transparency during spectro obs; floats [0-1] 'MOONFRAC': array of moonfraction values on a tile. 'SEEING': array of FWHM seeing during spectroscopic observation on a tile. parameter_filename: yaml file with quickcat parameters @@ -164,26 +166,19 @@ def get_redshift_efficiency(simtype, targets, truth, targets_in_tile, obsconditi n = len(targetid) try: - if 'DECAM_FLUX' in targets.dtype.names : - true_gflux = targets['DECAM_FLUX'][:, 1] - true_rflux = targets['DECAM_FLUX'][:, 2] - else: true_gflux = targets['FLUX_G'] true_rflux = targets['FLUX_R'] + true_zflux = targets['FLUX_Z'] except: raise Exception('Missing photometry needed to estimate redshift efficiency!') a_small_flux=1e-40 true_gflux[true_gflux 0.0), "{} {}".format(objtype, zerr[ii].min()) + zout[ii] += np.random.normal(scale=zerr[ii]) - + log.info("ELG sigma={:6.5f} index={:3.2f} median zerr={:6.5f}".format(sigma,powerlawindex,np.median(zerr[ii]))) elif objtype == 'LRG' : + sigma = np.float(params["LRG"]["UNCERTAINTY"]["SIGMA_17"]) + powerlawindex = params["LRG"]["UNCERTAINTY"]["POWER_LAW_INDEX"] + zerr[ii] = sigma/(1.e-9 + true_rflux[ii]**powerlawindex)*(1. + truez[ii]) - sigma = params["LRG"]["UNCERTAINTY"]["SIGMA_17"] - powerlawindex = params["LRG"]["UNCERTAINTY"]["POWER_LAW_INDEX"] + assert np.all(zerr[ii] > 0.0), "{} {} {} {} {}".format(objtype, zerr[ii].min(), sigma, true_rflux[ii].min(), truez[ii].min()) + + zout[ii] += np.random.normal(scale=zerr[ii]) - zerr[ii] = sigma/(1.e-9+true_rflux[ii]**powerlawindex)*(1.+truez[ii]) - zout[ii] += np.random.normal(scale=zerr[ii]) - log.info("LRG sigma={:6.5f} index={:3.2f} median zerr={:6.5f}".format(sigma,powerlawindex,np.median(zerr[ii]))) elif objtype == 'QSO' : - - zsplit = params['QSO_ZSPLIT'] + zsplit = params['QSO_ZSPLIT'] sigma = params["LOWZ_QSO"]["UNCERTAINTY"]["SIGMA_17"] powerlawindex = params["LOWZ_QSO"]["UNCERTAINTY"]["POWER_LAW_INDEX"] jj=ii&(truth['TRUEZ']<=zsplit) zerr[jj] = sigma/(1.e-9+(true_rflux[jj])**powerlawindex)*(1.+truez[jj]) + assert np.all(zerr[jj] > 0.0), "{} {}".format(objtype, zerr[ii].min()) + log.info("LOWZ QSO sigma={:6.5f} index={:3.2f} median zerr={:6.5f}".format(sigma,powerlawindex,np.median(zerr[jj]))) sigma = params["LYA_QSO"]["UNCERTAINTY"]["SIGMA_17"] powerlawindex = params["LYA_QSO"]["UNCERTAINTY"]["POWER_LAW_INDEX"] jj=ii&(truth['TRUEZ']>zsplit) zerr[jj] = sigma/(1.e-9+(true_rflux[jj])**powerlawindex)*(1.+truez[jj]) + + assert np.all(zerr[jj] > 0.0), "{} {}".format(objtype, zerr[ii].min()) log.info("LYA QSO sigma={:6.5f} index={:3.2f} median zerr={:6.5f}".format(sigma,powerlawindex,np.median(zerr[jj]))) @@ -444,7 +438,11 @@ def get_observed_redshifts(targets, truth, targets_in_tile, obsconditions, param log.info("{} use constant sigmav = {} km/s".format(objtype,_sigma_v[objtype])) ii = (simtype == objtype) zerr[ii] = _sigma_v[objtype] * (1+truez[ii]) / c + + assert np.all(zerr[ii] > 0.0), "{} {}".format(objtype, zerr[ii].min()) + zout[ii] += np.random.normal(scale=zerr[ii]) + else : log.info("{} no redshift error model, will use truth") @@ -507,7 +505,7 @@ def get_median_obsconditions(tileids): 'TILEID': array of tile IDs 'AIRMASS': array of airmass values on a tile 'EBMV': array of E(B-V) values on a tile - 'LINTRANS': array of atmospheric transparency during spectro obs; floats [0-1] + 'TRANSP': array of atmospheric transparency during spectro obs; floats [0-1] 'MOONFRAC': array of moonfraction values on a tile. 'SEEING': array of FWHM seeing during spectroscopic observation on a tile. """ @@ -539,7 +537,7 @@ def get_median_obsconditions(tileids): obsconditions['TILEID'] = tileids obsconditions['AIRMASS'] = tiles['AIRMASS'] obsconditions['EBMV'] = tiles['EBV_MED'] - obsconditions['LINTRANS'] = np.ones(n) + obsconditions['TRANSP'] = np.ones(n) obsconditions['SEEING'] = np.ones(n) * 1.1 #- Add lunar conditions, defaulting to dark time @@ -560,6 +558,70 @@ def get_median_obsconditions(tileids): return obsconditions +def exp_derivedprops(exposures, tiles=None): + ''' + Given RA, DEC, and MJD from exposures table, return derived proprs for Kitt Peak, e.g. MOONFRAC, MOONSEP etc. + Note: temporary copy of https://github.com/changhoonhahn/feasiBGS/blob/master/run/bright_exposure/surveysim_output.py, line 384. + ''' + import ephem + import desisurvey.config + import desisurvey.utils as dutils + import astropy.units as u + + from astropy.time import Time + + + config = desisurvey.config.Configuration() + + mayall = ephem.Observer() + mayall.lat = config.location.latitude().to(u.rad).value + mayall.lon = config.location.longitude().to(u.rad).value + mayall.elevation = config.location.elevation().to(u.m).value + + ## Observed time (MJD). + mjd = Time(exposures['MJD'].quantity.value, format='mjd') + + exposures['MOONALT'] = np.zeros(len(mjd)) + exposures['MOONRA'] = np.zeros(len(mjd)) + exposures['MOONDEC'] = np.zeros(len(mjd)) + exposures['MOONFRAC'] = np.zeros(len(mjd)) + + exposures['SUNALT'] = np.zeros(len(mjd)) + exposures['SUNRA'] = np.zeros(len(mjd)) + exposures['SUNDEC'] = np.zeros(len(mjd)) + + for i in range(len(mjd)): + mayall.date = mjd.datetime[i] + + _moon = ephem.Moon() + _moon.compute(mayall) + + _sun = ephem.Sun() + _sun.compute(mayall) + + exposures['MOONALT'][i] = 180./np.pi*_moon.alt + exposures['MOONRA'][i] = 180./np.pi*_moon.ra + exposures['MOONDEC'][i] = 180./np.pi*_moon.dec + exposures['MOONFRAC'][i] = _moon.moon_phase + + exposures['SUNALT'][i] = 180./np.pi*_sun.alt + exposures['SUNRA'][i] = 180./np.pi*_sun.ra + exposures['SUNDEC'][i] = 180./np.pi*_sun.dec + + if tiles is not None: + exposures = join(exposures, tiles, keys=['TILEID'], join_type='outer', table_names=['', '~'], uniq_col_name='{col_name}{table_name}') + exposures.remove_column('AIRMASS~') + + exposures['MOONSEP'] = np.diag(dutils.separation_matrix(exposures['MOONRA'].quantity * u.deg, exposures['MOONDEC'].quantity * u.deg,\ + np.atleast_1d(exposures['RA'].quantity * u.deg), np.atleast_1d(exposures['DEC'].quantity * u.deg))) + + exposures['SUNSEP'] = np.diag(dutils.separation_matrix(exposures['SUNRA'].quantity * u.deg, exposures['SUNDEC'].quantity * u.deg,\ + np.atleast_1d(exposures['RA'].quantity * u.deg), np.atleast_1d(exposures['DEC'].quantity * u.deg))) + + exposures.remove_rows(np.where(exposures['EXPTIME'].mask == True)) + + return exposures + def quickcat(tilefiles, targets, truth, zcat=None, obsconditions=None, perfect=False): """ Generates quick output zcatalog @@ -588,7 +650,7 @@ def quickcat(tilefiles, targets, truth, zcat=None, obsconditions=None, perfect=F tileids = list() for infile in tilefiles: - fibassign, header = fits.getdata(infile, 'FIBERASSIGN', header=True) + fibassign, header = fits.getdata(infile, 'FASSIGN', header=True) # hack needed here rnc 7/26/18 if 'TILEID' in header: @@ -607,7 +669,7 @@ def quickcat(tilefiles, targets, truth, zcat=None, obsconditions=None, perfect=F targets_in_tile[tileidnew] = fibassign['TARGETID'][ii] #- Trim obsconditions to just the tiles that were observed - if obsconditions is not None: + if obsconditions is not None and len(tileids) > 0: ii = np.in1d(obsconditions['TILEID'], tileids) if np.any(ii == False): obsconditions = obsconditions[ii] @@ -617,12 +679,25 @@ def quickcat(tilefiles, targets, truth, zcat=None, obsconditions=None, perfect=F #- This might not be needed, but is fast for O(20k) tiles and may #- prevent future surprises if code expects them to be row aligned tileids = np.array(tileids) + + if (obsconditions is not None): + ## Handle repeated TILEIDs for Cosmic splits. + order = [list(np.where(obsconditions['TILEID'].quantity == x)[0]) for x in tileids] + order = [item for sublist in order for item in sublist] + + obsconditions = obsconditions[order] + + + if (obsconditions is not None) and \ (np.any(tileids != obsconditions['TILEID'])): - i = np.argsort(tileids) - j = np.argsort(obsconditions['TILEID']) - k = np.argsort(i) - obsconditions = obsconditions[j[k]] +# i = np.argsort(tileids) +# j = np.argsort(obsconditions['TILEID']) +# k = np.argsort(i) +# obsconditions = obsconditions[j[k]] + match_tiles = np.intersect1d(obsconditions['TILEID'],tileids,return_indices=True)[1] + order_tiles = np.argsort(tileids) + obsconditions = obsconditions[match_tiles][order_tiles] assert np.all(tileids == obsconditions['TILEID']) #- Trim truth down to just ones that have already been observed @@ -644,6 +719,13 @@ def quickcat(tilefiles, targets, truth, zcat=None, obsconditions=None, perfect=F #- Copy TRUETYPE -> SPECTYPE so that we can change without altering original newzcat['SPECTYPE'] = truth['TRUESPECTYPE'].copy() + #- Propagate target/truth information to output catalog + newzcat['RA'] = targets['RA'].copy() + newzcat['DEC'] = targets['DEC'].copy() + newzcat['TRUEZ'] = truth['TRUEZ'].copy() + newzcat['FLUX_R'] = truth['FLUX_R'].copy() + newzcat['OIIFLUX'] = truth['OIIFLUX'].copy() + #- Add ZERR and ZWARN log.info('{} QC Adding ZERR and ZWARN'.format(asctime())) nz = len(newzcat) diff --git a/py/desisim/quicksurvey.py b/py/desisim/quicksurvey.py index 709e6d4b1..0a6ee9478 100644 --- a/py/desisim/quicksurvey.py +++ b/py/desisim/quicksurvey.py @@ -19,13 +19,17 @@ import shutil import glob import subprocess +import desimodel +import astropy.io.fits as fits from astropy.table import Table, Column +from astropy.io import fits import os.path from collections import Counter from time import time, asctime +from astropy.time import Time import fitsio import desitarget.mtl -from desisim.quickcat import quickcat +from desisim.quickcat import quickcat, exp_derivedprops from astropy.table import join from desitarget.targetmask import desi_mask @@ -41,7 +45,7 @@ class SimSetup(object): n_epochs (int): number of epochs to be simulated. """ - def __init__(self, output_path, targets_path, fiberassign, exposures, fiberassign_dates): + def __init__(self, output_path, targets_path, fiberassign, exposures, fiberassign_dates, footprint=None, status=None): """Initializes all the paths, filenames and numbers describing DESI survey. Args: @@ -59,12 +63,13 @@ def __init__(self, output_path, targets_path, fiberassign, exposures, fiberassig self.tmp_output_path = os.path.join(self.output_path, 'tmp/') self.tmp_fiber_path = os.path.join(self.tmp_output_path, 'fiberassign/') + self.footprint = footprint self.surveyfile = os.path.join(self.tmp_output_path, 'survey_list.txt') self.skyfile = os.path.join(self.targets_path,'sky.fits') - self.stdfile = os.path.join(self.targets_path,'std.fits') + self.stdfile = os.path.join(self.targets_path,'standards.fits') self.truthfile = os.path.join(self.targets_path,'truth.fits') self.targetsfile = os.path.join(self.targets_path,'targets.fits') - self.fibstatusfile = os.path.join(self.targets_path,'fiberstatus.ecsv') + self.status = status self.zcat_file = None self.mtl_file = None @@ -75,8 +80,11 @@ def __init__(self, output_path, targets_path, fiberassign, exposures, fiberassig self.epochs_list = list() self.n_epochs = 0 self.start_epoch = 0 + + if self.footprint == None: + ## Default to the nominal DESI footprint a la fiberassign. + self.footprint = desimodel.io.findfile('footprint/desi-tiles.fits') - dateobs = np.core.defchararray.decode(self.exposures['NIGHT']) dates = list() with open(fiberassign_dates) as fx: for line in fx: @@ -85,26 +93,29 @@ def __init__(self, output_path, targets_path, fiberassign, exposures, fiberassig continue yearmmdd = line.replace('-', '') year_mm_dd = yearmmdd[0:4]+yearmmdd[4:6]+yearmmdd[6:8] - dates.append(year_mm_dd) - - #- add pre- and post- dates for date range bookkeeping - if dates[0] < min(dateobs[0]): - dates.insert(0, dateobs[0]) - - dates.append('9999-99-99') - print(dates) + dates.append(np.int(year_mm_dd)) self.n_epochs = len(dates) - 1 + ## Convert MJD present in exposures to 'NIGHT' in yyyymmdd format. + iso = Time(self.exposures['MJD'], format='mjd', scale='utc').iso + dateobs = np.array([x.split(' ')[0].replace('-', '') for x in iso]).astype(np.int) + + #- add pre- and post- dates for date range bookkeeping + if dates[0] < min(dateobs): + dates.insert(19000131, dateobs[0]) + dates.append(30000131) + + print(dates) for i in range(len(dates)-1): - ii = (dateobs >= dates[i]) & (dateobs < dates[i+1]) - epoch_tiles = list() + ii = (dateobs >= dates[i]) & (dateobs < dates[i+1]) + epoch_tiles = list(np.unique(self.exposures['TILEID'][ii])) for tile in self.exposures['TILEID'][ii]: if tile not in epoch_tiles: epoch_tiles.append(tile) self.epoch_tiles.append(epoch_tiles) - print('tiles in epoch {} [{} to {}]: {}'.format(i,dates[i], dates[i+1], len(self.epoch_tiles[i]))) + print('tiles in epoch {} [{} to {}]: {}'.format(i, dates[i], dates[i+1], len(epoch_tiles))) def create_directories(self): @@ -189,19 +200,15 @@ def update_observed_tiles(self, epoch): self.tilefiles = list() tiles = self.epoch_tiles[epoch] for i in tiles: - tilename = os.path.join(self.tmp_fiber_path, 'tile-%06d.fits'%(i)) - # retain ability to use previous version of tile files - oldtilename = os.path.join(self.tmp_fiber_path, 'tile_%05d.fits'%(i)) + tilename = os.path.join(self.tmp_fiber_path, 'fba-{:06d}.fits'.format(i)) if os.path.isfile(tilename): self.tilefiles.append(tilename) - elif os.path.isfile(oldtilename): - self.tilefiles.append(oldtilename) - #else: - # print('Suggested but does not exist {}'.format(tilename)) + else: + print('Suggested but does not exist {}'.format(tilename)) print("{} {} tiles to gather in zcat".format(asctime(), len(self.tilefiles))) - def simulate_epoch(self, epoch, truth, targets, perfect=False, zcat=None): + def simulate_epoch(self, epoch, truth, targets, obscon, perfect=False, zcat=None): """Core routine simulating a DESI epoch, Args: @@ -211,6 +218,7 @@ def simulate_epoch(self, epoch, truth, targets, perfect=False, zcat=None): False: redshifts include uncertainties. truth (Table): Truth data targets (Table): Targets data + obscon (str): A combination of strings that are in the desitarget bitmask yaml file, e.g. "DARK|GRAY". Governs behavior of how priorities are set. zcat (Table): Redshift Catalog Data Notes: This routine simulates three steps: @@ -223,9 +231,9 @@ def simulate_epoch(self, epoch, truth, targets, perfect=False, zcat=None): print("{} Starting MTL".format(asctime())) self.mtl_file = os.path.join(self.tmp_output_path, 'mtl.fits') if zcat is None: - mtl = desitarget.mtl.make_mtl(targets) + mtl = desitarget.mtl.make_mtl(targets, obscon) else: - mtl = desitarget.mtl.make_mtl(targets, zcat) + mtl = desitarget.mtl.make_mtl(targets, obscon, zcat) mtl.write(self.mtl_file, overwrite=True) del mtl @@ -233,7 +241,7 @@ def simulate_epoch(self, epoch, truth, targets, perfect=False, zcat=None): print("{} Finished MTL".format(asctime())) # clean files and prepare fiberasign inputs - tilefiles = sorted(glob.glob(self.tmp_fiber_path+'/tile*.fits')) + tilefiles = sorted(glob.glob(self.tmp_fiber_path+'/fba*.fits')) if tilefiles: for tilefile in tilefiles: os.remove(tilefile) @@ -244,18 +252,19 @@ def simulate_epoch(self, epoch, truth, targets, perfect=False, zcat=None): # launch fiberassign print("{} Launching fiberassign".format(asctime())) - f = open('fiberassign.log','a') - - p = subprocess.call([self.fiberassign, - '--mtl', os.path.join(self.tmp_output_path, 'mtl.fits'), - '--stdstar', self.stdfile, - '--sky', self.skyfile, - '--surveytiles', self.surveyfile, - '--outdir',os.path.join(self.tmp_output_path, 'fiberassign'), - '--fibstatusfile', self.fibstatusfile], - stdout=f) + f = open('fiberassign.log','a') + + # Assume stds within mtl. {'--stdstar', self.stdfile} otherwise. + call = [self.fiberassign, '--mtl', os.path.join(self.tmp_output_path, 'mtl.fits'), '--sky', self.skyfile,\ + '--surveytiles', self.surveyfile, '--footprint', self.footprint, '--outdir', os.path.join(self.tmp_output_path, 'fiberassign'), '--overwrite'] + if self.status is not None: + call.append('--status') + call.append(self.status) + ## Doesn't catch no overwrite return on fiberassign. + p = subprocess.call(call, stdout=f) + print("{} Finished fiberassign".format(asctime())) f.close() @@ -266,7 +275,33 @@ def simulate_epoch(self, epoch, truth, targets, perfect=False, zcat=None): # progress_data = Table.read(self.progress_files[epoch + 1]) # ii = np.in1d(progress_data['TILEID'], self.observed_tiles) # obsconditions = progress_data[ii] - obsconditions = None + + #- Use median conditions per tile or get observing info from exposures file + median_conditions = False + if median_conditions: + obsconditions = None + else: + # update obsconditions, must contain: TILEID, AIRMASS, EBMV, TRANS, MOONFRAC, SEEING. + exposures = Table(self.exposures) + tiles = Table(fits.open(self.footprint)[1].data) + + ii = np.in1d(exposures['TILEID'].quantity, self.epoch_tiles[epoch]) + + obsconditions = exp_derivedprops(exposures[ii], tiles) + + print('\n\nObs conditions for epoch: {}.'.format(epoch)) + print(obsconditions) + + print('Solving for {} tilefiles.'.format(len(self.tilefiles))) + + obsconditions = Table() + obsconditions['TILEID'] = self.exposures['TILEID'] + obsconditions['AIRMASS'] = self.exposures['AIRMASS'] + obsconditions['SEEING'] = self.exposures['SEEING'] + obsconditions['TRANSP'] = self.exposures['TRANSP'] + obsconditions['MOONFRAC'] = self.exposures['MOONFRAC'] + obsconditions['MOONALT'] = self.exposures['MOONALT'] + print('tilefiles', len(self.tilefiles)) # write the zcat, it uses the tilesfiles constructed in the last step @@ -287,8 +322,18 @@ def simulate(self): """ self.create_directories() - truth = Table.read(os.path.join(self.targets_path,'truth.fits')) - targets = Table.read(os.path.join(self.targets_path,'targets.fits')) + truth = Table.read(self.truthfile,hdu='TRUTH') + targets = Table.read(self.targetsfile) + obscon = fits.open(self.truthfile)['TRUTH'].header['OBSCON'] + + #- Add OII flux from ELG HDU to truth table + truth_elg = Table.read(self.truthfile,hdu='TRUTH_ELG') + truthid = truth['TARGETID'] + elgid = truth_elg['TARGETID'] + match_elg = np.intersect1d(truthid,elgid,return_indices=True)[1] + oiiflux = np.zeros(truthid.shape) + oiiflux[match_elg] = truth_elg['OIIFLUX'] + truth['OIIFLUX'] = oiiflux print(truth.keys()) #- Drop columns that aren't needed to save memory while manipulating @@ -317,7 +362,7 @@ def simulate(self): zcat = Table.read(os.path.join(epochdir, 'zcat.fits')) # Update mtl and zcat - self.simulate_epoch(epoch, truth, targets, perfect=True, zcat=zcat) + self.simulate_epoch(epoch, truth, targets, obscon, perfect=True, zcat=zcat) # copy mtl and zcat to epoch directory self.backup_epoch_data(epoch_id=epoch) diff --git a/py/desisim/test/test_quickcat.py b/py/desisim/test/test_quickcat.py index ff01ea62c..691a4e8a2 100644 --- a/py/desisim/test/test_quickcat.py +++ b/py/desisim/test/test_quickcat.py @@ -1,10 +1,10 @@ import os import numpy as np import unittest -from astropy.table import Table, Column -from astropy.io import fits -from desisim.quickcat import quickcat -from desitarget.targetmask import desi_mask, bgs_mask, mws_mask +from astropy.table import Table, Column +from astropy.io import fits +from desisim.quickcat import quickcat +from desitarget.targetmask import desi_mask, bgs_mask, mws_mask import desimodel.io class TestQuickCat(unittest.TestCase): @@ -20,10 +20,12 @@ def setUpClass(cls): cls.nspec = n = 5000 targets = Table() - targets['TARGETID'] = np.random.randint(0,2**60, size=n) + targets['TARGETID'] = np.random.randint(0,2**60, size=n) + targets['RA'] = np.random.uniform(0,360., size=n) + targets['DEC'] = np.random.uniform(-90.,90., size=n) targets['DESI_TARGET'] = 2**np.random.randint(0,3,size=n) - targets['BGS_TARGET'] = np.zeros(n, dtype=int) - targets['MWS_TARGET'] = np.zeros(n, dtype=int) + targets['BGS_TARGET'] = np.zeros(n, dtype=int) + targets['MWS_TARGET'] = np.zeros(n, dtype=int) isLRG = (targets['DESI_TARGET'] & desi_mask.LRG) != 0 isELG = (targets['DESI_TARGET'] & desi_mask.ELG) != 0 isQSO = (targets['DESI_TARGET'] & desi_mask.QSO) != 0 @@ -44,16 +46,22 @@ def setUpClass(cls): targets['MWS_TARGET'][isMWS] = mws_mask.MWS_MAIN #- Add some fake photometry; no attempt to get colors right + #- https://portal.nersc.gov/project/desi/users/adamyers/0.31.1/desitargetQA-dr8-0.31.1/LRG.html flux = np.zeros((n, 6)) #- ugrizY; DESI has grz - flux[isLRG, 1] = np.random.uniform(0, 1.0, np.count_nonzero(isLRG)) - flux[isLRG, 2] = np.random.uniform(0, 5.0, np.count_nonzero(isLRG)) - flux[isLRG, 4] = np.random.uniform(0, 5.0, np.count_nonzero(isLRG)) + + #- Assumed nanomaggies. + flux[isLRG, 1] = np.random.uniform(19., 20., np.count_nonzero(isLRG)) + flux[isLRG, 2] = np.random.uniform(19., 20., np.count_nonzero(isLRG)) + flux[isLRG, 4] = np.random.uniform(19., 20., np.count_nonzero(isLRG)) + flux[isELG, 1] = np.random.uniform(0, 4.0, np.count_nonzero(isELG)) flux[isELG, 2] = np.random.uniform(0, 4.0, np.count_nonzero(isELG)) flux[isELG, 4] = np.random.uniform(0, 10.0, np.count_nonzero(isELG)) + flux[isQSO, 1] = np.random.uniform(0, 4.0, np.count_nonzero(isQSO)) flux[isQSO, 2] = np.random.uniform(0, 4.0, np.count_nonzero(isQSO)) flux[isQSO, 4] = np.random.uniform(0, 6.0, np.count_nonzero(isQSO)) + # isBGS and isMWS are arrays of indices, not arrays of booleans flux[isBGS, 1] = np.random.uniform(10, 600, isBGS.size) flux[isBGS, 2] = np.random.uniform(15, 1000, isBGS.size) @@ -61,12 +69,28 @@ def setUpClass(cls): flux[isMWS, 1] = np.random.uniform(10, 150, isMWS.size) flux[isMWS, 2] = np.random.uniform(15, 350, isMWS.size) flux[isMWS, 4] = np.random.uniform(10, 1500, isMWS.size) - targets['DECAM_FLUX'] = flux + targets['FLUX_G'] = flux[:,1] + targets['FLUX_R'] = flux[:,2] + targets['FLUX_Z'] = flux[:,4] + truth = Table() truth['TARGETID'] = targets['TARGETID'].copy() + + truth['FLUX_G'] = targets['FLUX_G'].copy() + truth['FLUX_R'] = targets['FLUX_R'].copy() + truth['FLUX_Z'] = targets['FLUX_Z'].copy() + truth['TRUEZ'] = np.random.uniform(0, 1.5, size=n) truth['TRUESPECTYPE'] = np.zeros(n, dtype=(str, 10)) + + truth['TEMPLATETYPE'] = np.zeros(n, dtype=(str, 10)) + truth['TEMPLATETYPE'][isLRG] = 'LRG' + truth['TEMPLATETYPE'][isELG] = 'ELG' + truth['TEMPLATETYPE'][isQSO] = 'QSO' + truth['TEMPLATETYPE'][isBGS] = 'BGS' + truth['TEMPLATETYPE'][isMWS] = 'MWS' + truth['GMAG'] = np.random.uniform(18.0, 24.0, size=n) ii = (targets['DESI_TARGET'] & desi_mask.mask('LRG|ELG|BGS_ANY')) != 0 truth['TRUESPECTYPE'][ii] = 'GALAXY' @@ -87,7 +111,7 @@ def setUpClass(cls): fiberassign = truth['TARGETID',] fiberassign['RA'] = np.random.uniform(0,5, size=n) fiberassign['DEC'] = np.random.uniform(0,5, size=n) - fiberassign.meta['EXTNAME'] = 'FIBERASSIGN' + fiberassign.meta['EXTNAME'] = 'FASSIGN' nx = cls.nspec // cls.ntiles cls.targets_in_tile = dict() for i, filename in enumerate(cls.tilefiles): @@ -132,7 +156,7 @@ def test_quickcat(self): zcat2 = quickcat(self.tilefiles[0:2], self.targets, truth=self.truth, perfect=False) zcat2_sorted = zcat2.copy() zcat2_sorted.sort(keys='TARGETID') - self.assertTrue(np.all(zcat2_sorted['TARGETID'] == truth_01['TARGETID'])) + self.assertTrue(np.all(zcat2_sorted['TARGETID'] == truth_01['TARGETID'])) self.assertTrue(np.all(zcat2_sorted['Z'] != truth_01['TRUEZ'])) self.assertTrue(np.any(zcat2_sorted['ZWARN'] != 0))