diff --git a/bin/build-fsps-templates b/bin/build-fsps-templates old mode 100644 new mode 100755 diff --git a/bin/fastqa b/bin/fastqa new file mode 100755 index 00000000..22b114df --- /dev/null +++ b/bin/fastqa @@ -0,0 +1,9 @@ +#!/usr/bin/env python +"""fastspecfit QA + +fastqa /global/cfs/cdirs/desi/spectro/fastspecfit/blanc/tiles/80606/deep/fastphot-0-80606-deep.fits --outdir . --ntargets 1 --firsttarget 50 + +""" +if __name__ == '__main__': + from fastspecfit.qa import fastqa + fastqa(comm=None) diff --git a/bin/fastspecfit-qa b/bin/fastspecfit-qa index f09836de..b07da155 100755 --- a/bin/fastspecfit-qa +++ b/bin/fastspecfit-qa @@ -1,279 +1,6 @@ #!/usr/bin/env python -"""fastspecfit QA - -fastspecfit-qa --fastphotfile /global/cfs/cdirs/desi/spectro/fastspecfit/blanc/tiles/80606/deep/fastphot-0-80606-deep.fits --outdir . --ntargets 20 --firsttarget 50 - -""" -import pdb # for debugging -import os, sys, time -import numpy as np - -from desiutil.log import get_logger -log = get_logger() - -def parse(options=None): - """Parse input arguments. - - """ - import sys, argparse - - parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) - - parser.add_argument('--healpix', default=None, type=str, nargs='*', help="""Generate QA for all objects - with this healpixels (only defined for coadd-type 'healpix').""") - parser.add_argument('--tile', default=None, type=str, nargs='*', help='Generate QA for all objects on this tile.') - parser.add_argument('--night', default=None, type=str, nargs='*', help="""Generate QA for all objects observed on this - night (only defined for coadd-type 'pernight' and 'perexp').""") - parser.add_argument('--redux_dir', default=None, type=str, help='Top-level path to the reduced spectra.') - parser.add_argument('--redrockfiles', nargs='*', help='Optional full path to redrock file(s).') - parser.add_argument('--redrockfile-prefix', type=str, default='redrock-', help='Prefix of the input Redrock file name(s).') - parser.add_argument('--specfile-prefix', type=str, default='coadd-', help='Prefix of the spectral file(s).') - parser.add_argument('--qnfile-prefix', type=str, default='qso_qn-', help='Prefix of the QuasarNet afterburner file(s).') - parser.add_argument('--mapdir', type=str, default=None, help='Optional directory name for the dust maps.') - parser.add_argument('--fphotodir', type=str, default=None, help='Top-level location of the source photometry.') - parser.add_argument('--fphotoinfo', type=str, default=None, help='Photometric information file.') - - parser.add_argument('--emline_snrmin', type=float, default=0.0, help='Minimum emission-line S/N to be displayed.') - parser.add_argument('--nsmoothspec', type=int, default=0, help='Smoothing pixel value.') - - parser.add_argument('--minspecwave', type=float, default=3500., help='Minimum spectral wavelength (Angstrom).') - parser.add_argument('--maxspecwave', type=float, default=9900., help='Maximum spectral wavelength (Angstrom).') - parser.add_argument('--minphotwave', type=float, default=0.1, help='Minimum photometric wavelength (micron).') - parser.add_argument('--maxphotwave', type=float, default=35., help='Maximum photometric wavelength (micron).') - - parser.add_argument('--targetids', type=str, default=None, help='Comma-separated list of target IDs to process.') - parser.add_argument('-n', '--ntargets', type=int, help='Number of targets to process in each file.') - parser.add_argument('--firsttarget', type=int, default=0, help='Index of first object to to process in each file (0-indexed).') - parser.add_argument('--mp', type=int, default=1, help='Number of multiprocessing processes per MPI rank or node.') - parser.add_argument('--nophoto', action='store_true', help='Do not include the photometry in the model fitting.') - parser.add_argument('--overwrite', action='store_true', help='Overwrite existing files.') - - parser.add_argument('--imf', type=str, default='chabrier', help='Initial mass function.') - parser.add_argument('--templateversion', type=str, default='1.0.0', help='Template version number.') - parser.add_argument('--templates', type=str, default=None, help='Optional full path and filename to the templates.') - - parser.add_argument('--outprefix', default=None, type=str, help='Optional prefix for output filename.') - parser.add_argument('-o', '--outdir', default='.', type=str, help='Full path to desired output directory.') - - parser.add_argument('fastfitfile', nargs=1, help='Full path to fastspec or fastphot fitting results.') - - if options is None: - args = parser.parse_args() - log.info(' '.join(sys.argv)) - else: - args = parser.parse_args(options) - log.info('fastspecfit-qa {}'.format(' '.join(options))) - - return args - -def main(args=None, comm=None): - """Main module. - - """ - import fitsio - from astropy.table import Table - from fastspecfit.fastspecfit import _desiqa_one, desiqa_one - from fastspecfit.io import (DESISpectra, get_templates_filename, get_qa_filename, - read_fastspecfit, DESI_ROOT_NERSC, select) - - if isinstance(args, (list, tuple, type(None))): - args = parse(args) - - if args.redux_dir is None: - args.redux_dir = os.path.join(os.environ.get('DESI_ROOT', DESI_ROOT_NERSC), 'spectro', 'redux') - if not os.path.isdir(args.redux_dir): - errmsg = 'Data reduction directory {} not found.'.format(args.redux_dir) - log.critical(errmsg) - raise IOError(errmsg) - - # Read the fitting results. - if not os.path.isfile(args.fastfitfile[0]): - log.warning('File {} not found.'.format(args.fastfitfile[0])) - return - - fastfit, metadata, coadd_type, fastphot = read_fastspecfit(args.fastfitfile[0]) - - # parse the targetids optional input - if args.targetids: - targetids = [int(x) for x in args.targetids.split(',')] - keep = np.where(np.isin(fastfit['TARGETID'], targetids))[0] - if len(keep) == 0: - log.warning('No matching targetids found!') - return - fastfit = fastfit[keep] - metadata = metadata[keep] - - if args.ntargets is not None: - keep = np.arange(args.ntargets) + args.firsttarget - log.info('Keeping {} targets.'.format(args.ntargets)) - fastfit = fastfit[keep] - metadata = metadata[keep] - - fastfit, metadata = select(fastfit, metadata, coadd_type, healpixels=args.healpix, - tiles=args.tile, nights=args.night) - - if coadd_type == 'custom' and args.redrockfiles is None: - errmsg = 'redrockfiles input is required if coadd_type==custom' - log.critical(errmsg) - raise IOError(errmsg) - - # check for the input redshifts header card - hdr = fitsio.read_header(args.fastfitfile[0]) - if 'INPUTZ' in hdr and hdr['INPUTZ']: - inputz = True - else: - inputz = False - - if 'NOSCORR' in hdr and hdr['NOSCORR']: - no_smooth_continuum = True - else: - no_smooth_continuum = False - - if args.outdir: - if not os.path.isdir(args.outdir): - os.makedirs(args.outdir, exist_ok=True) - - pngfile = get_qa_filename(metadata, coadd_type, outprefix=args.outprefix, - outdir=args.outdir, fastphot=fastphot, log=log) - - # are we overwriting? - if args.overwrite is False: - I = np.array([not os.path.isfile(_pngfile) for _pngfile in pngfile]) - J = ~I - if np.sum(J) > 0: - log.info('Skipping {} existing QA files.'.format(np.sum(J))) - fastfit = fastfit[I] - metadata = metadata[I] - - if len(metadata) == 0: - log.info('Done making all QA files!') - return - - log.info('Building QA for {} objects.'.format(len(metadata))) - - # Initialize the I/O class. - - Spec = DESISpectra(redux_dir=args.redux_dir, fphotodir=args.fphotodir, - fphotoinfo=args.fphotoinfo, mapdir=args.mapdir) - - templates = get_templates_filename(templateversion=args.templateversion, imf=args.imf) - - def _wrap_qa(redrockfile, indx=None, stackfit=False): - if indx is None: - indx = np.arange(len(fastfit)) - - if stackfit: - stackids = fastfit['STACKID'][indx] - data = Spec.read_stacked(redrockfile, stackids=stackids, mp=args.mp) - args.nophoto = True # force nophoto=True - - minspecwave = np.min(data[0]['coadd_wave']) - 20 - maxspecwave = np.max(data[0]['coadd_wave']) + 20 - else: - targetids = fastfit['TARGETID'][indx] - if inputz: - input_redshifts = fastfit['Z'][indx] - else: - input_redshifts = None - Spec.select(redrockfiles=redrockfile, targetids=targetids, - input_redshifts=input_redshifts, - redrockfile_prefix=args.redrockfile_prefix, - specfile_prefix=args.specfile_prefix, - qnfile_prefix=args.qnfile_prefix) - data = Spec.read_and_unpack(fastphot=fastphot, synthphot=True, mp=args.mp) - - minspecwave = args.minspecwave - maxspecwave = args.maxspecwave - - qaargs = [(data[igal], fastfit[indx[igal]], metadata[indx[igal]], - templates, coadd_type, Spec.fphoto, minspecwave, maxspecwave, - args.minphotwave, args.maxphotwave, args.emline_snrmin, - args.nsmoothspec, fastphot, args.nophoto, stackfit, inputz, - no_smooth_continuum, args.outdir, args.outprefix, log) - for igal in np.arange(len(indx))] - if args.mp > 1: - import multiprocessing - with multiprocessing.Pool(args.mp) as P: - P.map(_desiqa_one, qaargs) - else: - [desiqa_one(*_qaargs) for _qaargs in qaargs] - - t0 = time.time() - if coadd_type == 'healpix': - if args.redrockfiles is not None: - for redrockfile in args.redrockfiles: - _wrap_qa(redrockfile) - else: - allspecprods = metadata['SPECPROD'].data - allsurveys = metadata['SURVEY'].data - allprograms = metadata['PROGRAM'].data - allpixels = metadata['HEALPIX'].data - for specprod in set(allspecprods): - for survey in set(allsurveys): - for program in set(allprograms): - for pixel in set(allpixels): - indx = np.where((specprod == allspecprods) * (survey == allsurveys) * - (program == allprograms) * (pixel == allpixels))[0] - if len(indx) == 0: - #log.warning('No object found with specprod={}, survey={}, program={}, and healpixel={}!'.format( - # specprod, survey, program, pixel)) - continue - redrockfile = os.path.join(args.redux_dir, specprod, 'healpix', str(survey), str(program), str(pixel // 100), - str(pixel), 'redrock-{}-{}-{}.fits'.format(survey, program, pixel)) - _wrap_qa(redrockfile, indx) - elif coadd_type == 'custom': - for redrockfile in args.redrockfiles: - _wrap_qa(redrockfile) - elif coadd_type == 'stacked': - for redrockfile in args.redrockfiles: - _wrap_qa(redrockfile, stackfit=True) - else: - allspecprods = metadata['SPECPROD'].data - alltiles = metadata['TILEID'].astype(str).data - allnights = metadata['NIGHT'].astype(str).data - allpetals = metadata['FIBER'].data // 500 - if coadd_type == 'cumulative': - for specprod in set(allspecprods): - for tile in set(alltiles): - for petal in set(allpetals): - indx = np.where((specprod == allspecprods) * (tile == alltiles) * (petal == allpetals))[0] - if len(indx) == 0: - #log.warning('No object found with tileid={} and petal={}!'.format( - # tile, petal)) - continue - redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'cumulative', str(tile), allnights[indx[0]], - 'redrock-{}-{}-thru{}.fits'.format(petal, tile, allnights[indx[0]])) - _wrap_qa(redrockfile, indx) - elif coadd_type == 'pernight': - for specprod in set(allspecprods): - for night in set(allnights): - for tile in set(alltiles): - for petal in set(allpetals): - indx = np.where((specprod == allspecprods) * (night == allnights) * - (tile == alltiles) * (petal == allpetals))[0] - if len(indx) == 0: - continue - redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'pernight', str(tile), str(night), - 'redrock-{}-{}-{}.fits'.format(petal, tile, night)) - _wrap_qa(redrockfile, indx) - elif coadd_type == 'perexp': - allexpids = metadata['EXPID'].data - for specprod in set(allspecprods): - for night in set(allnights): - for expid in set(allexpids): - for tile in set(alltiles): - for petal in set(allpetals): - indx = np.where((specprod == allspecprods) * (night == allnights) * - (expid == allexpids) * (tile == alltiles) * - (petal == allpetals))[0] - if len(indx) == 0: - continue - redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'perexp', str(tile), '{:08d}'.format(expid), - 'redrock-{}-{}-exp{:08d}.fits'.format(petal, tile, expid)) - _wrap_qa(redrockfile, indx) - - log.info('QA for everything took: {:.2f} sec'.format(time.time()-t0)) if __name__ == '__main__': - main() - + from desiutil.log import get_logger + log = get_logger() + log.warning('fastspecfit-qa is deprecated; please use fastqa') diff --git a/bin/fastspecfit-setup.sh b/bin/fastspecfit-setup.sh index fac96b39..3dc747f0 100755 --- a/bin/fastspecfit-setup.sh +++ b/bin/fastspecfit-setup.sh @@ -32,7 +32,7 @@ elif [[ $1 == "env" ]]; then export DESI_ROOT='/global/cfs/cdirs/desi' export DUST_DIR='/global/cfs/cdirs/cosmo/data/dust/v0_1' - export DR9_DIR='/global/cfs/cdirs/desi/external/legacysurvey/dr9' + export FPHOTO_DIR='/global/cfs/cdirs/desi/external/legacysurvey/dr9' export FTEMPLATES_DIR='/global/cfs/cdirs/desi/external/templates/fastspecfit' export OMP_NUM_THREADS=1 diff --git a/bin/get-cutouts b/bin/get-cutouts index 32b9c1ce..d93c22d7 100755 --- a/bin/get-cutouts +++ b/bin/get-cutouts @@ -135,10 +135,9 @@ def plan(comm=None, specprod=None, coadd_type='healpix', else: rank, size = comm.rank, comm.size - desi_root = '/global/cfs/cdirs/desi' - #desi_root = os.environ.get('DESI_ROOT', DESI_ROOT_NERSC) - # look for data in the standard location + desi_root = os.environ.get('DESI_ROOT', '/dvs_ro/cfs/cdirs/desi') + # look for data in the standard location if coadd_type == 'healpix': subdir = 'healpix' if healpix is not None: @@ -352,7 +351,7 @@ def main(): parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument('--coadd-type', type=str, default='healpix', choices=['healpix', 'cumulative', 'pernight', 'perexp'], help='Specify which type of spectra/zbest files to process.') - parser.add_argument('--specprod', type=str, default='fuji', #choices=['fuji', 'guadalupe', 'iron'], + parser.add_argument('--specprod', type=str, default='iron', #choices=['fuji', 'guadalupe', 'iron'], help='Spectroscopic production to process.') parser.add_argument('--healpix', type=str, default=None, help='Comma-separated list of healpixels to process.') @@ -367,13 +366,14 @@ def main(): parser.add_argument('--nompi', action='store_true', help='Do not use MPI parallelism.') parser.add_argument('--dry-run', action='store_true', help='Generate but do not run commands.') - parser.add_argument('--outdir-html', default='/global/cfs/cdirs/desi/users/ioannis/fastspecfit', - type=str, help='Base output HTML directory.') - parser.add_argument('--outdir-data', default='/global/cfs/cdirs/desi/spectro/fastspecfit', - type=str, help='Base output data directory.') + parser.add_argument('--outdir-html', default='$PSCRATCH/fastspecfit/html', type=str, help='Base output HTML directory.') + parser.add_argument('--outdir-data', default='$PSCRATCH/fastspecfit/data', type=str, help='Base output data directory.') args = parser.parse_args() + outdir_data = os.path.expandvars(args.outdir_data) + outdir_html = os.path.expandvars(args.outdir_html) + if args.nompi: comm = None else: @@ -397,11 +397,10 @@ def main(): plan(comm=comm, specprod=args.specprod, coadd_type=args.coadd_type, survey=args.survey, program=args.program, healpix=args.healpix, tile=args.tile, night=args.night, - outdir_data=args.outdir_data, outdir_html=args.outdir_html, + outdir_data=outdir_data, outdir_html=outdir_html, overwrite=args.overwrite, mp=args.mp) else: - do_cutouts(args, comm=comm, outdir_data=args.outdir_data, - outdir_html=args.outdir_html) + do_cutouts(args, comm=comm, outdir_data=outdir_data, outdir_html=outdir_html) if __name__ == '__main__': main() diff --git a/bin/get-cutouts.sh b/bin/get-cutouts.sh deleted file mode 100755 index afefe32c..00000000 --- a/bin/get-cutouts.sh +++ /dev/null @@ -1,44 +0,0 @@ -#! /bin/bash - -#salloc -N 4 -C cpu -A desi -L cfs -t 04:00:00 --qos interactive --image=dstndstn/viewer-cutouts -#srun -n 4 -c 128 shifter /global/homes/i/ioannis/code/desihub/fastspecfit/bin/get-cutouts.sh fuji 128 > /global/cfs/cdirs/desi/spectro/fastspecfit/fuji/v2.0/logs/cutouts-fuji-01.log 2>&1 & - -codedir=/global/homes/i/ioannis/code/desihub -mpiscript=$codedir/fastspecfit/bin/get-cutouts - -outdir_data=/global/cfs/cdirs/desi/spectro/fastspecfit -outdir_html=/global/cfs/cdirs/desi/users/ioannis/fastspecfit - -export OMP_NUM_THREADS=1 -export MKL_NUM_THREADS=1 -export KMP_AFFINITY=disabled -export MPICH_GNI_FORK_MODE=FULLCOPY - -specprod=$1 -mp=$2 -coadd_type=$3 -survey=$4 -program=$5 - -args="--outdir-data $outdir_data --outdir-html $outdir_html" - -if [[ $specprod != " " ]] && [[ $specprod != "" ]] && [[ $specprod != "-" ]]; then - args=$args" --specprod $specprod" -fi -if [[ $mp != " " ]] && [[ $mp != "" ]] && [[ $mp != "-" ]]; then - args=$args" --mp $mp" -fi -if [[ $coadd_type != " " ]] && [[ $coadd_type != "" ]] && [[ $coadd_type != "-" ]]; then - args=$args" --coadd-type $coadd_type" -else - args=$args" --coadd-type healpix" -fi -if [[ $survey != " " ]] && [[ $survey != "" ]] && [[ $survey != "-" ]]; then - args=$args" --survey $survey" -fi -if [[ ! -z $program ]] && [[ $program != "" ]] && [[ $program != "-" ]]; then - args=$args" --program $program" -fi - -echo $mpiscript $args -time $mpiscript $args diff --git a/bin/get-cutouts.slurm b/bin/get-cutouts.slurm deleted file mode 100644 index 950ba655..00000000 --- a/bin/get-cutouts.slurm +++ /dev/null @@ -1,15 +0,0 @@ -#! /bin/bash -#SBATCH -A desi -#SBATCH -C cpu -#SBATCH -o /global/cfs/cdirs/desi/spectro/fastspecfit/fuji/v2.0/logs/cutouts-fuji-07.log.%j -#SBATCH --mail-user=jmoustakas@siena.edu -#SBATCH --image=dstndstn/viewer-cutouts -#SBATCH --mail-type=ALL -#SBATCH -q regular -#SBATCH -N 2 -#SBATCH -n 2 -#SBATCH -t 04:00:00 - -# sbatch /global/homes/i/ioannis/code/desihub/fastspecfit/bin/get-cutouts.slurm - -time srun -n 2 -c 128 /global/homes/i/ioannis/code/desihub/fastspecfit/bin/get-cutouts.sh fuji 128 healpix sv3 dark diff --git a/bin/mpi-fastspecfit b/bin/mpi-fastspecfit index 78ef48a6..249e6a18 100755 --- a/bin/mpi-fastspecfit +++ b/bin/mpi-fastspecfit @@ -26,8 +26,9 @@ log = get_logger() def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, makeqa=False, outdir_data='.', outdir_html='.'): - import sys, subprocess - from fastspecfit.mpi import backup_logs, plan + import sys + from desispec.parallel import stdouterr_redirected + from fastspecfit.mpi import plan if comm is None: rank, size = 0, 1 @@ -36,7 +37,6 @@ def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, t0 = time.time() if rank == 0: - #log.info('Starting at {}'.format(time.asctime())) _, zbestfiles, outfiles, groups, ntargets = plan( comm=comm, specprod=args.specprod, specprod_dir=specprod_dir, coadd_type=args.coadd_type, survey=args.survey, program=args.program, @@ -54,8 +54,6 @@ def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, groups = comm.bcast(groups, root=0) ntargets = comm.bcast(ntargets, root=0) - #if comm: - # comm.barrier() sys.stdout.flush() # all done @@ -65,45 +63,51 @@ def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, assert(len(groups) == size) assert(len(np.concatenate(groups)) == len(zbestfiles)) - #pixels = np.array([int(os.path.basename(os.path.dirname(x))) for x in zbestfiles]) for ii in groups[rank]: - log.debug('Rank {} started at {}'.format(rank, time.asctime())) + log.debug(f'Rank {rank} started at {time.asctime()}') sys.stdout.flush() # With --makeqa the desired output directories are in the 'zbestfiles'. if args.makeqa: - cmd = 'fastspecfit-qa {} -o {} --mp {}'.format(outfiles[ii], zbestfiles[ii], args.mp) - if args.ntargets: - cmd += ' --ntargets {}'.format(args.ntargets) + from fastspecfit.qa import fastqa as fast + cmd = 'fastqa' + cmdargs = f'{outfiles[ii]} -o {zbestfiles[ii]} --mp {args.mp}' else: if fastphot: - cmd = 'fastphot {} -o {} --mp {}'.format(zbestfiles[ii], outfiles[ii], args.mp) + from fastspecfit.fastspecfit import fastphot as fast + cmd = 'fastphot' else: - cmd = 'fastspec {} -o {} --mp {}'.format(zbestfiles[ii], outfiles[ii], args.mp) + from fastspecfit.fastspecfit import fastspec as fast + cmd = 'fastspec' + cmdargs = f'{zbestfiles[ii]} -o {outfiles[ii]} --mp {args.mp}' - if args.ntargets is not None: - cmd += ' --ntargets {}'.format(args.ntargets) + if args.ignore_quasarnet: + cmdargs += ' --ignore-quasarnet' - if args.targetids is not None: - cmd += ' --targetids {}'.format(args.targetids) + if args.ntargets is not None: + cmdargs += f' --ntargets {args.ntargets}' - if args.ignore_quasarnet: - cmd += ' --ignore-quasarnet' + if args.firsttarget is not None: + cmdargs += f' --firsttarget {args.firsttarget}' + + if args.targetids is not None: + cmdargs += f' --targetids {args.targetids}' if args.makeqa: logfile = os.path.join(zbestfiles[ii], os.path.basename(outfiles[ii]).replace('.gz', '').replace('.fits', '.log')) else: logfile = outfiles[ii].replace('.gz', '').replace('.fits', '.log') - assert(logfile != outfiles[ii]) + + if logfile == outfiles[ii]: + errmsg = f'Log file {logfile} cannot be the same as outfile!' + self.log.critical(errmsg) + raise ValueError(errmsg) if args.makeqa and args.plan: - log.info('Still need to make {} QA figure(s) across {} input file(s)'.format( - np.sum(ntargets), len(groups))) + log.info(f'Still need to make {np.sum(ntargets)} QA figure(s) across {len(groups)} input file(s)') return - log.info('Rank {}, ntargets={}: {}'.format(rank, ntargets[ii], cmd)) - #log.info(' rank {}: {}'.format(rank, cmd)) - #log.info('LOGGING to {}'.format(logfile)) + log.info(f'Rank {rank}, ntargets={ntargets[ii]}: {cmd} {cmdargs}') sys.stdout.flush() if args.dry_run: @@ -111,30 +115,25 @@ def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, try: t1 = time.time() - if os.path.exists(logfile) and not args.overwrite: - backup_logs(logfile) - # memory leak? Try making system call instead outdir = os.path.dirname(logfile) if not os.path.isdir(outdir): os.makedirs(outdir, exist_ok=True) if args.nolog: - err = subprocess.call(cmd.split()) + fast(args=cmdargs.split()) else: - with open(logfile, 'w') as mylog: - err = subprocess.call(cmd.split(), stdout=mylog, stderr=mylog) + with stdouterr_redirected(to=logfile, overwrite=args.overwrite): + fast(args=cmdargs.split()) dt1 = time.time() - t1 - if err == 0: - log.info(' rank {} done in {:.2f} sec'.format(rank, dt1)) - if not os.path.exists(outfiles[ii]): - log.warning(' rank {} missing {}'.format(rank, outfiles[ii])) - else: - log.warning(' rank {} broke after {:.1f} sec with error code {}'.format(rank, dt1, err)) - except Exception as err: - log.warning(' rank {} raised an exception'.format(rank)) + log.info(' rank {} done in {:.2f} sec'.format(rank, dt1)) + if not os.path.exists(outfiles[ii]): + log.warning(f' rank {rank} missing {outfiles[ii]}') + raise IOError + except: + log.warning(f' rank {rank} raised an exception') import traceback traceback.print_exc() - log.debug(' rank {} is done'.format(rank)) + log.debug(f' rank {rank} is done') sys.stdout.flush() if comm is not None: @@ -143,9 +142,9 @@ def run_fastspecfit(args, comm=None, fastphot=False, specprod_dir=None, if rank == 0 and not args.dry_run: for outfile in outfiles: if not os.path.exists(outfile): - log.warning('Missing {}'.format(outfile)) + log.warning(f'Missing {outfile}') - log.info('All done at {}'.format(time.asctime())) + log.info(f'All done at {time.asctime()}') def main(): """Main wrapper on fastphot and fastspec. @@ -159,7 +158,7 @@ def main(): parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument('--coadd-type', type=str, default='healpix', choices=['healpix', 'cumulative', 'pernight', 'perexp'], help='Specify which type of spectra/zbest files to process.') - parser.add_argument('--specprod', type=str, default='fuji', #choices=['fuji', 'guadalupe', 'iron'], + parser.add_argument('--specprod', type=str, default='iron', #choices=['fuji', 'guadalupe', 'iron'], help='Spectroscopic production to process.') parser.add_argument('--healpix', type=str, default=None, help='Comma-separated list of healpixels to process.') @@ -170,6 +169,7 @@ def main(): parser.add_argument('--mp', type=int, default=1, help='Number of multiprocessing processes per MPI rank or node.') parser.add_argument('-n', '--ntargets', type=int, help='Number of targets to process in each file.') + parser.add_argument('--firsttarget', type=int, default=0, help='Index of first object to to process in each file, zero-indexed.') parser.add_argument('--targetids', type=str, default=None, help='Comma-separated list of TARGETIDs to process.') parser.add_argument('--fastphot', action='store_true', help='Fit the broadband photometry.') @@ -189,15 +189,16 @@ def main(): parser.add_argument('--nolog', action='store_true', help='Do not write to the log file.') parser.add_argument('--dry-run', action='store_true', help='Generate but do not run commands.') - parser.add_argument('--outdir-html', default='/global/cfs/cdirs/desi/users/ioannis/fastspecfit', - type=str, help='Base output HTML directory.') - parser.add_argument('--outdir-data', default='/global/cfs/cdirs/desi/spectro/fastspecfit', - type=str, help='Base output data directory.') + parser.add_argument('--outdir-html', default='$PSCRATCH/fastspecfit/html', type=str, help='Base output HTML directory.') + parser.add_argument('--outdir-data', default='$PSCRATCH/fastspecfit/data', type=str, help='Base output data directory.') specprod_dir = None args = parser.parse_args() + outdir_data = os.path.expandvars(args.outdir_data) + outdir_html = os.path.expandvars(args.outdir_html) + if args.merge or args.mergeall or args.nompi: comm = None else: @@ -207,6 +208,11 @@ def main(): except ImportError: comm = None + # https://docs.nersc.gov/development/languages/python/parallel-python/#use-the-spawn-start-method + if args.mp > 1 and 'NERSC_HOST' in os.environ: + import multiprocessing + multiprocessing.set_start_method('spawn') + if args.coadd_type == 'healpix': args.survey = args.survey.split(',') args.program = args.program.split(',') @@ -215,7 +221,7 @@ def main(): if args.ntest: from astropy.table import Table rand = np.random.RandomState(seed=1) - tilepix = Table.read('/global/cfs/cdirs/desi/spectro/redux/{}/healpix/tilepix.fits'.format(args.specprod)) + tilepix = Table.read(f'$DESI_ROOT/spectro/redux/{args.specprod}/healpix/tilepix.fits') # trim to survey/program tilepix = tilepix[np.unique(np.hstack([np.where(tilepix['SURVEY'] == survey)[0] for survey in args.survey]))] @@ -230,14 +236,14 @@ def main(): survprog = np.array([ss+pp for ss, pp in zip(tilepix['SURVEY'], tilepix['PROGRAM'])]) I = rand.choice(len(tilepix), args.ntest, replace=False) args.healpix = ','.join(tilepix['HEALPIX'][I].astype(str)) - log.info('Selecting {} test healpixels: {}'.format(len(I), args.healpix)) + log.info(f'Selecting {len(I)} test healpixels: {args.healpix}') print(tilepix[I]) if args.merge or args.mergeall: from fastspecfit.mpi import merge_fastspecfit merge_fastspecfit(specprod=args.specprod, specprod_dir=specprod_dir, coadd_type=args.coadd_type, survey=args.survey, program=args.program, healpix=args.healpix, - tile=args.tile, night=args.night, outdir_data=args.outdir_data, + tile=args.tile, night=args.night, outdir_data=outdir_data, outsuffix=args.merge_suffix, mergedir=args.mergedir, overwrite=args.overwrite, fastphot=args.fastphot, supermerge=args.mergeall, mp=args.mp) return @@ -252,12 +258,11 @@ def main(): coadd_type=args.coadd_type, survey=args.survey, program=args.program, healpix=args.healpix, tile=args.tile, night=args.night, makeqa=args.makeqa, mp=args.mp, - fastphot=args.fastphot, outdir_data=args.outdir_data, - outdir_html=args.outdir_html, overwrite=args.overwrite) + fastphot=args.fastphot, outdir_data=outdir_data, + outdir_html=outdir_html, overwrite=args.overwrite) else: run_fastspecfit(args, comm=comm, fastphot=args.fastphot, specprod_dir=specprod_dir, - makeqa=args.makeqa, outdir_data=args.outdir_data, - outdir_html=args.outdir_html) + makeqa=args.makeqa, outdir_data=outdir_data, outdir_html=outdir_html) if __name__ == '__main__': main() diff --git a/bin/mpi-fastspecfit.sh b/bin/mpi-fastspecfit.sh deleted file mode 100755 index 2b5bde45..00000000 --- a/bin/mpi-fastspecfit.sh +++ /dev/null @@ -1,117 +0,0 @@ -#! /bin/bash - -# Shell script for running mpi-fastspecfit with MPI+shifter at NERSC. Required arguments: -# {1} stage [fastspec, fastphot, test] -# {2} mp [should match the resources requested.] - -# Example: build the coadds using 16 MPI tasks with 8 cores per node (and therefore 16*8/32=4 nodes) - -# Perlmutter w/o shifter -#salloc -N 4 -C cpu -A desi -L cfs -t 04:00:00 --qos interactive - -# Testing-- -#srun -n 4 -c 128 /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh fastspec-test iron 128 > /global/cfs/cdirs/desi/spectro/fastspecfit/iron/logs/fastspec-test.log.1 2>&1 & - -# Perlmutter w/ shifter -#salloc -N 4 -C cpu -A desi -L cfs -t 04:00:00 --qos interactive --image=docker:desihub/fastspecfit:2.1.1 -#srun -n 4 -c 128 shifter --module=mpich /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh fastspec iron 128 healpix cmx,sv1,sv2,sv3,special > /global/cfs/cdirs/desi/spectro/fastspecfit/iron/logs/fastspec-iron-sv1.log.1 2>&1 & - -#srun -n 16 -c 32 shifter --module=mpich /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh fastspec 32 sv2 dark > /global/cfs/cdirs/desi/spectro/fastspecfit/iron/logs/fastspec-iron-sv2-dark.log.1 2>&1 & -#srun -n 16 -c 32 shifter --module=mpich /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh fastphot 32 sv2 dark > /global/cfs/cdirs/desi/spectro/fastspecfit/iron/logs/fastphot-iron-sv2-dark.log.1 2>&1 & -#srun -n 16 -c 32 shifter --module=mpich /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh qafastspec 32 > qafastspec-iron.log.1 2>&1 & -#srun -n 16 -c 32 shifter --module=mpich /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh qafastphot 32 > qafastphot-iron.log.1 2>&1 & - -# Cori -#salloc -N 32 -C haswell -A desi -L cfs -t 04:00:00 --qos interactive --image=docker:desihub/fastspecfit:2.1.1 -#srun -n 32 -c 32 shifter --module=mpich /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh fastspec iron 32 healpix sv1 > /global/cfs/cdirs/desi/spectro/fastspecfit/iron/logs/fastspec-iron-sv2sv3.log.1 2>&1 & - -#srun -n 16 -c 32 shifter --module=mpich-cle6 /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh fastspec iron 32 - sv1 > /global/cfs/cdirs/desi/spectro/fastspecfit/iron/logs/fastspec-iron-sv1.log.1 2>&1 & -#srun -n 16 -c 32 shifter --module=mpich-cle6 /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh fastphot 32 sv2 dark > /global/cfs/cdirs/desi/spectro/fastspecfit/iron/logs/fastphot-iron-sv2-dark.log.1 2>&1 & -#srun -n 16 -c 32 shifter --module=mpich-cle6 /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh qafastspec 32 > qafastspec-iron.log.1 2>&1 & -#srun -n 16 -c 32 shifter --module=mpich-cle6 /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh qafastphot 32 > qafastphot-iron.log.1 2>&1 & - -codedir=/global/homes/i/ioannis/code/desihub -#codedir=/usr/local/bin -mpiscript=$codedir/fastspecfit/bin/mpi-fastspecfit - -echo 'Loading DESI software stack 23.1' -source /global/cfs/cdirs/desi/software/desi_environment.sh 23.1 -module swap desispec/0.57.0 -module load fastspecfit/2.1.1 - -#for package in desispec fastspecfit; do -# echo Loading local check-out of $package -# export PATH=$codedir/$package/bin:$PATH -# export PYTHONPATH=$codedir/$package/py:$PYTHONPATH -#done - -outdir_data=/global/cfs/cdirs/desi/spectro/fastspecfit -outdir_html=/global/cfs/cdirs/desi/users/ioannis/fastspecfit - -export DESI_ROOT='/global/cfs/cdirs/desi' -export DUST_DIR='/global/cfs/cdirs/cosmo/data/dust/v0_1' -export DR9_DIR='/global/cfs/cdirs/desi/external/legacysurvey/dr9' -export FTEMPLATES_DIR='/global/cfs/cdirs/desi/external/templates/fastspecfit' - -export TMPCACHE=$(mktemp -d) -export MPLCONFIGDIR=$TMPCACHE/matplotlib -mkdir $MPLCONFIGDIR -cp -r $HOME/.config/matplotlib $MPLCONFIGDIR - -export OMP_NUM_THREADS=1 -export MKL_NUM_THREADS=1 -export KMP_AFFINITY=disabled -export MPICH_GNI_FORK_MODE=FULLCOPY - -stage=$1 -specprod=$2 -mp=$3 -coadd_type=$4 -survey=$5 -program=$6 - -#echo stage=$stage -#echo specprod=$specprod -#echo mp=$mp -#echo coadd_type=$coadd_type -#echo survey=$survey -#echo program=$program - -# petal 8 on tiles 80613, 80606, 80607 -#args="--outdir-data $outdir_data --outdir-html $outdir_html --healpix 17684,17685,17686,17687,17692,7015,7020,7021,7022,7023,7026,7032,7015,7020,7021,7022,7023,7032" - -args="--outdir-data $outdir_data --outdir-html $outdir_html" - -if [[ $stage == "fastspec-test" ]]; then - args=$args" --ntest 12" -fi -if [[ $stage == "fastphot-test" ]]; then - args=$args" --fastphot --ntest 100" -fi -if [[ $stage == "fastphot" ]]; then - args=$args" --fastphot" -fi -if [[ $stage == "makeqa" ]]; then - args=$args" --makeqa" -fi - -if [[ $specprod != " " ]] && [[ $specprod != "" ]] && [[ $specprod != "-" ]]; then - args=$args" --specprod $specprod" -fi -if [[ $mp != " " ]] && [[ $mp != "" ]] && [[ $mp != "-" ]]; then - args=$args" --mp $mp" -fi -if [[ $coadd_type != " " ]] && [[ $coadd_type != "" ]] && [[ $coadd_type != "-" ]]; then - args=$args" --coadd-type $coadd_type" -else - args=$args" --coadd-type healpix" -fi -if [[ $survey != " " ]] && [[ $survey != "" ]] && [[ $survey != "-" ]]; then - args=$args" --survey $survey" -fi -if [[ ! -z $program ]] && [[ $program != "" ]] && [[ $program != "-" ]]; then - args=$args" --program $program" -fi - -echo python $mpiscript $args -time python $mpiscript $args diff --git a/bin/mpi-fastspecfit.slurm b/bin/mpi-fastspecfit.slurm deleted file mode 100644 index 12690545..00000000 --- a/bin/mpi-fastspecfit.slurm +++ /dev/null @@ -1,15 +0,0 @@ -#! /bin/bash -#SBATCH -A desi -#SBATCH -C cpu -#SBATCH -o /global/cfs/cdirs/desi/spectro/fastspecfit/fuji/v2.0/logs/fastspec-fuji-qa-03.log.%j -#SBATCH --mail-user=jmoustakas@siena.edu -#SBATCH --mail-type=ALL -#SBATCH -q regular -#SBATCH -N 32 -#SBATCH -n 32 -#SBATCH -t 02:00:00 - -# sbatch /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.slurm - -time srun -n 32 -c 128 /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh makeqa fuji 128 healpix -#time srun -n 32 -c 128 /global/homes/i/ioannis/code/desihub/fastspecfit/bin/mpi-fastspecfit.sh fastspec iron 128 healpix main bright diff --git a/doc/api.rst b/doc/api.rst index 98f844d3..0d703180 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -25,6 +25,9 @@ API .. automodule:: fastspecfit.mpi :members: +.. automodule:: fastspecfit.qa + :members: + .. automodule:: fastspecfit.templates.qa :members: diff --git a/doc/changes.rst b/doc/changes.rst index 8a732d55..dd4645b3 100644 --- a/doc/changes.rst +++ b/doc/changes.rst @@ -6,8 +6,10 @@ Change Log ------------------------ * Support non-DR9 photometry via a new photometric configuration file [`PR #133`_]. +* Miscellaneous updates and bug fixes ahead of next version of VACs [`PR #145`_]. .. _`PR #133`: https://github.com/desihub/fastspecfit/pull/133 +.. _`PR #145`: https://github.com/desihub/fastspecfit/pull/145 2.2.0 (2023-08-02) ------------------ diff --git a/doc/fastphot.rst b/doc/fastphot.rst index c8585a54..3c979a11 100644 --- a/doc/fastphot.rst +++ b/doc/fastphot.rst @@ -107,39 +107,49 @@ Name Type Units Description FLUX_SYNTH_PHOTMODEL_W2 float32 nmgy W2-band flux synthesized from the best-fitting photometric continuum model. FLUX_SYNTH_PHOTMODEL_W3 float32 nmgy W3-band flux synthesized from the best-fitting photometric continuum model. FLUX_SYNTH_PHOTMODEL_W4 float32 nmgy W4-band flux synthesized from the best-fitting photometric continuum model. - KCORR_U float32 mag K-correction used to derive ABSMAG_U band-shifted to z=0.0 assuming h=1.0. - ABSMAG_U float32 mag Absolute magnitude in Johnson/Cousins U-band band-shifted to z=0.0 assuming h=1.0. - ABSMAG_IVAR_U float32 1 / mag2 Inverse variance corresponding to ABSMAG_U. - KCORR_B float32 mag Like KCORR_U but for Johnson/Cousins B-band. - ABSMAG_B float32 mag Like ABSMAG_U but for Johnson/Cousins B-band. - ABSMAG_IVAR_B float32 1 / mag2 Like ABSMAG_IVAR_U but for Johnson/Cousins B-band. - KCORR_V float32 mag Like KCORR_U but for Johnson/Cousins V-band. - ABSMAG_V float32 mag Like ABSMAG_U but for Johnson/Cousins V-band. - ABSMAG_IVAR_V float32 1 / mag2 Like ABSMAG_IVAR_U but for Johnson/Cousins V-band. - KCORR_SDSS_U float32 mag K-correction used to derive ABSMAG_SDSS_U band-shifted to z=0.1 assuming h=1.0. - ABSMAG_SDSS_U float32 mag Absolute magnitude in SDSS u-band band-shifted to z=0.1 assuming h=1.0. - ABSMAG_IVAR_SDSS_U float32 1 / mag2 Inverse variance corresponding to ABSMAG_SDSS_U. - KCORR_SDSS_G float32 mag Like KCORR_SDSS_U but for SDSS g-band. - ABSMAG_SDSS_G float32 mag Like ABSMAG_SDSS_U but for SDSS g-band. - ABSMAG_IVAR_SDSS_G float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS g-band. - KCORR_SDSS_R float32 mag Like KCORR_SDSS_U but for SDSS r-band. - ABSMAG_SDSS_R float32 mag Like ABSMAG_SDSS_U but for SDSS r-band. - ABSMAG_IVAR_SDSS_R float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS r-band. - KCORR_SDSS_I float32 mag Like KCORR_SDSS_U but for SDSS i-band. - ABSMAG_SDSS_I float32 mag Like ABSMAG_SDSS_U but for SDSS i-band. - ABSMAG_IVAR_SDSS_I float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS i-band. - KCORR_SDSS_Z float32 mag Like KCORR_SDSS_U but for SDSS z-band. - ABSMAG_SDSS_Z float32 mag Like ABSMAG_SDSS_U but for SDSS z-band. - ABSMAG_IVAR_SDSS_Z float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS z-band. - KCORR_W1 float32 mag K-correction used to derive ABSMAG_W1 band-shifted to z=0.0 assuming h=1.0. - ABSMAG_W1 float32 mag Absolute magnitude in WISE W1-band band-shifted to z=0.0 assuming h=1.0. - ABSMAG_IVAR_W1 float32 1 / mag2 Inverse variance corresponding to ABSMAG_W1. - KCORR_W2 float32 mag K-correction used to derive ABSMAG_W2 band-shifted to z=0.0 assuming h=1.0. - ABSMAG_W2 float32 mag Absolute magnitude in WISE W2-band band-shifted to z=0.0 assuming h=1.0. - ABSMAG_IVAR_W2 float32 1 / mag2 Inverse variance corresponding to ABSMAG_W2. + ABSMAG10_DECAM_G [4]_ float32 mag Absolute magnitude in DECam g-band band-shifted to z=1.0 assuming h=1.0. + ABSMAG10_IVAR_DECAM_G float32 1 / mag2 Inverse variance corresponding to ABSMAG10_DECAM_G. + KCORR10_DECAM_G float32 mag K-correction used to derive ABSMAG10_DECAM_G band-shifted to z=1.0. + ABSMAG10_DECAM_R [4]_ float32 mag Absolute magnitude in DECam r-band band-shifted to z=1.0 assuming h=1.0. + ABSMAG10_IVAR_DECAM_R float32 1 / mag2 Inverse variance corresponding to ABSMAG10_DECAM_R. + KCORR10_DECAM_R float32 mag K-correction used to derive ABSMAG10_DECAM_R band-shifted to z=1.0. + ABSMAG10_DECAM_Z [4]_ float32 mag Absolute magnitude in DECam z-band band-shifted to z=1.0 assuming h=1.0. + ABSMAG10_IVAR_DECAM_Z float32 1 / mag2 Inverse variance corresponding to ABSMAG10_DECAM_Z. + KCORR10_DECAM_Z float32 mag K-correction used to derive ABSMAG10_DECAM_Z band-shifted to z=1.0. + ABSMAG00_U [4]_ float32 mag Absolute magnitude in Johnson/Cousins U-band band-shifted to z=0.0 assuming h=1.0. + ABSMAG00_IVAR_U float32 1 / mag2 Inverse variance corresponding to ABSMAG_U. + KCORR00_U float32 mag K-correction used to derive ABSMAG_U band-shifted to z=0.0. + ABSMAG00_B [4]_ float32 mag Like ABSMAG_U but for Johnson/Cousins B-band. + ABSMAG00_IVAR_B float32 1 / mag2 Like ABSMAG_IVAR_U but for Johnson/Cousins B-band. + KCORR00_B float32 mag Like KCORR_U but for Johnson/Cousins B-band. + ABSMAG00_V [4]_ float32 mag Like ABSMAG_U but for Johnson/Cousins V-band. + ABSMAG00_IVAR_V float32 1 / mag2 Like ABSMAG_IVAR_U but for Johnson/Cousins V-band. + KCORR00_V float32 mag Like KCORR_U but for Johnson/Cousins V-band. + ABSMAG01_SDSS_U [4]_ float32 mag Absolute magnitude in SDSS u-band band-shifted to z=0.1 assuming h=1.0. + ABSMAG01_IVAR_SDSS_U float32 1 / mag2 Inverse variance corresponding to ABSMAG_SDSS_U. + KCORR01_SDSS_U float32 mag K-correction used to derive ABSMAG_SDSS_U band-shifted to z=0.1. + ABSMAG01_SDSS_G [4]_ float32 mag Like ABSMAG_SDSS_U but for SDSS g-band. + ABSMAG01_IVAR_SDSS_G float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS g-band. + KCORR01_SDSS_G float32 mag Like KCORR_SDSS_U but for SDSS g-band. + ABSMAG01_SDSS_R [4]_ float32 mag Like ABSMAG_SDSS_U but for SDSS r-band. + ABSMAG01_IVAR_SDSS_R float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS r-band. + KCORR01_SDSS_R float32 mag Like KCORR_SDSS_U but for SDSS r-band. + ABSMAG01_SDSS_I [4]_ float32 mag Like ABSMAG_SDSS_U but for SDSS i-band. + ABSMAG01_IVAR_SDSS_I float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS i-band. + KCORR01_SDSS_I float32 mag Like KCORR_SDSS_U but for SDSS i-band. + ABSMAG01_SDSS_Z [4]_ float32 mag Like ABSMAG_SDSS_U but for SDSS z-band. + ABSMAG01_IVAR_SDSS_Z float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS z-band. + KCORR01_SDSS_Z float32 mag Like KCORR_SDSS_U but for SDSS z-band. + ABSMAG00_W1 [4]_ float32 mag Absolute magnitude in WISE W1-band band-shifted to z=0.0 assuming h=1.0. + ABSMAG00_IVAR_W1 float32 1 / mag2 Inverse variance corresponding to ABSMAG_W1. + KCORR00_W1 float32 mag K-correction used to derive ABSMAG_W1 band-shifted to z=0.0. LOGLNU_1500 float32 1e-28 erg / (s Hz) Monochromatic luminosity at 1500 A in the rest-frame. LOGLNU_2800 float32 1e-28 erg / (s Hz) Monochromatic luminosity at 2800 A in the rest-frame. - LOGL_5100 float32 1e+10 Lsun Total luminosity at 5100 A in the rest-frame. + LOGL_1450 float32 1e+10 Lsun Integrated luminosity at 1450 A in the rest-frame. + LOGL_1700 float32 1e+10 Lsun Integrated luminosity at 1700 A in the rest-frame. + LOGL_3000 float32 1e+10 Lsun Integrated luminosity at 3000 A in the rest-frame. + LOGL_5100 float32 1e+10 Lsun Integrated luminosity at 5100 A in the rest-frame. + FLYA_1215_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 1215.67 A in the rest-frame. FOII_3727_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 3728.483 A in the rest-frame. FHBETA_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 4862.683 A in the rest-frame. FOIII_5007_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 5008.239 A in the rest-frame. @@ -195,18 +205,18 @@ Name Type Units Description BGS_TARGET int64 BGS targeting bit. MWS_TARGET int64 MWS targeting bit. SCND_TARGET int64 Secondary target targeting bit. - SV1_DESI_TARGET [4]_ int64 SV1 DESI targeting bit. - SV1_BGS_TARGET [4]_ int64 SV1 BGS targeting bit. - SV1_MWS_TARGET [4]_ int64 SV1 MWS targeting bit. - SV2_DESI_TARGET [4]_ int64 SV2 DESI targeting bit. - SV2_BGS_TARGET [4]_ int64 SV2 BGS targeting bit. - SV2_MWS_TARGET [4]_ int64 SV2 MWS targeting bit. - SV3_DESI_TARGET [4]_ int64 SV3 DESI targeting bit. - SV3_BGS_TARGET [4]_ int64 SV3 BGS targeting bit. - SV3_MWS_TARGET [4]_ int64 SV3 MWS targeting bit. - SV1_SCND_TARGET [4]_ int64 SV1 secondary targeting bit. - SV2_SCND_TARGET [4]_ int64 SV2 secondary targeting bit. - SV3_SCND_TARGET [4]_ int64 SV3 secondary targeting bit. + SV1_DESI_TARGET [5]_ int64 SV1 DESI targeting bit. + SV1_BGS_TARGET [5]_ int64 SV1 BGS targeting bit. + SV1_MWS_TARGET [5]_ int64 SV1 MWS targeting bit. + SV2_DESI_TARGET [5]_ int64 SV2 DESI targeting bit. + SV2_BGS_TARGET [5]_ int64 SV2 BGS targeting bit. + SV2_MWS_TARGET [5]_ int64 SV2 MWS targeting bit. + SV3_DESI_TARGET [5]_ int64 SV3 DESI targeting bit. + SV3_BGS_TARGET [5]_ int64 SV3 BGS targeting bit. + SV3_MWS_TARGET [5]_ int64 SV3 MWS targeting bit. + SV1_SCND_TARGET [5]_ int64 SV1 secondary targeting bit. + SV2_SCND_TARGET [5]_ int64 SV2 secondary targeting bit. + SV3_SCND_TARGET [5]_ int64 SV3 secondary targeting bit. Z float64 Redshift based on Redrock or QuasarNet (for QSO targets only). ZWARN int8 Redrock zwarning bit. DELTACHI2 float64 Redrock delta-chi-squared. @@ -255,4 +265,11 @@ Name Type Units Description .. [3] Column only present when fitting per-exposure tile-based coadds. -.. [4] Column only present in Commissioning and Survey Validation spectroscopic observations. +.. [4] Only observed photometry with a minimum signal-to-noise ratio of two is + used to compute K-corrections. Absolute magnitudes are estimated using + from the (S/N>2) observed-frame bandpass nearest in wavelength to the + desired band-shifted rest-frame bandpass. If no observed-frame photometry + is available, then the absolute magnitude is synthesized from the + best-fitting model and the corresponding inverse variance is set to zero. + +.. [5] Column only present in Commissioning and Survey Validation spectroscopic observations. diff --git a/doc/fastspec.rst b/doc/fastspec.rst index 8ea56942..4433947e 100644 --- a/doc/fastspec.rst +++ b/doc/fastspec.rst @@ -126,39 +126,49 @@ Name Type Units Descript FLUX_SYNTH_PHOTMODEL_W2 float32 nmgy W2-band flux synthesized from the best-fitting photometric continuum model. FLUX_SYNTH_PHOTMODEL_W3 float32 nmgy W3-band flux synthesized from the best-fitting photometric continuum model. FLUX_SYNTH_PHOTMODEL_W4 float32 nmgy W4-band flux synthesized from the best-fitting photometric continuum model. - KCORR_U float32 mag K-correction used to derive ABSMAG_U band-shifted to z=0.0 assuming h=1.0. - ABSMAG_U float32 mag Absolute magnitude in Johnson/Cousins U-band band-shifted to z=0.0 assuming h=1.0. - ABSMAG_IVAR_U float32 1 / mag2 Inverse variance corresponding to ABSMAG_U. - KCORR_B float32 mag Like KCORR_U but for Johnson/Cousins B-band. - ABSMAG_B float32 mag Like ABSMAG_U but for Johnson/Cousins B-band. - ABSMAG_IVAR_B float32 1 / mag2 Like ABSMAG_IVAR_U but for Johnson/Cousins B-band. - KCORR_V float32 mag Like KCORR_U but for Johnson/Cousins V-band. - ABSMAG_V float32 mag Like ABSMAG_U but for Johnson/Cousins V-band. - ABSMAG_IVAR_V float32 1 / mag2 Like ABSMAG_IVAR_U but for Johnson/Cousins V-band. - KCORR_SDSS_U float32 mag K-correction used to derive ABSMAG_SDSS_U band-shifted to z=0.1 assuming h=1.0. - ABSMAG_SDSS_U float32 mag Absolute magnitude in SDSS u-band band-shifted to z=0.1 assuming h=1.0. - ABSMAG_IVAR_SDSS_U float32 1 / mag2 Inverse variance corresponding to ABSMAG_SDSS_U. - KCORR_SDSS_G float32 mag Like KCORR_SDSS_U but for SDSS g-band. - ABSMAG_SDSS_G float32 mag Like ABSMAG_SDSS_U but for SDSS g-band. - ABSMAG_IVAR_SDSS_G float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS g-band. - KCORR_SDSS_R float32 mag Like KCORR_SDSS_U but for SDSS r-band. - ABSMAG_SDSS_R float32 mag Like ABSMAG_SDSS_U but for SDSS r-band. - ABSMAG_IVAR_SDSS_R float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS r-band. - KCORR_SDSS_I float32 mag Like KCORR_SDSS_U but for SDSS i-band. - ABSMAG_SDSS_I float32 mag Like ABSMAG_SDSS_U but for SDSS i-band. - ABSMAG_IVAR_SDSS_I float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS i-band. - KCORR_SDSS_Z float32 mag Like KCORR_SDSS_U but for SDSS z-band. - ABSMAG_SDSS_Z float32 mag Like ABSMAG_SDSS_U but for SDSS z-band. - ABSMAG_IVAR_SDSS_Z float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS z-band. - KCORR_W1 float32 mag K-correction used to derive ABSMAG_W1 band-shifted to z=0.0 assuming h=1.0. - ABSMAG_W1 float32 mag Absolute magnitude in WISE W1-band band-shifted to z=0.0 assuming h=1.0. - ABSMAG_IVAR_W1 float32 1 / mag2 Inverse variance corresponding to ABSMAG_W1. - KCORR_W2 float32 mag K-correction used to derive ABSMAG_W2 band-shifted to z=0.0 assuming h=1.0. - ABSMAG_W2 float32 mag Absolute magnitude in WISE W2-band band-shifted to z=0.0 assuming h=1.0. - ABSMAG_IVAR_W2 float32 1 / mag2 Inverse variance corresponding to ABSMAG_W2. + ABSMAG10_DECAM_G [4]_ float32 mag Absolute magnitude in DECam g-band band-shifted to z=1.0 assuming h=1.0. + ABSMAG10_IVAR_DECAM_G float32 1 / mag2 Inverse variance corresponding to ABSMAG10_DECAM_G. + KCORR10_DECAM_G float32 mag K-correction used to derive ABSMAG10_DECAM_G band-shifted to z=1.0. + ABSMAG10_DECAM_R [4]_ float32 mag Absolute magnitude in DECam r-band band-shifted to z=1.0 assuming h=1.0. + ABSMAG10_IVAR_DECAM_R float32 1 / mag2 Inverse variance corresponding to ABSMAG10_DECAM_R. + KCORR10_DECAM_R float32 mag K-correction used to derive ABSMAG10_DECAM_R band-shifted to z=1.0. + ABSMAG10_DECAM_Z [4]_ float32 mag Absolute magnitude in DECam z-band band-shifted to z=1.0 assuming h=1.0. + ABSMAG10_IVAR_DECAM_Z float32 1 / mag2 Inverse variance corresponding to ABSMAG10_DECAM_Z. + KCORR10_DECAM_Z float32 mag K-correction used to derive ABSMAG10_DECAM_Z band-shifted to z=1.0. + ABSMAG00_U [4]_ float32 mag Absolute magnitude in Johnson/Cousins U-band band-shifted to z=0.0 assuming h=1.0. + ABSMAG00_IVAR_U float32 1 / mag2 Inverse variance corresponding to ABSMAG_U. + KCORR00_U float32 mag K-correction used to derive ABSMAG_U band-shifted to z=0.0. + ABSMAG00_B [4]_ float32 mag Like ABSMAG_U but for Johnson/Cousins B-band. + ABSMAG00_IVAR_B float32 1 / mag2 Like ABSMAG_IVAR_U but for Johnson/Cousins B-band. + KCORR00_B float32 mag Like KCORR_U but for Johnson/Cousins B-band. + ABSMAG00_V [4]_ float32 mag Like ABSMAG_U but for Johnson/Cousins V-band. + ABSMAG00_IVAR_V float32 1 / mag2 Like ABSMAG_IVAR_U but for Johnson/Cousins V-band. + KCORR00_V float32 mag Like KCORR_U but for Johnson/Cousins V-band. + ABSMAG01_SDSS_U [4]_ float32 mag Absolute magnitude in SDSS u-band band-shifted to z=0.1 assuming h=1.0. + ABSMAG01_IVAR_SDSS_U float32 1 / mag2 Inverse variance corresponding to ABSMAG_SDSS_U. + KCORR01_SDSS_U float32 mag K-correction used to derive ABSMAG_SDSS_U band-shifted to z=0.1. + ABSMAG01_SDSS_G [4]_ float32 mag Like ABSMAG_SDSS_U but for SDSS g-band. + ABSMAG01_IVAR_SDSS_G float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS g-band. + KCORR01_SDSS_G float32 mag Like KCORR_SDSS_U but for SDSS g-band. + ABSMAG01_SDSS_R [4]_ float32 mag Like ABSMAG_SDSS_U but for SDSS r-band. + ABSMAG01_IVAR_SDSS_R float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS r-band. + KCORR01_SDSS_R float32 mag Like KCORR_SDSS_U but for SDSS r-band. + ABSMAG01_SDSS_I [4]_ float32 mag Like ABSMAG_SDSS_U but for SDSS i-band. + ABSMAG01_IVAR_SDSS_I float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS i-band. + KCORR01_SDSS_I float32 mag Like KCORR_SDSS_U but for SDSS i-band. + ABSMAG01_SDSS_Z [4]_ float32 mag Like ABSMAG_SDSS_U but for SDSS z-band. + ABSMAG01_IVAR_SDSS_Z float32 1 / mag2 Like ABSMAG_IVAR_SDSS_U but for SDSS z-band. + KCORR01_SDSS_Z float32 mag Like KCORR_SDSS_U but for SDSS z-band. + ABSMAG01_W1 [4]_ float32 mag Absolute magnitude in WISE W1-band band-shifted to z=0.0 assuming h=1.0. + ABSMAG01_IVAR_W1 float32 1 / mag2 Inverse variance corresponding to ABSMAG_W1. + KCORR01_W1 float32 mag K-correction used to derive ABSMAG_W1 band-shifted to z=0.0. LOGLNU_1500 float32 1e-28 erg / (s Hz) Monochromatic luminosity at 1500 A in the rest-frame. LOGLNU_2800 float32 1e-28 erg / (s Hz) Monochromatic luminosity at 2800 A in the rest-frame. - LOGL_5100 float32 1e+10 Lsun Total luminosity at 5100 A in the rest-frame. + LOGL_1450 float32 1e+10 Lsun Integrated luminosity at 1450 A in the rest-frame. + LOGL_1700 float32 1e+10 Lsun Integrated luminosity at 1700 A in the rest-frame. + LOGL_3000 float32 1e+10 Lsun Integrated luminosity at 3000 A in the rest-frame. + LOGL_5100 float32 1e+10 Lsun Integrated luminosity at 5100 A in the rest-frame. + FLYA_1215_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 1215.67 A in the rest-frame. FOII_3727_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 3728.483 A in the rest-frame. FHBETA_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 4862.683 A in the rest-frame. FOIII_5007_CONT float32 1e-17 erg / (Angstrom cm2 s) Continuum flux at 5008.239 A in the rest-frame. @@ -926,23 +936,23 @@ Name Type Units Description RA float64 deg Right ascension from target catalog. DEC float64 deg Declination from target catalog. COADD_FIBERSTATUS int64 Fiber status bit. - CMX_TARGET [4]_ int64 Commissioning (CMX) targeting bit. + CMX_TARGET [5]_ int64 Commissioning (CMX) targeting bit. DESI_TARGET int64 DESI targeting bit. BGS_TARGET int64 BGS targeting bit. MWS_TARGET int64 MWS targeting bit. SCND_TARGET int64 Secondary target targeting bit. - SV1_DESI_TARGET [4]_ int64 SV1 DESI targeting bit. - SV1_BGS_TARGET [4]_ int64 SV1 BGS targeting bit. - SV1_MWS_TARGET [4]_ int64 SV1 MWS targeting bit. - SV2_DESI_TARGET [4]_ int64 SV2 DESI targeting bit. - SV2_BGS_TARGET [4]_ int64 SV2 BGS targeting bit. - SV2_MWS_TARGET [4]_ int64 SV2 MWS targeting bit. - SV3_DESI_TARGET [4]_ int64 SV3 DESI targeting bit. - SV3_BGS_TARGET [4]_ int64 SV3 BGS targeting bit. - SV3_MWS_TARGET [4]_ int64 SV3 MWS targeting bit. - SV1_SCND_TARGET [4]_ int64 SV1 secondary targeting bit. - SV2_SCND_TARGET [4]_ int64 SV2 secondary targeting bit. - SV3_SCND_TARGET [4]_ int64 SV3 secondary targeting bit. + SV1_DESI_TARGET [5]_ int64 SV1 DESI targeting bit. + SV1_BGS_TARGET [5]_ int64 SV1 BGS targeting bit. + SV1_MWS_TARGET [5]_ int64 SV1 MWS targeting bit. + SV2_DESI_TARGET [5]_ int64 SV2 DESI targeting bit. + SV2_BGS_TARGET [5]_ int64 SV2 BGS targeting bit. + SV2_MWS_TARGET [5]_ int64 SV2 MWS targeting bit. + SV3_DESI_TARGET [5]_ int64 SV3 DESI targeting bit. + SV3_BGS_TARGET [5]_ int64 SV3 BGS targeting bit. + SV3_MWS_TARGET [5]_ int64 SV3 MWS targeting bit. + SV1_SCND_TARGET [5]_ int64 SV1 secondary targeting bit. + SV2_SCND_TARGET [5]_ int64 SV2 secondary targeting bit. + SV3_SCND_TARGET [5]_ int64 SV3 secondary targeting bit. Z float64 Redshift based on Redrock or QuasarNet (for QSO targets only). ZWARN int8 Redrock zwarning bit. DELTACHI2 float64 Redrock delta-chi-squared. @@ -1026,4 +1036,11 @@ Data: FITS image [int32, 7781x3,338] .. [3] Column only present when fitting per-exposure tile-based coadds. -.. [4] Column only present in Commissioning and Survey Validation spectroscopic observations. +.. [4] Only observed photometry with a minimum signal-to-noise ratio of two is + used to compute K-corrections. Absolute magnitudes are estimated using + from the (S/N>2) observed-frame bandpass nearest in wavelength to the + desired band-shifted rest-frame bandpass. If no observed-frame photometry + is available, then the absolute magnitude is synthesized from the + best-fitting model and the corresponding inverse variance is set to zero. + +.. [5] Column only present in Commissioning and Survey Validation spectroscopic observations. diff --git a/doc/vacs.rst b/doc/vacs.rst index a9c3a942..017eaa1e 100644 --- a/doc/vacs.rst +++ b/doc/vacs.rst @@ -182,6 +182,7 @@ contributions to the VACs presented herein from the following individuals: * Jeffrey Newman (University of Pittsburgh) * Ragadeepika Pucha (University of Arizona) * Anand Raichoor (Lawrence Berkeley National Lab) +* Dirk Scholte (University College London) * Khaled Said (Australian National University) * David Setton (University of Pittsburgh) * Benjamin Weaver (NSF's NOIRLab) diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py index 9ae994af..ca805680 100644 --- a/py/fastspecfit/continuum.py +++ b/py/fastspecfit/continuum.py @@ -543,8 +543,14 @@ def __init__(self, nophoto=False, fphoto=None, load_filters=True): self.synth_bands = np.array(['g', 'r', 'z']) # for synthesized photometry self.fiber_bands = np.array(['g', 'r', 'z']) # for fiber fluxes - self.absmag_bands = ['U', 'B', 'V', 'sdss_u', 'sdss_g', 'sdss_r', 'sdss_i', 'sdss_z', 'W1', 'W2'] - self.band_shift = [0.0, 0.0, 0.0, 0.1, 0.1, 0.1, 0.1, 0.1, 0.0, 0.0] + self.absmag_bands = ['decam2014_g', 'decam2014_r', 'decam2014_z', + 'U', 'B', 'V', + 'sdss_u', 'sdss_g', 'sdss_r', 'sdss_i', 'sdss_z', + 'W1']#, 'W2'] + self.band_shift = [1.0, 1.0, 1.0, + 0.0, 0.0, 0.0, + 0.1, 0.1, 0.1, 0.1, 0.1, + 0.1]#, 0.1] if load_filters: self.filters = {'N': filters.load_filters('BASS-g', 'BASS-r', 'MzLS-z', @@ -556,11 +562,17 @@ def __init__(self, nophoto=False, fphoto=None, load_filters=True): self.fiber_filters = self.synth_filters self.absmag_filters = filters.FilterSequence(( + filters.load_filter('decam2014-g'), filters.load_filter('decam2014-r'), filters.load_filter('decam2014-z'), filters.load_filter('bessell-U'), filters.load_filter('bessell-B'), filters.load_filter('bessell-V'), filters.load_filter('sdss2010-u'), filters.load_filter('sdss2010-g'), filters.load_filter('sdss2010-r'), filters.load_filter('sdss2010-i'), filters.load_filter('sdss2010-z'), - filters.load_filter('wise2010-W1'), filters.load_filter('wise2010-W2'))) + filters.load_filter('wise2010-W1')))#, filters.load_filter('wise2010-W2'))) + if len(self.absmag_bands) != len(self.band_shift): + errmsg = 'absmag_bands and band_shift must have the same number of elements.' + log.critical(errmsg) + raise ValueError(errmsg) + if self.photounits != 'nanomaggies': errmsg = 'nanomaggies is the only currently supported photometric unit!' log.critical(errmsg) @@ -1748,8 +1760,8 @@ def _kcorr_and_absmag(filters_out, band_shift): dfactor = (1 + redshift) * 4.0 * np.pi * (3.08567758e24 * dlum)**2 / FLUXNORM lums = {} - cwaves = [1500.0, 2800.0, 5100.0] - labels = ['LOGLNU_1500', 'LOGLNU_2800', 'LOGL_5100'] + cwaves = [1500.0, 2800.0, 1450., 1700., 3000., 5100.] + labels = ['LOGLNU_1500', 'LOGLNU_2800', 'LOGL_1450', 'LOGL_1700', 'LOGL_3000', 'LOGL_5100'] norms = [1e28, 1e28, 1e10] for cwave, norm, label in zip(cwaves, norms, labels): J = (templatewave > cwave-500) * (templatewave < cwave+500) @@ -1758,7 +1770,7 @@ def _kcorr_and_absmag(filters_out, band_shift): clipflux, _, _ = sigmaclip(smooth[I], low=1.5, high=3) cflux = np.median(clipflux) # [flux in 10**-17 erg/s/cm2/A] cflux *= dfactor # [monochromatic luminosity in erg/s/A] - if label == 'LOGL_5100': + if 'LOGL_' in label: cflux *= cwave / 3.846e33 / norm # [luminosity in 10**10 L_sun] else: # Convert the UV fluxes to rest-frame luminosity in erg/s/Hz. This @@ -1769,8 +1781,8 @@ def _kcorr_and_absmag(filters_out, band_shift): lums[label] = np.log10(cflux) # * u.erg/(u.second*u.Hz) cfluxes = {} - cwaves = [3728.483, 4862.683, 5008.239, 6564.613] - labels = ['FOII_3727_CONT', 'FHBETA_CONT', 'FOIII_5007_CONT', 'FHALPHA_CONT'] + cwaves = [1215.67, 3728.483, 4862.683, 5008.239, 6564.613] + labels = ['FLYA_1215_CONT', 'FOII_3727_CONT', 'FHBETA_CONT', 'FOIII_5007_CONT', 'FHALPHA_CONT'] for cwave, label in zip(cwaves, labels): J = (templatewave > cwave-500) * (templatewave < cwave+500) I = (templatewave[J] > cwave-20) * (templatewave[J] < cwave+20) @@ -2181,10 +2193,10 @@ def _younger_than_universe(age, tuniv, agepad=0.5): #result['FAGN'] = fagn result['DN4000_MODEL'] = dn4000_model - for iband, band in enumerate(CTools.absmag_bands): - result['KCORR_{}'.format(band.upper())] = kcorr[iband] # * u.mag - result['ABSMAG_{}'.format(band.upper())] = absmag[iband] # * u.mag - result['ABSMAG_IVAR_{}'.format(band.upper())] = ivarabsmag[iband] # / (u.mag**2) + for iband, (band, shift) in enumerate(zip(CTools.absmag_bands, CTools.band_shift)): + result['KCORR{:02d}_{}'.format(int(10*shift), band.upper())] = kcorr[iband] # * u.mag + result['ABSMAG{:02d}_{}'.format(int(10*shift), band.upper())] = absmag[iband] # * u.mag + result['ABSMAG{:02d}_IVAR_{}'.format(int(10*shift), band.upper())] = ivarabsmag[iband] # / (u.mag**2) for iband, band in enumerate(CTools.bands): result['FLUX_SYNTH_PHOTMODEL_{}'.format(band.upper())] = 1e9 * synth_bestmaggies[iband] # * u.nanomaggy if bool(lums): diff --git a/py/fastspecfit/data/legacysurvey-dr9.yaml b/py/fastspecfit/data/legacysurvey-dr9.yaml index b33e651f..ff5780e3 100644 --- a/py/fastspecfit/data/legacysurvey-dr9.yaml +++ b/py/fastspecfit/data/legacysurvey-dr9.yaml @@ -40,10 +40,17 @@ fluxcols: ['FLUX_G', 'FLUX_R', 'FLUX_Z', fluxivarcols: ['FLUX_IVAR_G', 'FLUX_IVAR_R', 'FLUX_IVAR_Z', 'FLUX_IVAR_W1', 'FLUX_IVAR_W2', 'FLUX_IVAR_W3', 'FLUX_IVAR_W4'] -absmag_bands: ['U', 'B', 'V', 'sdss_u', 'sdss_g', 'sdss_r', 'sdss_i', 'sdss_z', 'W1', 'W2'] - -band_shift: [0.0, 0.0, 0.0, 0.1, 0.1, 0.1, 0.1, 0.1, 0.0, 0.0] - -absmag_filters: ['bessell-U', 'bessell-B', 'bessell-V', +absmag_bands: ['decam_g', 'decam_r', 'decam_z', + 'U', 'B', 'V', + 'sdss_u', 'sdss_g', 'sdss_r', 'sdss_i', 'sdss_z', + 'W1'] + +band_shift: [1.0, 1.0, 1.0, + 0.0, 0.0, 0.0, + 0.1, 0.1, 0.1, 0.1, 0.1, + 0.1] + +absmag_filters: ['decam2014-g', 'decam2014-r', 'decam2014-z', + 'bessell-U', 'bessell-B', 'bessell-V', 'sdss2010-u', 'sdss2010-g', 'sdss2010-r', 'sdss2010-i', 'sdss2010-z', - 'wise2010-W1', 'wise2010-W2'] + 'wise2010-W1'] diff --git a/py/fastspecfit/emlines.py b/py/fastspecfit/emlines.py index 2837d2d7..70b91447 100644 --- a/py/fastspecfit/emlines.py +++ b/py/fastspecfit/emlines.py @@ -1096,7 +1096,7 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, """ from scipy.stats import sigmaclip - from fastspecfit.emlines import build_emline_model + from fastspecfit.util import centers2edges for param in finalfit: val = param['value'] @@ -1121,7 +1121,8 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, else: result[param['param_name'].upper()] = val - dpixwave = emlinewave[1]-emlinewave[0] # pixel size [Angstrom] + emlinewave_edges = centers2edges(emlinewave) + dpixwave = np.diff(emlinewave_edges)[0] # pixel size [Angstrom] # zero out all out-of-range lines for oneline in self.fit_linetable[~self.fit_linetable['inrange']]: @@ -1187,10 +1188,11 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, log.critical(errmsg) raise ValueError(errmsg) - # boxcar integration of the flux; should we weight by the line-profile??? - boxflux = np.trapz(emlineflux[lineindx], x=emlinewave[lineindx]) - boxflux_ivar = 1 / np.trapz(1 / emlineivar[lineindx], x=emlinewave[lineindx]) - + # boxcar integration of the flux + dwave = np.median(np.abs(np.diff(emlinewave_edges[lineindx]))) + boxflux = np.sum(emlineflux[lineindx] * dwave) + boxflux_ivar = 1 / np.sum((1 / emlineivar[lineindx]) * dwave**2) + result['{}_BOXFLUX'.format(linename)] = boxflux # * u.erg/(u.second*u.cm**2) result['{}_BOXFLUX_IVAR'.format(linename)] = boxflux_ivar # * u.second**2*u.cm**4/u.erg**2 @@ -1210,21 +1212,22 @@ def populate_emtable(self, result, finalfit, finalmodel, emlinewave, emlineflux, linenorm = np.sqrt(2.0 * np.pi) * linesigma_ang # * u.Angstrom result['{}_FLUX'.format(linename)] = result['{}_MODELAMP'.format(linename)] * linenorm - ## weight by the per-pixel inverse variance line-profile - #lineprofile = build_emline_model(self.dlog10wave, redshift, np.array([result['{}_MODELAMP'.format(linename)]]), - # np.array([result['{}_VSHIFT'.format(linename)]]), np.array([result['{}_SIGMA'.format(linename)]]), - # np.array([oneline['restwave']]), emlinewave, resolution_matrix, camerapix) - # - #weight = np.sum(lineprofile[lineindx])#**2 - #if weight == 0.0: - # errmsg = 'Line-profile should never sum to zero!' - # log.critical(errmsg) - # raise ValueError(errmsg) - # - #flux_ivar = weight / np.sum(lineprofile[lineindx] / emlineivar[lineindx]) - #result['{}_FLUX_IVAR'.format(linename)] = flux_ivar # * u.second**2*u.cm**4/u.erg**2 + # weight by the per-pixel inverse variance line-profile + lineprofile = build_emline_model(self.dlog10wave, redshift, + np.array([result['{}_MODELAMP'.format(linename)]]), + np.array([result['{}_VSHIFT'.format(linename)]]), + np.array([result['{}_SIGMA'.format(linename)]]), + np.array([oneline['restwave']]), emlinewave, + resolution_matrix, camerapix) - result['{}_FLUX_IVAR'.format(linename)] = boxflux_ivar # * u.second**2*u.cm**4/u.erg**2 + lineweight = np.sum(lineprofile[lineindx])**2 + if lineweight == 0.0: + errmsg = 'Line-profile should never sum to zero!' + log.critical(errmsg) + raise ValueError(errmsg) + + flux_ivar = lineweight / np.sum(lineprofile[lineindx]**2 / emlineivar[lineindx]) + result['{}_FLUX_IVAR'.format(linename)] = flux_ivar # * u.second**2*u.cm**4/u.erg**2 result['{}_CHI2'.format(linename)] = np.sum(emlineivar[lineindx] * (emlineflux[lineindx] - finalmodel[lineindx])**2) diff --git a/py/fastspecfit/fastspecfit.py b/py/fastspecfit/fastspecfit.py index f0b253ab..758995ef 100644 --- a/py/fastspecfit/fastspecfit.py +++ b/py/fastspecfit/fastspecfit.py @@ -15,10 +15,6 @@ def _fastspec_one(args): """Multiprocessing wrapper.""" return fastspec_one(*args) -def _desiqa_one(args): - """Multiprocessing wrapper.""" - return desiqa_one(*args) - def _assign_units_to_columns(fastfit, metadata, Spec, templates, fastphot, stackfit, log): """Assign astropy units to output tables. @@ -78,32 +74,6 @@ def fastspec_one(iobj, data, out, meta, fphoto, templates, log=None, return out, meta, emmodel -def desiqa_one(data, fastfit, metadata, templates, coadd_type, fphoto, - minspecwave=3500., maxspecwave=9900., minphotwave=0.1, - maxphotwave=35., emline_snrmin=0.0, nsmoothspec=1, - fastphot=False, nophoto=False, stackfit=False, inputz=False, - no_smooth_continuum=False, outdir=None, outprefix=None, log=None): - """Multiprocessing wrapper to generate QA for a single object. - - """ - from fastspecfit.io import cache_templates - templatecache = cache_templates(templates, log=log, mintemplatewave=450.0, - maxtemplatewave=40e4, fastphot=fastphot) - - if inputz: - from fastspecfit.util import TabulatedDESI - cosmo = TabulatedDESI() - else: - cosmo = None - - qa_fastspec(data, templatecache, fastfit, metadata, coadd_type=coadd_type, - spec_wavelims=(minspecwave, maxspecwave), - phot_wavelims=(minphotwave, maxphotwave), - no_smooth_continuum=no_smooth_continuum, - emline_snrmin=emline_snrmin, nsmoothspec=nsmoothspec, - fastphot=fastphot, fphoto=fphoto, stackfit=stackfit, - outprefix=outprefix, outdir=outdir, log=log, cosmo=cosmo) - def parse(options=None, log=None): """Parse input arguments to fastspec and fastphot scripts. @@ -112,7 +82,7 @@ def parse(options=None, log=None): parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('redrockfiles', nargs='*', help='Full path to input redrock file(s).') + parser.add_argument('redrockfiles', nargs='+', help='Full path to input redrock file(s).') parser.add_argument('-o', '--outfile', type=str, required=True, help='Full path to output filename (required).') parser.add_argument('--mp', type=int, default=1, help='Number of multiprocessing threads per MPI rank.') parser.add_argument('-n', '--ntargets', type=int, help='Number of targets to process in each file.') @@ -304,1124 +274,3 @@ def stackfit(args=None, comm=None): """ fastspec(stackfit=True, args=args, comm=comm) - -def qa_fastspec(data, templatecache, fastspec, metadata, coadd_type='healpix', - spec_wavelims=(3550, 9900), phot_wavelims=(0.1, 35), - fastphot=False, fphoto=None, stackfit=False, outprefix=None, - no_smooth_continuum=False, emline_snrmin=0.0, nsmoothspec=1, - outdir=None, log=None, cosmo=None): - """QA plot the emission-line spectrum and best-fitting model. - - """ - import subprocess - from scipy.ndimage import median_filter - - import matplotlib.pyplot as plt - import matplotlib.ticker as ticker - from matplotlib import colors - from matplotlib.patches import Circle, Rectangle, ConnectionPatch - from matplotlib.lines import Line2D - import matplotlib.gridspec as gridspec - import matplotlib.image as mpimg - - import astropy.units as u - from astropy.io import fits - from astropy.wcs import WCS - import seaborn as sns - from PIL import Image, ImageDraw - - from fastspecfit.util import ivar2var, C_LIGHT - from fastspecfit.io import FLUXNORM, get_qa_filename - from fastspecfit.continuum import ContinuumTools - from fastspecfit.emlines import build_emline_model, EMFitTools - - if log is None: - from desiutil.log import get_logger - log = get_logger() - - Image.MAX_IMAGE_PIXELS = None - - sns.set(context='talk', style='ticks', font_scale=1.3)#, rc=rc) - - if stackfit: - col1 = [colors.to_hex(col) for col in ['violet']] - col2 = [colors.to_hex(col) for col in ['purple']] - else: - col1 = [colors.to_hex(col) for col in ['dodgerblue', 'darkseagreen', 'orangered']] - col2 = [colors.to_hex(col) for col in ['darkblue', 'darkgreen', 'darkred']] - - photcol1 = colors.to_hex('darkorange') - - @ticker.FuncFormatter - def major_formatter(x, pos): - if (x >= 0.01) and (x < 0.1): - return f'{x:.2f}' - elif (x >= 0.1) and (x < 1): - return f'{x:.1f}' - else: - return f'{x:.0f}' - - CTools = ContinuumTools(fphoto=fphoto, - continuum_pixkms=templatecache['continuum_pixkms'], - pixkms_wavesplit=templatecache['pixkms_wavesplit']) - if 'legacysurveydr' in fphoto.keys(): - layer = 'ls-{}'.format(fphoto['legacysurveydr']) - else: - layer = 'ls-dr9' - - if not fastphot: - EMFit = EMFitTools(minspecwave=spec_wavelims[0], maxspecwave=spec_wavelims[1]) - - filters = CTools.synth_filters[metadata['PHOTSYS']] - allfilters = CTools.filters[metadata['PHOTSYS']] - redshift = fastspec['Z'] - - # cosmo will be provided if input_redshifts are used - if cosmo is not None: - dlum = cosmo.luminosity_distance(redshift) - else: - dlum = data['dluminosity'] - - pngfile = get_qa_filename(metadata, coadd_type, outprefix=outprefix, - outdir=outdir, fastphot=fastphot, log=log) - - # some arrays to use for the legend - if coadd_type == 'healpix': - target = [ - 'Survey/Program/Healpix: {}/{}/{}'.format(metadata['SURVEY'], metadata['PROGRAM'], metadata['HEALPIX']), - 'TargetID: {}'.format(metadata['TARGETID']), - ] - elif coadd_type == 'cumulative': - target = [ - 'Tile/Night/Fiber: {}/{}/{}'.format(metadata['TILEID'], metadata['NIGHT'], metadata['FIBER']), - 'TargetID: {}'.format(metadata['TARGETID']), - ] - elif coadd_type == 'pernight': - target = [ - 'Tile/Night/Fiber: {}/{}/{}'.format(metadata['TILEID'], metadata['NIGHT'], metadata['FIBER']), - 'TargetID: {}'.format(metadata['TARGETID']), - ] - elif coadd_type == 'perexp': - target = [ - 'Tile/Night/Fiber: {}/{}/{}'.format(metadata['TILEID'], metadata['NIGHT'], metadata['FIBER']), - 'Night/Fiber: {}/{}'.format(metadata['NIGHT'], metadata['FIBER']), - 'TargetID: {}'.format(metadata['TARGETID']), - ] - elif coadd_type == 'custom': - target = [ - 'Survey/Program/Healpix: {}/{}/{}'.format(metadata['SURVEY'], metadata['PROGRAM'], metadata['HEALPIX']), - 'TargetID: {}'.format(metadata['TARGETID']), - ] - elif coadd_type == 'stacked': - target = [ - 'StackID: {}'.format(metadata['STACKID']), - '', - ] - else: - errmsg = 'Unrecognized coadd_type {}!'.format(coadd_type) - log.critical(errmsg) - raise ValueError(errmsg) - - leg = { - 'z': '$z={:.7f}$'.format(redshift), - 'rchi2_phot': r'$\chi^{2}_{\nu,\mathrm{phot}}=$'+r'${:.2f}$'.format(fastspec['RCHI2_PHOT']), - 'dn4000_model': r'$D_{n}(4000)_{\mathrm{model}}=$'+r'${:.3f}$'.format(fastspec['DN4000_MODEL']), - 'age': r'Age$={:.3f}$ Gyr'.format(fastspec['AGE']), - 'AV': r'$A_{V}=$'+r'${:.3f}$ mag'.format(fastspec['AV']), - 'mstar': r'$\log_{10}(M/M_{\odot})=$'+r'${:.3f}$'.format(fastspec['LOGMSTAR']), - 'sfr': r'$\mathrm{SFR}=$'+'${:.1f}$'.format(fastspec['SFR'])+' $M_{\odot}/\mathrm{yr}$', - 'zzsun': r'$Z/Z_{\odot}=$'+r'${:.3f}$'.format(fastspec['ZZSUN']), - } - - # try to figure out which absmags to display - gindx = np.argmin(np.abs(CTools.absmag_filters.effective_wavelengths.value - 4900)) - rindx = np.argmin(np.abs(CTools.absmag_filters.effective_wavelengths.value - 6500)) - zindx = np.argmin(np.abs(CTools.absmag_filters.effective_wavelengths.value - 9200)) - absmag_gband = CTools.absmag_bands[gindx] - absmag_rband = CTools.absmag_bands[rindx] - absmag_zband = CTools.absmag_bands[zindx] - - leg.update({'absmag_r': '$M_{{{}}}={:.2f}$'.format(absmag_rband.lower().replace('decam_', '').replace('sdss_', ''), - fastspec['ABSMAG_{}'.format(absmag_rband.upper())])}) - if gindx != rindx: - gr = fastspec['ABSMAG_{}'.format(absmag_gband.upper())] - fastspec['ABSMAG_{}'.format(absmag_rband.upper())] - leg.update({'absmag_gr': '$M_{{{}}}-M_{{{}}}={:.2f}$'.format(absmag_gband.lower(), absmag_rband.lower(), gr).replace('decam_', '').replace('sdss_', '')}) - if zindx != rindx: - rz = fastspec['ABSMAG_{}'.format(absmag_rband.upper())] - fastspec['ABSMAG_{}'.format(absmag_zband.upper())] - leg.update({'absmag_rz': '$M_{{{}}}-M_{{{}}}={:.2f}$'.format(absmag_rband.lower(), absmag_zband.lower(), rz).replace('decam_', '').replace('sdss_', '')}) - - #leg.update({ - # 'absmag_r': '$M_{{0.1r}}={:.2f}$'.format(fastspec['ABSMAG_SDSS_R']), - # 'absmag_gr': '$^{{0.1}}(g-r)={:.3f}$'.format(fastspec['ABSMAG_SDSS_G']-fastspec['ABSMAG_SDSS_R']), - # 'absmag_rz': '$^{{0.1}}(r-z)={:.3f}$'.format(fastspec['ABSMAG_SDSS_R']-fastspec['ABSMAG_SDSS_Z']), - # }) - - #leg['radec'] = '$(\\alpha,\\delta)=({:.7f},{:.6f})$'.format(metadata['RA'], metadata['DEC']) - #leg['zwarn'] = '$z_{{\\rm warn}}={}$'.format(metadata['ZWARN']) - - if fastphot: - leg['vdisp'] = r'$\sigma_{star}=$'+'{:g}'.format(fastspec['VDISP'])+' km/s' - else: - if fastspec['VDISP_IVAR'] > 0: - leg['vdisp'] = r'$\sigma_{{star}}={:.0f}\pm{:.0f}$ km/s'.format(fastspec['VDISP'], 1/np.sqrt(fastspec['VDISP_IVAR'])) - else: - leg['vdisp'] = r'$\sigma_{{star}}={:g}$ km/s'.format(fastspec['VDISP']) - - leg['rchi2'] = r'$\chi^{2}_{\nu,\mathrm{specphot}}$='+'{:.2f}'.format(fastspec['RCHI2']) - leg['rchi2_cont'] = r'$\chi^{2}_{\nu,\mathrm{cont}}$='+'{:.2f}'.format(fastspec['RCHI2_CONT']) - - if not stackfit: - if redshift != metadata['Z_RR']: - leg['zredrock'] = '$z_{\mathrm{Redrock}}=$'+r'${:.7f}$'.format(metadata['Z_RR']) - - if fastphot: - fontsize1 = 16 - fontsize2 = 22 - else: - if stackfit: - fontsize1 = 16 - fontsize2 = 22 - else: - fontsize1 = 18 # 24 - fontsize2 = 24 - - apercorr = fastspec['APERCORR'] - - if fastspec['DN4000_IVAR'] > 0: - leg['dn4000_spec'] = r'$D_{n}(4000)_{\mathrm{data}}=$'+r'${:.3f}$'.format(fastspec['DN4000']) - - # kinematics - if fastspec['NARROW_Z'] != redshift: - if fastspec['NARROW_ZRMS'] > 0: - leg['dv_narrow'] = r'$\Delta v_{\mathrm{narrow}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( - C_LIGHT*(fastspec['NARROW_Z']-redshift), C_LIGHT*fastspec['NARROW_ZRMS']) - else: - leg['dv_narrow'] = r'$\Delta v_{\mathrm{narrow}}=$'+r'${:.0f}$ km/s'.format( - C_LIGHT*(fastspec['NARROW_Z']-redshift)) - if fastspec['NARROW_SIGMA'] != 0.0: - if fastspec['NARROW_SIGMARMS'] > 0: - leg['sigma_narrow'] = r'$\sigma_{\mathrm{narrow}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( - fastspec['NARROW_SIGMA'], fastspec['NARROW_SIGMARMS']) - else: - leg['sigma_narrow'] = r'$\sigma_{\mathrm{narrow}}=$'+r'${:.0f}$ km/s'.format(fastspec['NARROW_SIGMA']) - - snrcut = 1.5 - leg_broad, leg_narrow, leg_uv = {}, {}, {} - - if fastspec['UV_Z'] != redshift: - if fastspec['UV_ZRMS'] > 0: - leg_uv['dv_uv'] = r'$\Delta v_{\mathrm{UV}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( - C_LIGHT*(fastspec['UV_Z']-redshift), C_LIGHT*fastspec['UV_ZRMS']) - else: - leg_uv['dv_uv'] = r'$\Delta v_{\mathrm{UV}}=$'+r'${:.0f}$ km/s'.format( - C_LIGHT*(fastspec['UV_Z']-redshift)) - if fastspec['UV_SIGMA'] != 0.0: - if fastspec['UV_SIGMARMS'] > 0: - leg_uv['sigma_uv'] = r'$\sigma_{\mathrm{UV}}$'+r'$={:.0f}\pm{:.0f}$ km/s'.format( - fastspec['UV_SIGMA'], fastspec['UV_SIGMARMS']) - else: - leg_uv['sigma_uv'] = r'$\sigma_{\mathrm{UV}}=$'+r'${:.0f}$ km/s'.format(fastspec['UV_SIGMA']) - if fastspec['BROAD_Z'] != redshift: - if fastspec['BROAD_ZRMS'] > 0: - leg_broad['dv_broad'] = r'$\Delta v_{\mathrm{broad}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( - C_LIGHT*(fastspec['BROAD_Z']-redshift), C_LIGHT*fastspec['BROAD_ZRMS']) - else: - leg_broad['dv_broad'] = r'$\Delta v_{\mathrm{broad}}=$'+r'${:.0f}$ km/s'.format( - C_LIGHT*(fastspec['BROAD_Z']-redshift)) - if fastspec['BROAD_SIGMA'] != 0.0: - if fastspec['BROAD_SIGMARMS'] > 0: - leg_broad['sigma_broad'] = r'$\sigma_{\mathrm{broad}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( - fastspec['BROAD_SIGMA'], fastspec['BROAD_SIGMARMS']) - else: - leg_broad['sigma_broad'] = r'$\sigma_{\mathrm{broad}}=$'+r'${:.0f}$ km/s'.format(fastspec['BROAD_SIGMA']) - - # emission lines - - # UV - if fastspec['LYALPHA_AMP']*np.sqrt(fastspec['LYALPHA_AMP_IVAR']) > snrcut: - leg_uv['ewlya'] = r'EW(Ly$\alpha)=$'+r'${:.1f}$'.format(fastspec['LYALPHA_EW'])+r' $\AA$' - if fastspec['CIV_1549_AMP']*np.sqrt(fastspec['CIV_1549_AMP_IVAR']) > snrcut: - leg_uv['ewciv'] = r'EW(CIV)=$'+r'${:.1f}$'.format(fastspec['CIV_1549_EW'])+r' $\AA$' - if fastspec['CIII_1908_AMP']*np.sqrt(fastspec['CIII_1908_AMP_IVAR']) > snrcut: - leg_uv['ewciii'] = r'EW(CIII])=$'+r'${:.1f}$'.format(fastspec['CIII_1908_EW'])+r' $\AA$' - if (fastspec['MGII_2796_AMP']*np.sqrt(fastspec['MGII_2796_AMP_IVAR']) > snrcut or - fastspec['MGII_2803_AMP']*np.sqrt(fastspec['MGII_2803_AMP_IVAR']) > snrcut): - leg_uv['ewmgii'] = r'EW(MgII)=$'+r'${:.1f}$'.format(fastspec['MGII_2796_EW']+fastspec['MGII_2803_EW'])+r' $\AA$' - leg_uv['mgii_doublet'] = r'MgII $\lambda2796/\lambda2803={:.3f}$'.format(fastspec['MGII_DOUBLET_RATIO']) - - leg_broad['linerchi2'] = r'$\chi^{2}_{\nu,\mathrm{line}}=$'+r'${:.2f}$'.format(fastspec['RCHI2_LINE']) - leg_broad['deltachi2'] = r'$\Delta\chi^{2}_{\mathrm{nobroad}}=$'+r'${:.2f}$'.format(fastspec['DELTA_LINECHI2']) - leg_broad['deltandof'] = r'$\Delta\nu_{\mathrm{nobroad}}=$'+r'${:.0f}$'.format(fastspec['DELTA_LINENDOF']) - - # choose one broad Balmer line - if fastspec['HALPHA_BROAD_AMP']*np.sqrt(fastspec['HALPHA_BROAD_AMP_IVAR']) > snrcut: - leg_broad['ewbalmer_broad'] = r'EW(H$\alpha)_{\mathrm{broad}}=$'+r'${:.1f}$'.format(fastspec['HALPHA_BROAD_EW'])+r' $\AA$' - elif fastspec['HBETA_BROAD_AMP']*np.sqrt(fastspec['HBETA_BROAD_AMP_IVAR']) > snrcut: - leg_broad['ewbalmer_broad'] = r'EW(H$\beta)_{\mathrm{broad}}=$'+r'${:.1f}$'.format(fastspec['HBETA_BROAD_EW'])+r' $\AA$' - elif fastspec['HGAMMA_BROAD_AMP']*np.sqrt(fastspec['HGAMMA_BROAD_AMP_IVAR']) > snrcut: - leg_broad['ewbalmer_broad'] = r'EW(H$\gamma)_{\mathrm{broad}}=$'+r'${:.1f}$'.format(fastspec['HGAMMA_BROAD_EW'])+r' $\AA$' - - if (fastspec['OII_3726_AMP']*np.sqrt(fastspec['OII_3726_AMP_IVAR']) > snrcut or - fastspec['OII_3729_AMP']*np.sqrt(fastspec['OII_3729_AMP_IVAR']) > snrcut): - leg_narrow['ewoii'] = r'EW([OII])'+r'$={:.1f}$'.format(fastspec['OII_3726_EW']+fastspec['OII_3729_EW'])+r' $\AA$' - - if fastspec['OIII_5007_AMP']*np.sqrt(fastspec['OIII_5007_AMP_IVAR']) > snrcut: - leg_narrow['ewoiii'] = r'EW([OIII])'+r'$={:.1f}$'.format(fastspec['OIII_5007_EW'])+r' $\AA$' - - # choose one Balmer line - if fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut: - leg_narrow['ewbalmer_narrow'] = r'EW(H$\alpha)=$'+r'${:.1f}$'.format(fastspec['HALPHA_EW'])+r' $\AA$' - elif fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut: - leg_narrow['ewbalmer_narrow'] = r'EW(H$\beta)=$'+r'${:.1f}$'.format(fastspec['HBETA_EW'])+r' $\AA$' - elif fastspec['HGAMMA_AMP']*np.sqrt(fastspec['HGAMMA_AMP_IVAR']) > snrcut: - leg_narrow['ewbalmer_narrow'] = r'EW(H$\gamma)=$'+r'${:.1f}$'.format(fastspec['HGAMMA_EW'])+r' $\AA$' - - if (fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut and - fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut): - leg_narrow['hahb'] = r'$\mathrm{H}\alpha/\mathrm{H}\beta=$'+r'${:.3f}$'.format(fastspec['HALPHA_FLUX']/fastspec['HBETA_FLUX']) - if 'hahb' not in leg_narrow.keys() and (fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut and - fastspec['HGAMMA_AMP']*np.sqrt(fastspec['HGAMMA_AMP_IVAR']) > snrcut): - leg_narrow['hbhg'] = r'$\mathrm{H}\beta/\mathrm{H}\gamma=$'+r'${:.3f}$'.format(fastspec['HBETA_FLUX']/fastspec['HGAMMA_FLUX']) - if (fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut and - fastspec['OIII_5007_AMP']*np.sqrt(fastspec['OIII_5007_AMP_IVAR']) > snrcut and - fastspec['HBETA_FLUX'] > 0 and fastspec['OIII_5007_FLUX'] > 0): - leg_narrow['oiiihb'] = r'$\log_{10}(\mathrm{[OIII]/H}\beta)=$'+r'${:.3f}$'.format(np.log10(fastspec['OIII_5007_FLUX']/fastspec['HBETA_FLUX'])) - if (fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut and - fastspec['NII_6584_AMP']*np.sqrt(fastspec['NII_6584_AMP_IVAR']) > snrcut and - fastspec['HALPHA_FLUX'] > 0 and fastspec['NII_6584_FLUX'] > 0): - leg_narrow['niiha'] = r'$\log_{10}(\mathrm{[NII]/H}}\alpha)=$'+r'${:.3f}$'.format(np.log10(fastspec['NII_6584_FLUX']/fastspec['HALPHA_FLUX'])) - - if (fastspec['OII_3726_AMP']*np.sqrt(fastspec['OII_3726_AMP_IVAR']) > snrcut or - fastspec['OII_3729_AMP']*np.sqrt(fastspec['OII_3729_AMP_IVAR']) > snrcut): - #if fastspec['OII_DOUBLET_RATIO'] != 0: - leg_narrow['oii_doublet'] = r'[OII] $\lambda3726/\lambda3729={:.3f}$'.format(fastspec['OII_DOUBLET_RATIO']) - - if fastspec['SII_6716_AMP']*np.sqrt(fastspec['SII_6716_AMP_IVAR']) > snrcut or fastspec['SII_6731_AMP']*np.sqrt(fastspec['SII_6731_AMP_IVAR']) > snrcut: - #if fastspec['SII_DOUBLET_RATIO'] != 0: - leg_narrow['sii_doublet'] = r'[SII] $\lambda6731/\lambda6716={:.3f}$'.format(fastspec['SII_DOUBLET_RATIO']) - - # rebuild the best-fitting broadband photometric model - if not stackfit: - sedmodel, sedphot = CTools.templates2data( - templatecache['templateflux'], templatecache['templatewave'], - redshift=redshift, dluminosity=dlum, photsys=metadata['PHOTSYS'], - synthphot=True, coeff=fastspec['COEFF'] * CTools.massnorm) - sedwave = templatecache['templatewave'] * (1 + redshift) - - phot = CTools.parse_photometry(CTools.bands, - maggies=np.array([metadata['FLUX_{}'.format(band.upper())] for band in CTools.bands]), - ivarmaggies=np.array([metadata['FLUX_IVAR_{}'.format(band.upper())] for band in CTools.bands]), - lambda_eff=allfilters.effective_wavelengths.value, - min_uncertainty=CTools.min_uncertainty) - #if hasattr(CTools, 'fiber_bands'): - # fiberphot = CTools.parse_photometry(CTools.fiber_bands, - # maggies=np.array([metadata['FIBERTOTFLUX_{}'.format(band.upper())] - # for band in CTools.fiber_bands]), - # lambda_eff=filters.effective_wavelengths.value) - #else: - # fiberphot = None - - indx_phot = np.where((sedmodel > 0) * (sedwave/1e4 > phot_wavelims[0]) * - (sedwave/1e4 < phot_wavelims[1]))[0] - sedwave = sedwave[indx_phot] - sedmodel = sedmodel[indx_phot] - - if not fastphot: - # Rebuild the best-fitting spectroscopic model; prefix "desi" means - # "per-camera" and prefix "full" has the cameras h-stacked. - fullwave = np.hstack(data['wave']) - - desicontinuum, _ = CTools.templates2data(templatecache['templateflux_nolines'], templatecache['templatewave'], - redshift=redshift, dluminosity=dlum, synthphot=False, - specwave=data['wave'], specres=data['res'], - specmask=data['mask'], cameras=data['cameras'], - vdisp=fastspec['VDISP'], - coeff=fastspec['COEFF']) - - # remove the aperture correction - desicontinuum = [_desicontinuum / apercorr for _desicontinuum in desicontinuum] - fullcontinuum = np.hstack(desicontinuum) - - # Need to be careful we don't pass a large negative residual where - # there are gaps in the data. - desiresiduals = [] - for icam in np.arange(len(data['cameras'])): - resid = data['flux'][icam] - desicontinuum[icam] - I = (data['flux'][icam] == 0.0) * (data['flux'][icam] == 0.0) - if np.any(I): - resid[I] = 0.0 - desiresiduals.append(resid) - - if np.all(fastspec['COEFF'] == 0): - fullsmoothcontinuum = np.zeros_like(fullwave) - else: - fullsmoothcontinuum, _ = CTools.smooth_continuum( - fullwave, np.hstack(desiresiduals), np.hstack(data['ivar']), - redshift=redshift, linemask=np.hstack(data['linemask'])) - if no_smooth_continuum: - fullsmoothcontinuum *= 0 - - desismoothcontinuum = [] - for campix in data['camerapix']: - desismoothcontinuum.append(fullsmoothcontinuum[campix[0]:campix[1]]) - - # full model spectrum + individual line-spectra - desiemlines = EMFit.emlinemodel_bestfit(data['wave'], data['res'], fastspec, snrcut=emline_snrmin) - - desiemlines_oneline = [] - inrange = ( (EMFit.linetable['restwave'] * (1+redshift) > np.min(fullwave)) * - (EMFit.linetable['restwave'] * (1+redshift) < np.max(fullwave)) ) - - for oneline in EMFit.linetable[inrange]: # for all lines in range - linename = oneline['name'].upper() - amp = fastspec['{}_MODELAMP'.format(linename)] - if amp != 0: - desiemlines_oneline1 = build_emline_model( - EMFit.dlog10wave, redshift, np.array([amp]), - np.array([fastspec['{}_VSHIFT'.format(linename)]]), - np.array([fastspec['{}_SIGMA'.format(linename)]]), - np.array([oneline['restwave']]), data['wave'], data['res']) - desiemlines_oneline.append(desiemlines_oneline1) - - # Grab the viewer cutout. - if not stackfit: - pixscale = 0.262 - width = int(30 / pixscale) # =1 arcmin - height = int(width / 1.3) # 3:2 aspect ratio - - hdr = fits.Header() - hdr['NAXIS'] = 2 - hdr['NAXIS1'] = width - hdr['NAXIS2'] = height - hdr['CTYPE1'] = 'RA---TAN' - hdr['CTYPE2'] = 'DEC--TAN' - hdr['CRVAL1'] = metadata['RA'] - hdr['CRVAL2'] = metadata['DEC'] - hdr['CRPIX1'] = width/2+0.5 - hdr['CRPIX2'] = height/2+0.5 - hdr['CD1_1'] = -pixscale/3600 - hdr['CD1_2'] = 0.0 - hdr['CD2_1'] = 0.0 - hdr['CD2_2'] = +pixscale/3600 - wcs = WCS(hdr) - - cutoutjpeg = os.path.join(outdir, 'tmp.'+os.path.basename(pngfile.replace('.png', '.jpeg'))) - if not os.path.isfile(cutoutjpeg): - wait = 5 # seconds - #cmd = 'timeout {wait} wget -q -o /dev/null -O {outfile} https://www.legacysurvey.org/viewer/jpeg-cutout?ra={ra}&dec={dec}&width={width}&height={height}&layer={layer}' - #cmd = cmd.format(wait=wait, outfile=cutoutjpeg, ra=metadata['RA'], - # dec=metadata['DEC'], width=width, height=height, - # layer=layer) - cmd = 'wget -q -O {outfile} https://www.legacysurvey.org/viewer/jpeg-cutout?ra={ra}&dec={dec}&width={width}&height={height}&layer={layer}' - cmd = cmd.format(outfile=cutoutjpeg, ra=metadata['RA'], dec=metadata['DEC'], - width=width, height=height, layer=layer) - log.info(cmd) - try: - subprocess.check_output(cmd.split(), stderr=subprocess.DEVNULL, timeout=wait) - #subprocess.check_call(cmd.split(), stderr=subprocess.DEVNULL) - except: - log.warning('No cutout from viewer after {} seconds; stopping wget'.format(wait)) - try: - img = mpimg.imread(cutoutjpeg) - except: - log.warning('Problem reading cutout for targetid {}.'.format(metadata['TARGETID'])) - img = np.zeros((height, width, 3)) - - if os.path.isfile(cutoutjpeg): - os.remove(cutoutjpeg) - - # QA choices - legxpos, legypos, legypos2, legfntsz1, legfntsz = 0.98, 0.94, 0.05, 16, 18 - bbox = dict(boxstyle='round', facecolor='lightgray', alpha=0.15) - bbox2 = dict(boxstyle='round', facecolor='lightgray', alpha=0.7) - - if fastphot: - fullheight = 9 # inches - fullwidth = 18 - - nrows = 3 - ncols = 8 - - fig = plt.figure(figsize=(fullwidth, fullheight)) - gs = fig.add_gridspec(nrows, ncols)#, width_ratios=width_ratios) - - cutax = fig.add_subplot(gs[0:2, 5:8], projection=wcs) # rows x cols - sedax = fig.add_subplot(gs[0:3, 0:5]) - elif stackfit: - fullheight = 14 # inches - fullwidth = 24 - - # 8 columns: 3 for the spectra, and 5 for the lines - # 8 rows: 4 for the SED, 2 each for the spectra, 1 gap, and 3 for the lines - nlinerows = 6 - nlinecols = 4 - nrows = nlinerows - ncols = 9 - - #height_ratios = np.hstack(([1.0]*3, 0.25, [1.0]*6)) - #width_ratios = np.hstack(([1.0]*5, [1.0]*3)) - - fig = plt.figure(figsize=(fullwidth, fullheight)) - gs = fig.add_gridspec(nrows, ncols)#, height_ratios=height_ratios)#, width_ratios=width_ratios) - - specax = fig.add_subplot(gs[0:4, 0:5]) - else: - fullheight = 18 # inches - fullwidth = 24 - - # 8 columns: 3 for the SED, 5 for the spectra, and 8 for the lines - # 8 rows: 4 for the SED, 2 each for the spectra, 1 gap, and 3 for the lines - ngaprows = 1 - nlinerows = 6 - nlinecols = 3 - nrows = 9 + ngaprows - ncols = 8 - - height_ratios = np.hstack(([1.0]*3, 0.25, [1.0]*6)) # small gap - #width_ratios = np.hstack(([1.0]*5, [1.0]*3)) - - fig = plt.figure(figsize=(fullwidth, fullheight)) - gs = fig.add_gridspec(nrows, ncols, height_ratios=height_ratios)#, width_ratios=width_ratios) - - cutax = fig.add_subplot(gs[0:3, 5:8], projection=wcs) # rows x cols - sedax = fig.add_subplot(gs[0:3, 0:5]) - specax = fig.add_subplot(gs[4:8, 0:5]) - - # viewer cutout - if not stackfit: - cutax.imshow(img, origin='lower')#, interpolation='nearest') - cutax.set_xlabel('RA [J2000]') - cutax.set_ylabel('Dec [J2000]') - - cutax.coords[1].set_ticks_position('r') - cutax.coords[1].set_ticklabel_position('r') - cutax.coords[1].set_axislabel_position('r') - - if metadata['DEC'] > 0: - sgn = '+' - else: - sgn = '' - - cutax.text(0.04, 0.95, '$(\\alpha,\\delta)$=({:.7f}, {}{:.6f})'.format(metadata['RA'], sgn, metadata['DEC']), - ha='left', va='top', color='k', fontsize=fontsize1, bbox=bbox2, - transform=cutax.transAxes) - - sz = img.shape - cutax.add_artist(Circle((sz[1] / 2, sz[0] / 2), radius=1.5/2/pixscale, facecolor='none', # DESI fiber=1.5 arcsec diameter - edgecolor='red', ls='-', alpha=0.8))#, label='3" diameter')) - cutax.add_artist(Circle((sz[1] / 2, sz[0] / 2), radius=10/2/pixscale, facecolor='none', - edgecolor='red', ls='--', alpha=0.8))#, label='15" diameter')) - handles = [Line2D([0], [0], color='red', lw=2, ls='-', label='1.5 arcsec'), - Line2D([0], [0], color='red', lw=2, ls='--', label='10 arcsec')] - - cutax.legend(handles=handles, loc='lower left', fontsize=fontsize1, facecolor='lightgray') - - if not fastphot: - # plot the full spectrum + best-fitting (total) model - specax.plot(fullwave/1e4, fullsmoothcontinuum, color='gray', alpha=0.4) - specax.plot(fullwave/1e4, fullcontinuum, color='k', alpha=0.6) - - spec_ymin, spec_ymax = 1e6, -1e6 - - desimodelspec = [] - for icam in np.arange(len(data['cameras'])): # iterate over cameras - wave = data['wave'][icam] - flux = data['flux'][icam] - modelflux = desiemlines[icam] + desicontinuum[icam] + desismoothcontinuum[icam] - - sigma, camgood = ivar2var(data['ivar'][icam], sigma=True, allmasked_ok=True, clip=0) - - wave = wave[camgood] - flux = flux[camgood] - sigma = sigma[camgood] - modelflux = modelflux[camgood] - - desimodelspec.append(apercorr * (desicontinuum[icam] + desiemlines[icam])) - - # get the robust range - filtflux = median_filter(flux, 51, mode='nearest') - if np.sum(camgood) > 0: - sigflux = np.diff(np.percentile(flux - modelflux, [25, 75]))[0] / 1.349 # robust - if -2 * sigflux < spec_ymin: - spec_ymin = -2 * sigflux - if 6 * sigflux > spec_ymax: - spec_ymax = 6 * sigflux - if np.max(filtflux) > spec_ymax: - #print(icam, spec_ymax, np.max(filtflux), np.max(filtflux) * 1.2) - spec_ymax = np.max(filtflux) * 1.25 - if np.max(modelflux) > spec_ymax: - spec_ymax = np.max(modelflux) * 1.25 - #print(spec_ymin, spec_ymax) - #pdb.set_trace() - - #specax.fill_between(wave, flux-sigma, flux+sigma, color=col1[icam], alpha=0.2) - if nsmoothspec > 1: - from scipy.ndimage import gaussian_filter - specax.plot(wave/1e4, gaussian_filter(flux, nsmoothspec), color=col1[icam], alpha=0.8) - specax.plot(wave/1e4, gaussian_filter(modelflux, nsmoothspec), color=col2[icam], lw=2, alpha=0.8) - else: - specax.plot(wave/1e4, flux, color=col1[icam], alpha=0.8) - specax.plot(wave/1e4, modelflux, color=col2[icam], lw=2, alpha=0.8) - - fullmodelspec = np.hstack(desimodelspec) - - if stackfit: - specax_twin = specax.twiny() - specax_twin.set_xlim(spec_wavelims[0]/(1+redshift)/1e4, spec_wavelims[1]/(1+redshift)/1e4) - specax_twin.xaxis.set_major_formatter(major_formatter) - restticks = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1]) - restticks = restticks[(restticks >= spec_wavelims[0]/(1+redshift)/1e4) * (restticks <= spec_wavelims[1]/(1+redshift)/1e4)] - specax_twin.set_xticks(restticks) - else: - specax.spines[['top']].set_visible(False) - - specax.set_xlim(spec_wavelims[0]/1e4, spec_wavelims[1]/1e4) - specax.set_ylim(spec_ymin, spec_ymax) - specax.set_xlabel(r'Observed-frame Wavelength ($\mu$m)') - #specax.set_xlabel(r'Observed-frame Wavelength ($\AA$)') - specax.set_ylabel(r'$F_{\lambda}\ (10^{-17}~{\rm erg}~{\rm s}^{-1}~{\rm cm}^{-2}~\AA^{-1})$') - - # photometric SED - if not stackfit: - abmag_good = phot['abmag_ivar'] > 0 - abmag_goodlim = phot['abmag_limit'] > 0 - - if len(sedmodel) == 0: - log.warning('Best-fitting photometric continuum is all zeros or negative!') - if np.sum(abmag_good) > 0: - medmag = np.median(phot['abmag'][abmag_good]) - elif np.sum(abmag_goodlim) > 0: - medmag = np.median(phot['abmag_limit'][abmag_goodlim]) - else: - medmag = 0.0 - sedmodel_abmag = np.zeros_like(templatecache['templatewave']) + medmag - else: - factor = 10**(0.4 * 48.6) * sedwave**2 / (C_LIGHT * 1e13) / FLUXNORM / CTools.massnorm # [erg/s/cm2/A --> maggies] - sedmodel_abmag = -2.5*np.log10(sedmodel * factor) - sedax.plot(sedwave / 1e4, sedmodel_abmag, color='slategrey', alpha=0.9, zorder=1) - - sedax.scatter(sedphot['lambda_eff']/1e4, sedphot['abmag'], marker='s', - s=400, color='k', facecolor='none', linewidth=2, alpha=1.0, zorder=2) - - if not fastphot: - #factor = 10**(0.4 * 48.6) * fullwave**2 / (C_LIGHT * 1e13) / FLUXNORM # [erg/s/cm2/A --> maggies] - #good = fullmodelspec > 0 - #sedax.plot(fullwave[good]/1e4, -2.5*np.log10(fullmodelspec[good]*factor[good]), color='purple', alpha=0.8) - for icam in np.arange(len(data['cameras'])): - factor = 10**(0.4 * 48.6) * data['wave'][icam]**2 / (C_LIGHT * 1e13) / FLUXNORM # [erg/s/cm2/A --> maggies] - good = desimodelspec[icam] > 0 - _wave = data['wave'][icam][good]/1e4 - _flux = -2.5*np.log10(desimodelspec[icam][good]*factor[good]) - sedax.plot(_wave, _flux, color=col2[icam], alpha=0.8) - #pdb.set_trace() - - # we have to set the limits *before* we call errorbar, below! - dm = 1.5 - sed_ymin = np.nanmax(sedmodel_abmag) + dm - sed_ymax = np.nanmin(sedmodel_abmag) - dm - if np.sum(abmag_good) > 0 and np.sum(abmag_goodlim) > 0: - sed_ymin = np.max((np.nanmax(phot['abmag'][abmag_good]), np.nanmax(phot['abmag_limit'][abmag_goodlim]), np.nanmax(sedmodel_abmag))) + dm - sed_ymax = np.min((np.nanmin(phot['abmag'][abmag_good]), np.nanmin(phot['abmag_limit'][abmag_goodlim]), np.nanmin(sedmodel_abmag))) - dm - elif np.sum(abmag_good) > 0 and np.sum(abmag_goodlim) == 0: - sed_ymin = np.max((np.nanmax(phot['abmag'][abmag_good]), np.nanmax(sedmodel_abmag))) + dm - sed_ymax = np.min((np.nanmin(phot['abmag'][abmag_good]), np.nanmin(sedmodel_abmag))) - dm - elif np.sum(abmag_good) == 0 and np.sum(abmag_goodlim) > 0: - sed_ymin = np.max((np.nanmax(phot['abmag_limit'][abmag_goodlim]), np.nanmax(sedmodel_abmag))) + dm - sed_ymax = np.min((np.nanmin(phot['abmag_limit'][abmag_goodlim]), np.nanmin(sedmodel_abmag))) - dm - else: - abmag_good = phot['abmag'] > 0 - abmag_goodlim = phot['abmag_limit'] > 0 - if np.sum(abmag_good) > 0 and np.sum(abmag_goodlim) > 0: - sed_ymin = np.max((np.nanmax(phot['abmag'][abmag_good]), np.nanmax(phot['abmag_limit'][abmag_goodlim]))) + dm - sed_ymax = np.min((np.nanmin(phot['abmag'][abmag_good]), np.nanmin(phot['abmag_limit'][abmag_goodlim]))) - dm - elif np.sum(abmag_good) > 0 and np.sum(abmag_goodlim) == 0: - sed_ymin = np.nanmax(phot['abmag'][abmag_good]) + dm - sed_ymax = np.nanmin(phot['abmag'][abmag_good]) - dm - elif np.sum(abmag_good) == 0 and np.sum(abmag_goodlim) > 0: - sed_ymin = np.nanmax(phot['abmag_limit'][abmag_goodlim]) + dm - sed_ymax = np.nanmin(phot['abmag_limit'][abmag_goodlim]) - dm - #else: - # sed_ymin, sed_ymax = [30, 20] - - if sed_ymin > 30: - sed_ymin = 30 - if np.isnan(sed_ymin) or np.isnan(sed_ymax): - raise('Problem here!') - #print(sed_ymin, sed_ymax) - - #sedax.set_xlabel(r'Observed-frame Wavelength ($\mu$m)') - sedax.set_xlim(phot_wavelims[0], phot_wavelims[1]) - sedax.set_xscale('log') - sedax.set_ylabel('AB mag') - #sedax.set_ylabel(r'Apparent Brightness (AB mag)') - sedax.set_ylim(sed_ymin, sed_ymax) - - sedax.xaxis.set_major_formatter(major_formatter) - obsticks = np.array([0.1, 0.2, 0.5, 1.0, 1.5, 3.0, 5.0, 10.0, 20.0]) - obsticks = obsticks[(obsticks >= phot_wavelims[0]) * (obsticks <= phot_wavelims[1])] - sedax.set_xticks(obsticks) - - if fastphot: - sedax.set_xlabel(r'Observed-frame Wavelength ($\mu$m)') - - sedax_twin = sedax.twiny() - sedax_twin.set_xlim(phot_wavelims[0]/(1+redshift), phot_wavelims[1]/(1+redshift)) - sedax_twin.set_xscale('log') - #sedax_twin.set_xlabel(r'Rest-frame Wavelength ($\mu$m)') - - sedax_twin.xaxis.set_major_formatter(major_formatter) - restticks = np.array([0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 1.5, 3.0, 5.0, 10.0, 15.0, 20.0]) - restticks = restticks[(restticks >= phot_wavelims[0]/(1+redshift)) * (restticks <= phot_wavelims[1]/(1+redshift))] - sedax_twin.set_xticks(restticks) - - # integrated flux / photometry - abmag = np.squeeze(phot['abmag']) - abmag_limit = np.squeeze(phot['abmag_limit']) - abmag_fainterr = np.squeeze(phot['abmag_fainterr']) - abmag_brighterr = np.squeeze(phot['abmag_brighterr']) - yerr = np.squeeze([abmag_brighterr, abmag_fainterr]) - - markersize = 14 - - dofit = np.where(CTools.bands_to_fit)[0] - if len(dofit) > 0: - good = np.where((abmag[dofit] > 0) * (abmag_limit[dofit] == 0))[0] - upper = np.where(abmag_limit[dofit] > 0)[0] - if len(good) > 0: - sedax.errorbar(phot['lambda_eff'][dofit][good]/1e4, abmag[dofit][good], - yerr=yerr[:, dofit[good]], - fmt='o', markersize=markersize, markeredgewidth=1, markeredgecolor='k', - markerfacecolor=photcol1, elinewidth=3, ecolor=photcol1, capsize=4, - label=r'$grz\,W_{1}W_{2}W_{3}W_{4}$', zorder=2, alpha=1.0) - if len(upper) > 0: - sedax.errorbar(phot['lambda_eff'][dofit][upper]/1e4, abmag_limit[dofit][upper], - lolims=True, yerr=0.75, - fmt='o', markersize=markersize, markeredgewidth=3, markeredgecolor='k', - markerfacecolor=photcol1, elinewidth=3, ecolor=photcol1, capsize=4, alpha=0.7) - - ignorefit = np.where(CTools.bands_to_fit == False)[0] - if len(ignorefit) > 0: - good = np.where((abmag[ignorefit] > 0) * (abmag_limit[ignorefit] == 0))[0] - upper = np.where(abmag_limit[ignorefit] > 0)[0] - if len(good) > 0: - sedax.errorbar(phot['lambda_eff'][ignorefit][good]/1e4, abmag[ignorefit][good], - yerr=yerr[:, ignorefit[good]], - fmt='o', markersize=markersize, markeredgewidth=3, markeredgecolor='k', - markerfacecolor='none', elinewidth=3, ecolor=photcol1, capsize=4, alpha=0.7) - if len(upper) > 0: - sedax.errorbar(phot['lambda_eff'][ignorefit][upper]/1e4, abmag_limit[ignorefit][upper], - lolims=True, yerr=0.75, fmt='o', markersize=markersize, markeredgewidth=3, - markeredgecolor='k', markerfacecolor='none', elinewidth=3, - ecolor=photcol1, capsize=5, alpha=0.7) - - if fastphot: - txt = leg['rchi2_phot'] - else: - # Label the DESI wavelength range and the aperture correction. - sedax.plot([np.min(fullwave)/1e4, np.max(fullwave)/1e4], [sed_ymin-1, sed_ymin-1], - lw=2, ls='-', color='gray', marker='s')#, alpha=0.5) - sedax.text(((np.max(fullwave)-np.min(fullwave))/2+np.min(fullwave)*0.8)/1e4, sed_ymin-1.7, - 'DESI x {:.2f}'.format(apercorr), ha='center', va='center', fontsize=16, - color='k') - - txt = '\n'.join((leg['rchi2_cont'], leg['rchi2_phot'], leg['rchi2'])) - - sedax.text(0.02, 0.94, txt, ha='left', va='top', - transform=sedax.transAxes, fontsize=legfntsz)#, bbox=bbox) - - txt = '\n'.join(( - #r'{}'.format(leg['fagn']), - r'{}'.format(leg['zzsun']), - r'{}'.format(leg['AV']), - r'{}'.format(leg['sfr']), - r'{}'.format(leg['age']), - r'{}'.format(leg['mstar']), - )) - sedax.text(legxpos, legypos2, txt, ha='right', va='bottom', - transform=sedax.transAxes, fontsize=legfntsz1, bbox=bbox) - - if not fastphot: - # draw lines connecting the SED and spectral plots - sedax.add_artist(ConnectionPatch(xyA=(spec_wavelims[0]/1e4, sed_ymin), - xyB=(spec_wavelims[0]/1e4, spec_ymax), - coordsA='data', coordsB='data', - axesA=sedax, axesB=specax, color='k')) - sedax.add_artist(ConnectionPatch(xyA=(spec_wavelims[1]/1e4, sed_ymin), - xyB=(spec_wavelims[1]/1e4, spec_ymax), - coordsA='data', coordsB='data', - axesA=sedax, axesB=specax, color='k')) - - # zoom in on individual emission lines - use linetable! - if not fastphot: - linetable = EMFit.linetable - inrange = (linetable['restwave'] * (1+redshift) > np.min(fullwave)) * (linetable['restwave'] * (1+redshift) < np.max(fullwave)) - linetable = linetable[inrange] - - nline = len(set(linetable['plotgroup'])) - - plotsig_default = 200.0 # [km/s] - plotsig_default_balmer = 500.0 # [km/s] - plotsig_default_broad = 2000.0 # [km/s] - - minwaves, maxwaves, meanwaves, deltawaves, sigmas, linenames = [], [], [], [], [], [] - for plotgroup in set(linetable['plotgroup']): - I = np.where(plotgroup == linetable['plotgroup'])[0] - linename = linetable['nicename'][I[0]].replace('-', ' ') - linenames.append(linename) - minwaves.append(np.min(linetable['restwave'][I])) - maxwaves.append(np.max(linetable['restwave'][I])) - meanwaves.append(np.mean(linetable['restwave'][I])) - deltawaves.append((np.max(linetable['restwave'][I]) - np.min(linetable['restwave'][I])) / 2) - - sigmas1 = np.array([fastspec['{}_SIGMA'.format(line.upper())] for line in linetable[I]['name']]) - sigmas1 = sigmas1[sigmas1 > 0] - if len(sigmas1) > 0: - plotsig = 1.5*np.mean(sigmas1) - if plotsig < 50: - plotsig = 50.0 - else: - if np.any(linetable['isbroad'][I]): - if np.any(linetable['isbalmer'][I]): - plotsig = fastspec['BROAD_SIGMA'] - if plotsig < 50: - plotsig = fastspec['NARROW_SIGMA'] - if plotsig < 50: - plotsig = plotsig_default - #plotsig = plotsig_default_broad - else: - plotsig = fastspec['UV_SIGMA'] - if plotsig < 50: - plotsig = plotsig_default_broad - else: - plotsig = fastspec['NARROW_SIGMA'] - if plotsig < 50: - plotsig = plotsig_default - sigmas.append(plotsig) - - if len(linetable) == 0: - ax = [] - else: - srt = np.argsort(meanwaves) - minwaves = np.hstack(minwaves)[srt] - maxwaves = np.hstack(maxwaves)[srt] - meanwaves = np.hstack(meanwaves)[srt] - deltawaves = np.hstack(deltawaves)[srt] - sigmas = np.hstack(sigmas)[srt] - linenames = np.hstack(linenames)[srt] - - # Add the linenames to the spectrum plot. - for meanwave, linename in zip(meanwaves*(1+redshift), linenames): - #print(meanwave, ymax_spec) - if meanwave > spec_wavelims[0] and meanwave < spec_wavelims[1]: - if 'SiIII' in linename: - thislinename = '\n'+linename.replace('+', '+\n ') - elif '4363' in linename: - thislinename = linename+'\n' - else: - thislinename = linename - if stackfit: - specax.text(meanwave/1e4, spec_ymax*0.97, thislinename, ha='center', va='top', - rotation=270, fontsize=12, alpha=0.5) - else: - specax.text(meanwave/1e4, spec_ymax, thislinename, ha='center', va='top', - rotation=270, fontsize=12, alpha=0.5) - - removelabels = np.ones(nline, bool) - line_ymin, line_ymax = np.zeros(nline)+1e6, np.zeros(nline)-1e6 - - if stackfit: - ax, irow, colshift = [], 0, 5 - else: - ax, irow, colshift = [], 4, 5 # skip the gap row - - for iax, (minwave, maxwave, meanwave, deltawave, sig, linename) in enumerate( - zip(minwaves, maxwaves, meanwaves, deltawaves, sigmas, linenames)): - icol = iax % nlinecols - icol += colshift - if iax > 0 and iax % nlinecols == 0: - irow += 1 - #print(iax, irow, icol) - - xx = fig.add_subplot(gs[irow, icol]) - ax.append(xx) - - wmin = (minwave - deltawave) * (1+redshift) - 6 * sig * minwave * (1+redshift) / C_LIGHT - wmax = (maxwave + deltawave) * (1+redshift) + 6 * sig * maxwave * (1+redshift) / C_LIGHT - #wmin = (meanwave - deltawave) * (1+redshift) - 6 * sig * meanwave * (1+redshift) / C_LIGHT - #wmax = (meanwave + deltawave) * (1+redshift) + 6 * sig * meanwave * (1+redshift) / C_LIGHT - - # iterate over cameras - for icam in np.arange(len(data['cameras'])): # iterate over cameras - emlinewave = data['wave'][icam] - emlineflux = data['flux'][icam] - desicontinuum[icam] - desismoothcontinuum[icam] - emlinemodel = desiemlines[icam] - - emlinesigma, good = ivar2var(data['ivar'][icam], sigma=True, allmasked_ok=True, clip=0) - emlinewave = emlinewave[good] - emlineflux = emlineflux[good] - emlinesigma = emlinesigma[good] - emlinemodel = emlinemodel[good] - - emlinemodel_oneline = [] - for desiemlines_oneline1 in desiemlines_oneline: - emlinemodel_oneline.append(desiemlines_oneline1[icam][good]) - - indx = np.where((emlinewave > wmin) * (emlinewave < wmax))[0] - if len(indx) > 1: - removelabels[iax] = False - xx.plot(emlinewave[indx]/1e4, emlineflux[indx], color=col1[icam], alpha=0.5) - #if nsmoothspec > 1: - # xx.plot(emlinewave[indx]/1e4, gaussian_filter(emlineflux[indx], nsmoothspec), color=col1[icam], alpha=0.5) - #else: - # xx.plot(emlinewave[indx]/1e4, emlineflux[indx], color=col1[icam], alpha=0.5) - for emlinemodel_oneline1 in emlinemodel_oneline: - if np.sum(emlinemodel_oneline1[indx]) > 0: - xx.plot(emlinewave[indx]/1e4, emlinemodel_oneline1[indx], lw=1, alpha=0.8, color=col2[icam]) - #if nsmoothspec > 1: - # xx.plot(emlinewave[indx]/1e4, gaussian_filter(emlinemodel_oneline1[indx], nsmoothspec), lw=1, alpha=0.8, color=col2[icam]) - #else: - # xx.plot(emlinewave[indx]/1e4, emlinemodel_oneline1[indx], lw=1, alpha=0.8, color=col2[icam]) - xx.plot(emlinewave[indx]/1e4, emlinemodel[indx], color=col2[icam], lw=2, alpha=0.8) - #if nsmoothspec > 1: - # xx.plot(emlinewave[indx]/1e4, emlinemodel[indx], color=col2[icam], lw=2, alpha=0.8) - #else: - # xx.plot(emlinewave[indx]/1e4, gaussian_filter(emlinemodel[indx], nsmoothspec), color=col2[icam], lw=2, alpha=0.8) - - #xx.plot(emlinewave[indx], emlineflux[indx]-emlinemodel[indx], color='gray', alpha=0.3) - #xx.axhline(y=0, color='gray', ls='--') - - # get the robust range - sigflux = np.std(emlineflux[indx]) - filtflux = median_filter(emlineflux[indx], 3, mode='nearest') - - #_line_ymin, _line_ymax = -1.5 * sigflux, 4 * sigflux - #if np.max(emlinemodel[indx]) > _line_ymax: - # _line_ymax = np.max(emlinemodel[indx]) * 1.3 - _line_ymin, _line_ymax = -1.5 * sigflux, np.max(emlinemodel[indx]) * 1.4 - if 4 * sigflux > _line_ymax: - _line_ymax = 4 * sigflux - if np.max(filtflux) > _line_ymax: - _line_ymax = np.max(filtflux) - if np.min(emlinemodel[indx]) < _line_ymin: - _line_ymin = 0.8 * np.min(emlinemodel[indx]) - if _line_ymax > line_ymax[iax]: - line_ymax[iax] = _line_ymax - if _line_ymin < line_ymin[iax]: - line_ymin[iax] = _line_ymin - #print(linename, line_ymin[iax], line_ymax[iax]) - #if linename == '[OII] $\lambda\lambda$3726,29': - # pdb.set_trace() - - xx.set_xlim(wmin/1e4, wmax/1e4) - - xx.text(0.03, 0.89, linename, ha='left', va='center', - transform=xx.transAxes, fontsize=12) - xx.tick_params(axis='x', labelsize=16) - xx.tick_params(axis='y', labelsize=16) - - for iax, xx in enumerate(ax): - if removelabels[iax]: - xx.set_ylim(0, 1) - xx.set_xticklabels([]) - xx.set_yticklabels([]) - else: - xx.set_yticklabels([]) - xx.set_ylim(line_ymin[iax], line_ymax[iax]) - xx_twin = xx.twinx() - xx_twin.set_ylim(line_ymin[iax], line_ymax[iax]) - xlim = xx.get_xlim() - xx.xaxis.set_major_locator(ticker.MaxNLocator(2)) - - if fastphot: - plt.subplots_adjust(wspace=0.6, top=0.85, bottom=0.13, left=0.07, right=0.86) - - else: - plt.subplots_adjust(wspace=0.4, top=0.9, bottom=0.1, left=0.07, right=0.92, hspace=0.33) - - # common axis labels - if len(ax) > 0: - if len(ax) == 2: - xx = fig.add_subplot(gs[irow, icol+1]) - xx.axis('off') - ax.append(xx) - elif len(ax) == 1: - xx = fig.add_subplot(gs[irow, icol+1]) - xx.axis('off') - ax.append(xx) - xx = fig.add_subplot(gs[irow, icol+2]) - xx.axis('off') - ax.append(xx) - - ulpos = ax[0].get_position() - lpos = ax[nline-1].get_position() - urpos = ax[2].get_position() - xpos = (urpos.x1 - ulpos.x0) / 2 + ulpos.x0# + 0.03 - - ypos = lpos.y0 - 0.04 - fig.text(xpos, ypos, r'Observed-frame Wavelength ($\mu$m)', - ha='center', va='center', fontsize=fontsize2) - - xpos = urpos.x1 + 0.05 - ypos = (urpos.y1 - lpos.y0) / 2 + lpos.y0# + 0.03 - fig.text(xpos, ypos, - r'$F_{\lambda}\ (10^{-17}~{\rm erg}~{\rm s}^{-1}~{\rm cm}^{-2}~\AA^{-1})$', - ha='center', va='center', rotation=270, fontsize=fontsize2) - - legkeys = leg.keys() - legfntsz = 17 - - # rest wavelength plus targeting information - if fastphot: - ppos = sedax.get_position() - xpos = (ppos.x1 - ppos.x0) / 2 + ppos.x0 - ypos = ppos.y1 + 0.07 - - fig.text(xpos, ypos, r'Rest-frame Wavelength ($\mu$m)', - ha='center', va='bottom', fontsize=fontsize2) - fig.text(0.58, 0.87, '\n'.join(target), ha='left', va='bottom', - fontsize=fontsize2, linespacing=1.4) - - toppos, leftpos = 0.25, 0.59 - - txt = [ - r'{}'.format(leg['z']), - ] - if 'zredrock' in legkeys: - txt += [r'{}'.format(leg['zredrock'])] - txt += [ - #r'{}'.format(leg['zwarn']), - #'', - r'{}'.format(leg['vdisp']), - '', - r'{}'.format(leg['dn4000_model']), - ] - - fig.text(leftpos, toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, - bbox=bbox, linespacing=1.4) - - txt = [r'{}'.format(leg['absmag_r'])] - if 'absmag_gr' in legkeys: - txt += [r'{}'.format(leg['absmag_gr'])] - if 'absmag_rz' in legkeys: - txt += [r'{}'.format(leg['absmag_rz'])] - - fig.text(leftpos+0.18, toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, - bbox=bbox, linespacing=1.4) - - else: - if stackfit: - ppos = specax.get_position() - xpos = (ppos.x1 - ppos.x0) / 2 + ppos.x0 - ypos = ppos.y1 + 0.05 - - fig.text(xpos, ypos, r'Rest-frame Wavelength ($\mu$m)', - ha='center', va='bottom', fontsize=fontsize2) - fig.text(0.62, 0.91, '\n'.join(target), ha='left', va='bottom', - fontsize=fontsize2, linespacing=1.4) - else: - ppos = sedax.get_position() - xpos = (ppos.x1 - ppos.x0) / 2 + ppos.x0 - ypos = ppos.y1 + 0.03 - - fig.text(xpos, ypos, r'Rest-frame Wavelength ($\mu$m)', - ha='center', va='bottom', fontsize=fontsize2) - fig.text(0.647, 0.925, '\n'.join(target), ha='left', va='bottom', - fontsize=fontsize2, linespacing=1.4) - - # add some key results about the object at the bottom of the figure - - if stackfit: - #toppos, startpos, deltapos = 0.21, 0.04, 0.13 - toppos, leftpos, rightpos, adjust = 0.27, 0.05, 0.62, 0.01 - else: - toppos, leftpos, rightpos, adjust = 0.21, 0.03, 0.62, 0.01 - - nbox = 2 + 1*bool(leg_narrow) + bool(leg_broad) - boxpos = np.arange(nbox) * (rightpos - leftpos)/nbox + leftpos - - txt = [ - r'{}'.format(leg['z']), - ] - if 'zredrock' in legkeys: - txt += [r'{}'.format(leg['zredrock'])] - - txt += [ - #r'{}'.format(leg['zwarn']), - r'{}'.format(leg['vdisp']), - ] - if 'dv_narrow' in legkeys or 'dv_uv' in leg_uv.keys() or 'dv_broad' in leg_broad.keys(): - txt += [''] - - if 'dv_narrow' in legkeys: - txt += [ - r'{}'.format(leg['sigma_narrow']), - r'{}'.format(leg['dv_narrow']), - ] - if 'dv_uv' in leg_uv.keys(): - txt += [ - r'{}'.format(leg_uv['sigma_uv']), - r'{}'.format(leg_uv['dv_uv']), - ] - _ = leg_uv.pop('sigma_uv') - _ = leg_uv.pop('dv_uv') - if 'dv_broad' in leg_broad.keys(): - txt += [ - r'{}'.format(leg_broad['sigma_broad']), - r'{}'.format(leg_broad['dv_broad']), - ] - _ = leg_broad.pop('sigma_broad') - _ = leg_broad.pop('dv_broad') - - ibox = 0 - #fig.text(startpos, toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, - fig.text(boxpos[ibox], toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, - bbox=bbox, linespacing=1.4) - ibox += 1 - - txt = [ - r'{}'.format(leg['absmag_r']), - r'{}'.format(leg['absmag_gr']), - r'{}'.format(leg['absmag_rz']), - '', - r'{}'.format(leg['dn4000_model']) - ] - if 'dn4000_spec' in legkeys: - txt += [r'{}'.format(leg['dn4000_spec'])] - - #fig.text(startpos+deltapos, toppos, '\n'.join(txt), ha='left', va='top', - fig.text(boxpos[ibox], toppos, '\n'.join(txt), ha='left', va='top', - fontsize=legfntsz, bbox=bbox, linespacing=1.4) - ibox += 1 - - #factor = 2 - if bool(leg_narrow): - txt = [] - for key in leg_narrow.keys(): - txt += [r'{}'.format(leg_narrow[key])] - #fig.text(startpos+deltapos*factor, toppos, '\n'.join(txt), ha='left', va='top', - fig.text(boxpos[ibox]-adjust*2, toppos, '\n'.join(txt), ha='left', va='top', - fontsize=legfntsz, bbox=bbox, linespacing=1.4) - ibox += 1 - #factor += 1.25 - - if bool(leg_broad): - txt = [] - for key in leg_broad.keys(): - txt += [r'{}'.format(leg_broad[key])] - - if bool(leg_uv): - if bool(leg_broad): - txt += [''] - else: - txt = [] - for key in leg_uv.keys(): - txt += [r'{}'.format(leg_uv[key])] - - if bool(leg_uv) or bool(leg_broad): - #fig.text(startpos+deltapos*factor, toppos, '\n'.join(txt), ha='left', va='top', - fig.text(boxpos[ibox]-adjust*1, toppos, '\n'.join(txt), ha='left', va='top', - fontsize=legfntsz, bbox=bbox, linespacing=1.4) - - log.info('Writing {}'.format(pngfile)) - fig.savefig(pngfile)#, dpi=150) - plt.close() diff --git a/py/fastspecfit/io.py b/py/fastspecfit/io.py index b94102ac..e838e3f9 100644 --- a/py/fastspecfit/io.py +++ b/py/fastspecfit/io.py @@ -19,10 +19,10 @@ log = get_logger() # Default environment variables. -DESI_ROOT_NERSC = '/dvs_ro/cfs/cdirs/desi' -DUST_DIR_NERSC = '/dvs_ro/cfs/cdirs/cosmo/data/dust/v0_1' -FPHOTO_DIR_NERSC = '/dvs_ro/cfs/cdirs/desi/external/legacysurvey/dr9' -FTEMPLATES_DIR_NERSC = '/dvs_ro/cfs/cdirs/desi/science/gqp/templates/fastspecfit' +DESI_ROOT_NERSC = '/global/cfs/cdirs/desi' +DUST_DIR_NERSC = '/global/cfs/cdirs/cosmo/data/dust/v0_1' +FPHOTO_DIR_NERSC = '/global/cfs/cdirs/desi/external/legacysurvey/dr9' +FTEMPLATES_DIR_NERSC = '/global/cfs/cdirs/desi/science/gqp/templates/fastspecfit' # list of all possible targeting bit columns TARGETINGBITS = { @@ -1656,16 +1656,19 @@ def init_fastspec_output(input_meta, specprod, fphoto=None, templates=None, for band in fphoto['bands']: out.add_column(Column(name='FLUX_SYNTH_PHOTMODEL_{}'.format(band.upper()), length=nobj, dtype='f4', unit='nanomaggies')) - for band in fphoto['absmag_bands']: - out.add_column(Column(name='KCORR_{}'.format(band.upper()), length=nobj, dtype='f4', unit=u.mag)) - out.add_column(Column(name='ABSMAG_{}'.format(band.upper()), length=nobj, dtype='f4', unit=u.mag)) # absolute magnitudes - out.add_column(Column(name='ABSMAG_IVAR_{}'.format(band.upper()), length=nobj, dtype='f4', unit=1/u.mag**2)) + for band, shift in zip(fphoto['absmag_bands'], fphoto['band_shift']): + out.add_column(Column(name='ABSMAG{:02d}_{}'.format(int(10*shift), band.upper()), length=nobj, dtype='f4', unit=u.mag)) # absolute magnitudes + out.add_column(Column(name='ABSMAG{:02d}_IVAR_{}'.format(int(10*shift), band.upper()), length=nobj, dtype='f4', unit=1/u.mag**2)) + out.add_column(Column(name='KCORR{:02d}_{}'.format(int(10*shift), band.upper()), length=nobj, dtype='f4', unit=u.mag)) for cflux in ['LOGLNU_1500', 'LOGLNU_2800']: out.add_column(Column(name=cflux, length=nobj, dtype='f4', unit=10**(-28)*u.erg/u.second/u.Hz)) + out.add_column(Column(name='LOGL_1450', length=nobj, dtype='f4', unit=10**(10)*u.solLum)) + out.add_column(Column(name='LOGL_1700', length=nobj, dtype='f4', unit=10**(10)*u.solLum)) + out.add_column(Column(name='LOGL_3000', length=nobj, dtype='f4', unit=10**(10)*u.solLum)) out.add_column(Column(name='LOGL_5100', length=nobj, dtype='f4', unit=10**(10)*u.solLum)) - for cflux in ['FOII_3727_CONT', 'FHBETA_CONT', 'FOIII_5007_CONT', 'FHALPHA_CONT']: + for cflux in ['FLYA_1215_CONT', 'FOII_3727_CONT', 'FHBETA_CONT', 'FOIII_5007_CONT', 'FHALPHA_CONT']: out.add_column(Column(name=cflux, length=nobj, dtype='f4', unit=10**(-17)*u.erg/(u.second*u.cm**2*u.Angstrom))) if not fastphot: diff --git a/py/fastspecfit/mpi.py b/py/fastspecfit/mpi.py index f60cbaf8..e8b0581b 100644 --- a/py/fastspecfit/mpi.py +++ b/py/fastspecfit/mpi.py @@ -44,109 +44,6 @@ def get_ntargets_one(specfile, htmldir_root, outdir_root, coadd_type='healpix', ntargets = np.sum(J) return ntargets -def weighted_partition(weights, n): - """ - Partition ``weights`` into ``n`` groups with approximately same sum(weights). - - Args: - weights: array-like weights - n: number of groups - - Returns (groups, groupweights): - groups: list of lists of indices of weights for each group - groupweights: sum of weights assigned to each group - - - """ - sumweights = np.zeros(n, dtype=float) - groups = list() - for i in range(n): - groups.append(list()) - weights = np.asarray(weights) - for i in np.argsort(-weights): - j = np.argmin(sumweights) - groups[j].append(i) - sumweights[j] += weights[i] - - return groups, np.array([np.sum(x) for x in sumweights]) - -def group_redrockfiles(specfiles, maxnodes=256, comm=None, makeqa=False): - ''' - Group redrockfiles to balance runtimes - - Args: - specfiles: list of spectra filepaths - - Options: - maxnodes: split the spectra into this number of nodes - comm: MPI communicator - - Returns (groups, ntargets, grouptimes): - * groups: list of lists of indices to specfiles - * list of number of targets per group - * grouptimes: list of expected runtimes for each group - - ''' - import fitsio - - if comm is None: - rank, size = 0, 1 - else: - rank, size = comm.rank, comm.size - - npix = len(specfiles) - pixgroups = np.array_split(np.arange(npix), size) - ntargets = np.zeros(len(pixgroups[rank]), dtype=int) - for i, j in enumerate(pixgroups[rank]): - if makeqa: - ntargets[i] = fitsio.FITS(specfiles[j])[1].get_nrows() - else: - from fastspecfit.io import ZWarningMask - #ntargets[i] = fitsio.FITS(specfiles[j])[1].get_nrows() - zb = fitsio.read(specfile, 'REDSHIFTS', columns=['Z', 'ZWARN']) - fm = fitsio.read(specfile, 'FIBERMAP', columns=['TARGETID', 'OBJTYPE']) - J = ((zb['Z'] > 0.001) * (fm['OBJTYPE'] == 'TGT') * (zb['ZWARN'] & ZWarningMask.NODATA == 0)) - nt = np.sum(J) - ntargets[i] = nt - log.debug(i, j, specfiles[j], ntargets[i]) - - if False: - if comm is not None: - ntargets = comm.gather(ntargets) - if rank == 0: - ntargets = np.concatenate(ntargets) - ntargets = comm.bcast(ntargets, root=0) - - sumntargets = np.sum(sumntargets) - runtimes = 30 + 0.4*sumntargets - - # Aim for 25 minutes, but don't exceed maxnodes number of nodes. - ntime = 25 - if comm is not None: - numnodes = comm.size - else: - numnodes = min(maxnodes, int(np.ceil(np.sum(runtimes)/(ntime*60)))) - - groups, grouptimes = weighted_partition(runtimes, numnodes) - ntargets = np.array([np.sum(ntargets[ii]) for ii in groups]) - else: - groups = pixgroups - grouptimes = None - - return groups, ntargets, grouptimes - -def backup_logs(logfile): - ''' - Move logfile -> logfile.0 or logfile.1 or logfile.n as needed - - TODO: make robust against logfile.abc also existing - ''' - logfiles = glob(logfile+'.*') - newlog = logfile+'.'+str(len(logfiles)) - assert not os.path.exists(newlog) - os.rename(logfile, newlog) - return newlog - def plan(comm=None, specprod=None, specprod_dir=None, coadd_type='healpix', survey=None, program=None, healpix=None, tile=None, night=None, outdir_data='.', outdir_html='.', mp=1, merge=False, makeqa=False, @@ -154,6 +51,7 @@ def plan(comm=None, specprod=None, specprod_dir=None, coadd_type='healpix', import fitsio from astropy.table import Table, vstack + from desispec.parallel import weighted_partition from fastspecfit.io import DESI_ROOT_NERSC t0 = time.time() @@ -322,23 +220,25 @@ def _findfiles(filedir, prefix='redrock', survey=None, program=None, healpix=Non itodo = np.where(ntargets > 0)[0] if len(itodo) > 0: ntargets = ntargets[itodo] - indices = np.arange(len(ntargets)) if redrockfiles is not None: redrockfiles = redrockfiles[itodo] if outfiles is not None: outfiles = outfiles[itodo] - - # Assign the sample to ranks to make the ntargets distribution per rank ~flat. - # https://stackoverflow.com/questions/33555496/split-array-into-equally-weighted-chunks-based-on-order - cumuweight = ntargets.cumsum() / ntargets.sum() - idx = np.searchsorted(cumuweight, np.linspace(0, 1, size, endpoint=False)[1:]) - if len(idx) < size: # can happen in corner cases or with 1 rank - groups = np.array_split(indices, size) # unweighted - else: - groups = np.array_split(indices, idx) # weighted - for ii in range(size): # sort by weight - srt = np.argsort(ntargets[groups[ii]]) - groups[ii] = groups[ii][srt] + + groups = weighted_partition(ntargets, size) + + ## Assign the sample to ranks to make the ntargets distribution per rank ~flat. + ## https://stackoverflow.com/questions/33555496/split-array-into-equally-weighted-chunks-based-on-order + #indices = np.arange(len(ntargets)) + #cumuweight = ntargets.cumsum() / ntargets.sum() + #idx = np.searchsorted(cumuweight, np.linspace(0, 1, size, endpoint=False)[1:]) + #if len(idx) < size: # can happen in corner cases or with 1 rank + # groups = np.array_split(indices, size) # unweighted + #else: + # groups = np.array_split(indices, idx) # weighted + #for ii in range(size): # sort by weight + # srt = np.argsort(ntargets[groups[ii]]) + # groups[ii] = groups[ii][srt] else: if redrockfiles is not None: redrockfiles = [] @@ -360,11 +260,6 @@ def _findfiles(filedir, prefix='redrock', survey=None, program=None, healpix=Non # hack--build the output directories and pass them in the 'redrockfiles' # position! for coadd_type==cumulative, strip out the 'lastnight' argument if coadd_type == 'cumulative': - #redrockfiles = [] - #for outfile in outfiles: - # dd = os.path.split(outfile) - # redrockfiles.append(os.path.dirname(dd[0]).replace(outdir, htmldir)) - # os.path.dirname(dd[0]) redrockfiles = np.array([os.path.dirname(os.path.dirname(outfile)).replace(outdir, htmldir) for outfile in outfiles]) else: redrockfiles = np.array([os.path.dirname(outfile).replace(outdir, htmldir) for outfile in outfiles]) diff --git a/py/fastspecfit/qa.py b/py/fastspecfit/qa.py new file mode 100644 index 00000000..08f5b885 --- /dev/null +++ b/py/fastspecfit/qa.py @@ -0,0 +1,1431 @@ +#!/usr/bin/env python +""" +fastspecfit.qa +============== + +""" +import pdb # for debugging + +import os, time +import numpy as np + +from desiutil.log import get_logger +log = get_logger() + +def _desiqa_one(args): + """Multiprocessing wrapper.""" + return desiqa_one(*args) + +def desiqa_one(data, fastfit, metadata, templates, coadd_type, fphoto, + minspecwave=3500., maxspecwave=9900., minphotwave=0.1, + maxphotwave=35., emline_snrmin=0.0, nsmoothspec=1, + fastphot=False, nophoto=False, stackfit=False, inputz=False, + no_smooth_continuum=False, outdir=None, outprefix=None, log=None): + """Multiprocessing wrapper to generate QA for a single object. + + """ + from fastspecfit.io import cache_templates + templatecache = cache_templates(templates, log=log, mintemplatewave=450.0, + maxtemplatewave=40e4, fastphot=fastphot) + + if inputz: + from fastspecfit.util import TabulatedDESI + cosmo = TabulatedDESI() + else: + cosmo = None + + qa_fastspec(data, templatecache, fastfit, metadata, coadd_type=coadd_type, + spec_wavelims=(minspecwave, maxspecwave), + phot_wavelims=(minphotwave, maxphotwave), + no_smooth_continuum=no_smooth_continuum, + emline_snrmin=emline_snrmin, nsmoothspec=nsmoothspec, + fastphot=fastphot, fphoto=fphoto, stackfit=stackfit, + outprefix=outprefix, outdir=outdir, log=log, cosmo=cosmo) + +def qa_fastspec(data, templatecache, fastspec, metadata, coadd_type='healpix', + spec_wavelims=(3550, 9900), phot_wavelims=(0.1, 35), + fastphot=False, fphoto=None, stackfit=False, outprefix=None, + no_smooth_continuum=False, emline_snrmin=0.0, nsmoothspec=1, + outdir=None, log=None, cosmo=None): + """QA plot the emission-line spectrum and best-fitting model. + + """ + import subprocess + from scipy.ndimage import median_filter + + import matplotlib.pyplot as plt + import matplotlib.ticker as ticker + from matplotlib import colors + from matplotlib.patches import Circle, Rectangle, ConnectionPatch + from matplotlib.lines import Line2D + import matplotlib.gridspec as gridspec + import matplotlib.image as mpimg + + import astropy.units as u + from astropy.io import fits + from astropy.wcs import WCS + import seaborn as sns + from PIL import Image, ImageDraw + + from fastspecfit.util import ivar2var, C_LIGHT + from fastspecfit.io import FLUXNORM, get_qa_filename + from fastspecfit.continuum import ContinuumTools + from fastspecfit.emlines import build_emline_model, EMFitTools + + if log is None: + from desiutil.log import get_logger + log = get_logger() + + Image.MAX_IMAGE_PIXELS = None + + sns.set(context='talk', style='ticks', font_scale=1.3)#, rc=rc) + + if stackfit: + col1 = [colors.to_hex(col) for col in ['violet']] + col2 = [colors.to_hex(col) for col in ['purple']] + else: + col1 = [colors.to_hex(col) for col in ['dodgerblue', 'darkseagreen', 'orangered']] + col2 = [colors.to_hex(col) for col in ['darkblue', 'darkgreen', 'darkred']] + + photcol1 = colors.to_hex('darkorange') + + @ticker.FuncFormatter + def major_formatter(x, pos): + if (x >= 0.01) and (x < 0.1): + return f'{x:.2f}' + elif (x >= 0.1) and (x < 1): + return f'{x:.1f}' + else: + return f'{x:.0f}' + + CTools = ContinuumTools(fphoto=fphoto, + continuum_pixkms=templatecache['continuum_pixkms'], + pixkms_wavesplit=templatecache['pixkms_wavesplit']) + if 'legacysurveydr' in fphoto.keys(): + layer = 'ls-{}'.format(fphoto['legacysurveydr']) + else: + layer = 'ls-dr9' + + if not fastphot: + EMFit = EMFitTools(minspecwave=spec_wavelims[0], maxspecwave=spec_wavelims[1]) + + filters = CTools.synth_filters[metadata['PHOTSYS']] + allfilters = CTools.filters[metadata['PHOTSYS']] + redshift = fastspec['Z'] + + # cosmo will be provided if input_redshifts are used + if cosmo is not None: + dlum = cosmo.luminosity_distance(redshift) + else: + dlum = data['dluminosity'] + + pngfile = get_qa_filename(metadata, coadd_type, outprefix=outprefix, + outdir=outdir, fastphot=fastphot, log=log) + + # some arrays to use for the legend + if coadd_type == 'healpix': + target = [ + 'Survey/Program/Healpix: {}/{}/{}'.format(metadata['SURVEY'], metadata['PROGRAM'], metadata['HEALPIX']), + 'TargetID: {}'.format(metadata['TARGETID']), + ] + elif coadd_type == 'cumulative': + target = [ + 'Tile/Night/Fiber: {}/{}/{}'.format(metadata['TILEID'], metadata['NIGHT'], metadata['FIBER']), + 'TargetID: {}'.format(metadata['TARGETID']), + ] + elif coadd_type == 'pernight': + target = [ + 'Tile/Night/Fiber: {}/{}/{}'.format(metadata['TILEID'], metadata['NIGHT'], metadata['FIBER']), + 'TargetID: {}'.format(metadata['TARGETID']), + ] + elif coadd_type == 'perexp': + target = [ + 'Tile/Night/Fiber: {}/{}/{}'.format(metadata['TILEID'], metadata['NIGHT'], metadata['FIBER']), + 'Night/Fiber: {}/{}'.format(metadata['NIGHT'], metadata['FIBER']), + 'TargetID: {}'.format(metadata['TARGETID']), + ] + elif coadd_type == 'custom': + target = [ + 'Survey/Program/Healpix: {}/{}/{}'.format(metadata['SURVEY'], metadata['PROGRAM'], metadata['HEALPIX']), + 'TargetID: {}'.format(metadata['TARGETID']), + ] + elif coadd_type == 'stacked': + target = [ + 'StackID: {}'.format(metadata['STACKID']), + '', + ] + else: + errmsg = 'Unrecognized coadd_type {}!'.format(coadd_type) + log.critical(errmsg) + raise ValueError(errmsg) + + leg = { + 'z': '$z={:.7f}$'.format(redshift), + 'rchi2_phot': r'$\chi^{2}_{\nu,\mathrm{phot}}=$'+r'${:.2f}$'.format(fastspec['RCHI2_PHOT']), + 'dn4000_model': r'$D_{n}(4000)_{\mathrm{model}}=$'+r'${:.3f}$'.format(fastspec['DN4000_MODEL']), + 'age': r'Age$={:.3f}$ Gyr'.format(fastspec['AGE']), + 'AV': r'$A_{V}=$'+r'${:.3f}$ mag'.format(fastspec['AV']), + 'mstar': r'$\log_{10}(M/M_{\odot})=$'+r'${:.3f}$'.format(fastspec['LOGMSTAR']), + 'sfr': r'$\mathrm{SFR}=$'+'${:.1f}$'.format(fastspec['SFR'])+' $M_{\odot}/\mathrm{yr}$', + 'zzsun': r'$Z/Z_{\odot}=$'+r'${:.3f}$'.format(fastspec['ZZSUN']), + } + + # try to figure out which absmags to display + gindx = np.argmin(np.abs(CTools.absmag_filters.effective_wavelengths.value - 4900)) + rindx = np.argmin(np.abs(CTools.absmag_filters.effective_wavelengths.value - 6500)) + zindx = np.argmin(np.abs(CTools.absmag_filters.effective_wavelengths.value - 9200)) + absmag_gband = CTools.absmag_bands[gindx] + absmag_rband = CTools.absmag_bands[rindx] + absmag_zband = CTools.absmag_bands[zindx] + shift_gband = CTools.band_shift[gindx] + shift_rband = CTools.band_shift[rindx] + shift_zband = CTools.band_shift[zindx] + + leg.update({'absmag_r': '$M_{{{}}}={:.2f}$'.format( + absmag_rband.lower().replace('decam_', '').replace('sdss_', ''), + fastspec['ABSMAG{:02d}_{}'.format(int(10*shift_gband), absmag_rband.upper())])}) + if gindx != rindx: + gr = (fastspec['ABSMAG{:02d}_{}'.format(int(10*shift_gband), absmag_gband.upper())] - + fastspec['ABSMAG{:02d}_{}'.format(int(10*shift_rband), absmag_rband.upper())]) + leg.update({'absmag_gr': '$M_{{{}}}-M_{{{}}}={:.3f}$'.format( + absmag_gband.lower(), absmag_rband.lower(), gr).replace('decam_', '').replace('sdss_', '')}) + if zindx != rindx: + rz = (fastspec['ABSMAG{:02d}_{}'.format(int(10*shift_rband), absmag_rband.upper())] - + fastspec['ABSMAG{:02d}_{}'.format(int(10*shift_zband), absmag_zband.upper())]) + leg.update({'absmag_rz': '$M_{{{}}}-M_{{{}}}={:.3f}$'.format( + absmag_rband.lower(), absmag_zband.lower(), rz).replace('decam_', '').replace('sdss_', '')}) + + #leg['radec'] = '$(\\alpha,\\delta)=({:.7f},{:.6f})$'.format(metadata['RA'], metadata['DEC']) + #leg['zwarn'] = '$z_{{\\rm warn}}={}$'.format(metadata['ZWARN']) + + if fastphot: + leg['vdisp'] = r'$\sigma_{star}=$'+'{:g}'.format(fastspec['VDISP'])+' km/s' + else: + if fastspec['VDISP_IVAR'] > 0: + leg['vdisp'] = r'$\sigma_{{star}}={:.0f}\pm{:.0f}$ km/s'.format(fastspec['VDISP'], 1/np.sqrt(fastspec['VDISP_IVAR'])) + else: + leg['vdisp'] = r'$\sigma_{{star}}={:g}$ km/s'.format(fastspec['VDISP']) + + leg['rchi2'] = r'$\chi^{2}_{\nu,\mathrm{specphot}}$='+'{:.2f}'.format(fastspec['RCHI2']) + leg['rchi2_cont'] = r'$\chi^{2}_{\nu,\mathrm{cont}}$='+'{:.2f}'.format(fastspec['RCHI2_CONT']) + + if not stackfit: + if redshift != metadata['Z_RR']: + leg['zredrock'] = '$z_{\mathrm{Redrock}}=$'+r'${:.7f}$'.format(metadata['Z_RR']) + + if fastphot: + fontsize1 = 16 + fontsize2 = 22 + else: + if stackfit: + fontsize1 = 16 + fontsize2 = 22 + else: + fontsize1 = 18 # 24 + fontsize2 = 24 + + apercorr = fastspec['APERCORR'] + + if fastspec['DN4000_IVAR'] > 0: + leg['dn4000_spec'] = r'$D_{n}(4000)_{\mathrm{data}}=$'+r'${:.3f}$'.format(fastspec['DN4000']) + + # kinematics + if fastspec['NARROW_Z'] != redshift: + if fastspec['NARROW_ZRMS'] > 0: + leg['dv_narrow'] = r'$\Delta v_{\mathrm{narrow}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( + C_LIGHT*(fastspec['NARROW_Z']-redshift), C_LIGHT*fastspec['NARROW_ZRMS']) + else: + leg['dv_narrow'] = r'$\Delta v_{\mathrm{narrow}}=$'+r'${:.0f}$ km/s'.format( + C_LIGHT*(fastspec['NARROW_Z']-redshift)) + if fastspec['NARROW_SIGMA'] != 0.0: + if fastspec['NARROW_SIGMARMS'] > 0: + leg['sigma_narrow'] = r'$\sigma_{\mathrm{narrow}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( + fastspec['NARROW_SIGMA'], fastspec['NARROW_SIGMARMS']) + else: + leg['sigma_narrow'] = r'$\sigma_{\mathrm{narrow}}=$'+r'${:.0f}$ km/s'.format(fastspec['NARROW_SIGMA']) + + snrcut = 1.5 + leg_broad, leg_narrow, leg_uv = {}, {}, {} + + if fastspec['UV_Z'] != redshift: + if fastspec['UV_ZRMS'] > 0: + leg_uv['dv_uv'] = r'$\Delta v_{\mathrm{UV}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( + C_LIGHT*(fastspec['UV_Z']-redshift), C_LIGHT*fastspec['UV_ZRMS']) + else: + leg_uv['dv_uv'] = r'$\Delta v_{\mathrm{UV}}=$'+r'${:.0f}$ km/s'.format( + C_LIGHT*(fastspec['UV_Z']-redshift)) + if fastspec['UV_SIGMA'] != 0.0: + if fastspec['UV_SIGMARMS'] > 0: + leg_uv['sigma_uv'] = r'$\sigma_{\mathrm{UV}}$'+r'$={:.0f}\pm{:.0f}$ km/s'.format( + fastspec['UV_SIGMA'], fastspec['UV_SIGMARMS']) + else: + leg_uv['sigma_uv'] = r'$\sigma_{\mathrm{UV}}=$'+r'${:.0f}$ km/s'.format(fastspec['UV_SIGMA']) + if fastspec['BROAD_Z'] != redshift: + if fastspec['BROAD_ZRMS'] > 0: + leg_broad['dv_broad'] = r'$\Delta v_{\mathrm{broad}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( + C_LIGHT*(fastspec['BROAD_Z']-redshift), C_LIGHT*fastspec['BROAD_ZRMS']) + else: + leg_broad['dv_broad'] = r'$\Delta v_{\mathrm{broad}}=$'+r'${:.0f}$ km/s'.format( + C_LIGHT*(fastspec['BROAD_Z']-redshift)) + if fastspec['BROAD_SIGMA'] != 0.0: + if fastspec['BROAD_SIGMARMS'] > 0: + leg_broad['sigma_broad'] = r'$\sigma_{\mathrm{broad}}=$'+r'${:.0f}\pm{:.0f}$ km/s'.format( + fastspec['BROAD_SIGMA'], fastspec['BROAD_SIGMARMS']) + else: + leg_broad['sigma_broad'] = r'$\sigma_{\mathrm{broad}}=$'+r'${:.0f}$ km/s'.format(fastspec['BROAD_SIGMA']) + + # emission lines + + # UV + if fastspec['LYALPHA_AMP']*np.sqrt(fastspec['LYALPHA_AMP_IVAR']) > snrcut: + leg_uv['ewlya'] = r'EW(Ly$\alpha)=$'+r'${:.1f}$'.format(fastspec['LYALPHA_EW'])+r' $\AA$' + if fastspec['CIV_1549_AMP']*np.sqrt(fastspec['CIV_1549_AMP_IVAR']) > snrcut: + leg_uv['ewciv'] = r'EW(CIV)=$'+r'${:.1f}$'.format(fastspec['CIV_1549_EW'])+r' $\AA$' + if fastspec['CIII_1908_AMP']*np.sqrt(fastspec['CIII_1908_AMP_IVAR']) > snrcut: + leg_uv['ewciii'] = r'EW(CIII])=$'+r'${:.1f}$'.format(fastspec['CIII_1908_EW'])+r' $\AA$' + if (fastspec['MGII_2796_AMP']*np.sqrt(fastspec['MGII_2796_AMP_IVAR']) > snrcut or + fastspec['MGII_2803_AMP']*np.sqrt(fastspec['MGII_2803_AMP_IVAR']) > snrcut): + leg_uv['ewmgii'] = r'EW(MgII)=$'+r'${:.1f}$'.format(fastspec['MGII_2796_EW']+fastspec['MGII_2803_EW'])+r' $\AA$' + leg_uv['mgii_doublet'] = r'MgII $\lambda2796/\lambda2803={:.3f}$'.format(fastspec['MGII_DOUBLET_RATIO']) + + leg_broad['linerchi2'] = r'$\chi^{2}_{\nu,\mathrm{line}}=$'+r'${:.2f}$'.format(fastspec['RCHI2_LINE']) + leg_broad['deltachi2'] = r'$\Delta\chi^{2}_{\mathrm{nobroad}}=$'+r'${:.0f}$'.format(fastspec['DELTA_LINECHI2']) + leg_broad['deltandof'] = r'$\Delta\nu_{\mathrm{nobroad}}=$'+r'${:.0f}$'.format(fastspec['DELTA_LINENDOF']) + + # choose one broad Balmer line + if fastspec['HALPHA_BROAD_AMP']*np.sqrt(fastspec['HALPHA_BROAD_AMP_IVAR']) > snrcut: + leg_broad['ewbalmer_broad'] = r'EW(H$\alpha)_{\mathrm{broad}}=$'+r'${:.1f}$'.format(fastspec['HALPHA_BROAD_EW'])+r' $\AA$' + elif fastspec['HBETA_BROAD_AMP']*np.sqrt(fastspec['HBETA_BROAD_AMP_IVAR']) > snrcut: + leg_broad['ewbalmer_broad'] = r'EW(H$\beta)_{\mathrm{broad}}=$'+r'${:.1f}$'.format(fastspec['HBETA_BROAD_EW'])+r' $\AA$' + elif fastspec['HGAMMA_BROAD_AMP']*np.sqrt(fastspec['HGAMMA_BROAD_AMP_IVAR']) > snrcut: + leg_broad['ewbalmer_broad'] = r'EW(H$\gamma)_{\mathrm{broad}}=$'+r'${:.1f}$'.format(fastspec['HGAMMA_BROAD_EW'])+r' $\AA$' + + if (fastspec['OII_3726_AMP']*np.sqrt(fastspec['OII_3726_AMP_IVAR']) > snrcut or + fastspec['OII_3729_AMP']*np.sqrt(fastspec['OII_3729_AMP_IVAR']) > snrcut): + leg_narrow['ewoii'] = r'EW([OII])'+r'$={:.1f}$'.format(fastspec['OII_3726_EW']+fastspec['OII_3729_EW'])+r' $\AA$' + + if fastspec['OIII_5007_AMP']*np.sqrt(fastspec['OIII_5007_AMP_IVAR']) > snrcut: + leg_narrow['ewoiii'] = r'EW([OIII])'+r'$={:.1f}$'.format(fastspec['OIII_5007_EW'])+r' $\AA$' + + # choose one Balmer line + if fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut: + leg_narrow['ewbalmer_narrow'] = r'EW(H$\alpha)=$'+r'${:.1f}$'.format(fastspec['HALPHA_EW'])+r' $\AA$' + elif fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut: + leg_narrow['ewbalmer_narrow'] = r'EW(H$\beta)=$'+r'${:.1f}$'.format(fastspec['HBETA_EW'])+r' $\AA$' + elif fastspec['HGAMMA_AMP']*np.sqrt(fastspec['HGAMMA_AMP_IVAR']) > snrcut: + leg_narrow['ewbalmer_narrow'] = r'EW(H$\gamma)=$'+r'${:.1f}$'.format(fastspec['HGAMMA_EW'])+r' $\AA$' + + if (fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut and + fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut): + leg_narrow['hahb'] = r'$\mathrm{H}\alpha/\mathrm{H}\beta=$'+r'${:.3f}$'.format(fastspec['HALPHA_FLUX']/fastspec['HBETA_FLUX']) + if 'hahb' not in leg_narrow.keys() and (fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut and + fastspec['HGAMMA_AMP']*np.sqrt(fastspec['HGAMMA_AMP_IVAR']) > snrcut): + leg_narrow['hbhg'] = r'$\mathrm{H}\beta/\mathrm{H}\gamma=$'+r'${:.3f}$'.format(fastspec['HBETA_FLUX']/fastspec['HGAMMA_FLUX']) + if (fastspec['HBETA_AMP']*np.sqrt(fastspec['HBETA_AMP_IVAR']) > snrcut and + fastspec['OIII_5007_AMP']*np.sqrt(fastspec['OIII_5007_AMP_IVAR']) > snrcut and + fastspec['HBETA_FLUX'] > 0 and fastspec['OIII_5007_FLUX'] > 0): + leg_narrow['oiiihb'] = r'$\log_{10}(\mathrm{[OIII]/H}\beta)=$'+r'${:.3f}$'.format( + np.log10(fastspec['OIII_5007_FLUX']/fastspec['HBETA_FLUX'])) + if (fastspec['HALPHA_AMP']*np.sqrt(fastspec['HALPHA_AMP_IVAR']) > snrcut and + fastspec['NII_6584_AMP']*np.sqrt(fastspec['NII_6584_AMP_IVAR']) > snrcut and + fastspec['HALPHA_FLUX'] > 0 and fastspec['NII_6584_FLUX'] > 0): + leg_narrow['niiha'] = r'$\log_{10}(\mathrm{[NII]/H}\alpha)=$'+r'${:.3f}$'.format( + np.log10(fastspec['NII_6584_FLUX']/fastspec['HALPHA_FLUX'])) + + if (fastspec['OII_3726_AMP']*np.sqrt(fastspec['OII_3726_AMP_IVAR']) > snrcut or + fastspec['OII_3729_AMP']*np.sqrt(fastspec['OII_3729_AMP_IVAR']) > snrcut): + #if fastspec['OII_DOUBLET_RATIO'] != 0: + leg_narrow['oii_doublet'] = r'[OII] $\lambda3726/\lambda3729={:.3f}$'.format(fastspec['OII_DOUBLET_RATIO']) + + if fastspec['SII_6716_AMP']*np.sqrt(fastspec['SII_6716_AMP_IVAR']) > snrcut or fastspec['SII_6731_AMP']*np.sqrt(fastspec['SII_6731_AMP_IVAR']) > snrcut: + #if fastspec['SII_DOUBLET_RATIO'] != 0: + leg_narrow['sii_doublet'] = r'[SII] $\lambda6731/\lambda6716={:.3f}$'.format(fastspec['SII_DOUBLET_RATIO']) + + # rebuild the best-fitting broadband photometric model + if not stackfit: + sedmodel, sedphot = CTools.templates2data( + templatecache['templateflux'], templatecache['templatewave'], + redshift=redshift, dluminosity=dlum, photsys=metadata['PHOTSYS'], + synthphot=True, coeff=fastspec['COEFF'] * CTools.massnorm) + sedwave = templatecache['templatewave'] * (1 + redshift) + + phot = CTools.parse_photometry(CTools.bands, + maggies=np.array([metadata['FLUX_{}'.format(band.upper())] for band in CTools.bands]), + ivarmaggies=np.array([metadata['FLUX_IVAR_{}'.format(band.upper())] for band in CTools.bands]), + lambda_eff=allfilters.effective_wavelengths.value, + min_uncertainty=CTools.min_uncertainty) + #if hasattr(CTools, 'fiber_bands'): + # fiberphot = CTools.parse_photometry(CTools.fiber_bands, + # maggies=np.array([metadata['FIBERTOTFLUX_{}'.format(band.upper())] + # for band in CTools.fiber_bands]), + # lambda_eff=filters.effective_wavelengths.value) + #else: + # fiberphot = None + + indx_phot = np.where((sedmodel > 0) * (sedwave/1e4 > phot_wavelims[0]) * + (sedwave/1e4 < phot_wavelims[1]))[0] + sedwave = sedwave[indx_phot] + sedmodel = sedmodel[indx_phot] + + if not fastphot: + # Rebuild the best-fitting spectroscopic model; prefix "desi" means + # "per-camera" and prefix "full" has the cameras h-stacked. + fullwave = np.hstack(data['wave']) + + desicontinuum, _ = CTools.templates2data(templatecache['templateflux_nolines'], templatecache['templatewave'], + redshift=redshift, dluminosity=dlum, synthphot=False, + specwave=data['wave'], specres=data['res'], + specmask=data['mask'], cameras=data['cameras'], + vdisp=fastspec['VDISP'], + coeff=fastspec['COEFF']) + + # remove the aperture correction + desicontinuum = [_desicontinuum / apercorr for _desicontinuum in desicontinuum] + fullcontinuum = np.hstack(desicontinuum) + + # Need to be careful we don't pass a large negative residual where + # there are gaps in the data. + desiresiduals = [] + for icam in np.arange(len(data['cameras'])): + resid = data['flux'][icam] - desicontinuum[icam] + I = (data['flux'][icam] == 0.0) * (data['flux'][icam] == 0.0) + if np.any(I): + resid[I] = 0.0 + desiresiduals.append(resid) + + if np.all(fastspec['COEFF'] == 0): + fullsmoothcontinuum = np.zeros_like(fullwave) + else: + fullsmoothcontinuum, _ = CTools.smooth_continuum( + fullwave, np.hstack(desiresiduals), np.hstack(data['ivar']), + redshift=redshift, linemask=np.hstack(data['linemask'])) + if no_smooth_continuum: + fullsmoothcontinuum *= 0 + + desismoothcontinuum = [] + for campix in data['camerapix']: + desismoothcontinuum.append(fullsmoothcontinuum[campix[0]:campix[1]]) + + # full model spectrum + individual line-spectra + desiemlines = EMFit.emlinemodel_bestfit(data['wave'], data['res'], fastspec, snrcut=emline_snrmin) + + desiemlines_oneline = [] + inrange = ( (EMFit.linetable['restwave'] * (1+redshift) > np.min(fullwave)) * + (EMFit.linetable['restwave'] * (1+redshift) < np.max(fullwave)) ) + + for oneline in EMFit.linetable[inrange]: # for all lines in range + linename = oneline['name'].upper() + amp = fastspec['{}_MODELAMP'.format(linename)] + if amp != 0: + desiemlines_oneline1 = build_emline_model( + EMFit.dlog10wave, redshift, np.array([amp]), + np.array([fastspec['{}_VSHIFT'.format(linename)]]), + np.array([fastspec['{}_SIGMA'.format(linename)]]), + np.array([oneline['restwave']]), data['wave'], data['res']) + desiemlines_oneline.append(desiemlines_oneline1) + + # Grab the viewer cutout. + if not stackfit: + pixscale = 0.262 + width = int(30 / pixscale) # =1 arcmin + height = int(width / 1.3) # 3:2 aspect ratio + + hdr = fits.Header() + hdr['NAXIS'] = 2 + hdr['NAXIS1'] = width + hdr['NAXIS2'] = height + hdr['CTYPE1'] = 'RA---TAN' + hdr['CTYPE2'] = 'DEC--TAN' + hdr['CRVAL1'] = metadata['RA'] + hdr['CRVAL2'] = metadata['DEC'] + hdr['CRPIX1'] = width/2+0.5 + hdr['CRPIX2'] = height/2+0.5 + hdr['CD1_1'] = -pixscale/3600 + hdr['CD1_2'] = 0.0 + hdr['CD2_1'] = 0.0 + hdr['CD2_2'] = +pixscale/3600 + wcs = WCS(hdr) + + cutoutjpeg = os.path.join(outdir, 'tmp.'+os.path.basename(pngfile.replace('.png', '.jpeg'))) + if not os.path.isfile(cutoutjpeg): + wait = 5 # seconds + #cmd = 'timeout {wait} wget -q -o /dev/null -O {outfile} https://www.legacysurvey.org/viewer/jpeg-cutout?ra={ra}&dec={dec}&width={width}&height={height}&layer={layer}' + #cmd = cmd.format(wait=wait, outfile=cutoutjpeg, ra=metadata['RA'], + # dec=metadata['DEC'], width=width, height=height, + # layer=layer) + cmd = 'wget -q -O {outfile} https://www.legacysurvey.org/viewer/jpeg-cutout?ra={ra}&dec={dec}&width={width}&height={height}&layer={layer}' + cmd = cmd.format(outfile=cutoutjpeg, ra=metadata['RA'], dec=metadata['DEC'], + width=width, height=height, layer=layer) + log.info(cmd) + try: + subprocess.check_output(cmd.split(), stderr=subprocess.DEVNULL, timeout=wait) + #subprocess.check_call(cmd.split(), stderr=subprocess.DEVNULL) + except: + log.warning('No cutout from viewer after {} seconds; stopping wget'.format(wait)) + try: + img = mpimg.imread(cutoutjpeg) + except: + log.warning('Problem reading cutout for targetid {}.'.format(metadata['TARGETID'])) + img = np.zeros((height, width, 3)) + + if os.path.isfile(cutoutjpeg): + os.remove(cutoutjpeg) + + # QA choices + legxpos, legypos, legypos2, legfntsz1, legfntsz = 0.98, 0.94, 0.05, 16, 18 + bbox = dict(boxstyle='round', facecolor='lightgray', alpha=0.15) + bbox2 = dict(boxstyle='round', facecolor='lightgray', alpha=0.7) + + if fastphot: + fullheight = 9 # inches + fullwidth = 18 + + nrows = 3 + ncols = 8 + + fig = plt.figure(figsize=(fullwidth, fullheight)) + gs = fig.add_gridspec(nrows, ncols)#, width_ratios=width_ratios) + + cutax = fig.add_subplot(gs[0:2, 5:8], projection=wcs) # rows x cols + sedax = fig.add_subplot(gs[0:3, 0:5]) + elif stackfit: + fullheight = 14 # inches + fullwidth = 24 + + # 8 columns: 3 for the spectra, and 5 for the lines + # 8 rows: 4 for the SED, 2 each for the spectra, 1 gap, and 3 for the lines + nlinerows = 6 + nlinecols = 4 + nrows = nlinerows + ncols = 9 + + #height_ratios = np.hstack(([1.0]*3, 0.25, [1.0]*6)) + #width_ratios = np.hstack(([1.0]*5, [1.0]*3)) + + fig = plt.figure(figsize=(fullwidth, fullheight)) + gs = fig.add_gridspec(nrows, ncols)#, height_ratios=height_ratios)#, width_ratios=width_ratios) + + specax = fig.add_subplot(gs[0:4, 0:5]) + else: + fullheight = 18 # inches + fullwidth = 24 + + # 8 columns: 3 for the SED, 5 for the spectra, and 8 for the lines + # 8 rows: 4 for the SED, 2 each for the spectra, 1 gap, and 3 for the lines + ngaprows = 1 + nlinerows = 6 + nlinecols = 3 + nrows = 9 + ngaprows + ncols = 8 + + height_ratios = np.hstack(([1.0]*3, 0.25, [1.0]*6)) # small gap + #width_ratios = np.hstack(([1.0]*5, [1.0]*3)) + + fig = plt.figure(figsize=(fullwidth, fullheight)) + gs = fig.add_gridspec(nrows, ncols, height_ratios=height_ratios)#, width_ratios=width_ratios) + + cutax = fig.add_subplot(gs[0:3, 5:8], projection=wcs) # rows x cols + sedax = fig.add_subplot(gs[0:3, 0:5]) + specax = fig.add_subplot(gs[4:8, 0:5]) + + # viewer cutout + if not stackfit: + cutax.imshow(img, origin='lower')#, interpolation='nearest') + cutax.set_xlabel('RA [J2000]') + cutax.set_ylabel('Dec [J2000]') + + cutax.coords[1].set_ticks_position('r') + cutax.coords[1].set_ticklabel_position('r') + cutax.coords[1].set_axislabel_position('r') + + if metadata['DEC'] > 0: + sgn = '+' + else: + sgn = '' + + cutax.text(0.04, 0.95, '$(\\alpha,\\delta)$=({:.7f}, {}{:.6f})'.format(metadata['RA'], sgn, metadata['DEC']), + ha='left', va='top', color='k', fontsize=fontsize1, bbox=bbox2, + transform=cutax.transAxes) + + sz = img.shape + cutax.add_artist(Circle((sz[1] / 2, sz[0] / 2), radius=1.5/2/pixscale, facecolor='none', # DESI fiber=1.5 arcsec diameter + edgecolor='red', ls='-', alpha=0.8))#, label='3" diameter')) + cutax.add_artist(Circle((sz[1] / 2, sz[0] / 2), radius=10/2/pixscale, facecolor='none', + edgecolor='red', ls='--', alpha=0.8))#, label='15" diameter')) + handles = [Line2D([0], [0], color='red', lw=2, ls='-', label='1.5 arcsec'), + Line2D([0], [0], color='red', lw=2, ls='--', label='10 arcsec')] + + cutax.legend(handles=handles, loc='lower left', fontsize=fontsize1, facecolor='lightgray') + + if not fastphot: + # plot the full spectrum + best-fitting (total) model + specax.plot(fullwave/1e4, fullsmoothcontinuum, color='gray', alpha=0.4) + specax.plot(fullwave/1e4, fullcontinuum, color='k', alpha=0.6) + + spec_ymin, spec_ymax = 1e6, -1e6 + + desimodelspec = [] + for icam in np.arange(len(data['cameras'])): # iterate over cameras + wave = data['wave'][icam] + flux = data['flux'][icam] + modelflux = desiemlines[icam] + desicontinuum[icam] + desismoothcontinuum[icam] + + sigma, camgood = ivar2var(data['ivar'][icam], sigma=True, allmasked_ok=True, clip=0) + + wave = wave[camgood] + flux = flux[camgood] + sigma = sigma[camgood] + modelflux = modelflux[camgood] + + desimodelspec.append(apercorr * (desicontinuum[icam] + desiemlines[icam])) + + # get the robust range + filtflux = median_filter(flux, 51, mode='nearest') + if np.sum(camgood) > 0: + sigflux = np.diff(np.percentile(flux - modelflux, [25, 75]))[0] / 1.349 # robust + if -2 * sigflux < spec_ymin: + spec_ymin = -2 * sigflux + if 6 * sigflux > spec_ymax: + spec_ymax = 6 * sigflux + if np.max(filtflux) > spec_ymax: + #print(icam, spec_ymax, np.max(filtflux), np.max(filtflux) * 1.2) + spec_ymax = np.max(filtflux) * 1.25 + if np.max(modelflux) > spec_ymax: + spec_ymax = np.max(modelflux) * 1.25 + #print(spec_ymin, spec_ymax) + #pdb.set_trace() + + #specax.fill_between(wave, flux-sigma, flux+sigma, color=col1[icam], alpha=0.2) + if nsmoothspec > 1: + from scipy.ndimage import gaussian_filter + specax.plot(wave/1e4, gaussian_filter(flux, nsmoothspec), color=col1[icam], alpha=0.8) + specax.plot(wave/1e4, gaussian_filter(modelflux, nsmoothspec), color=col2[icam], lw=2, alpha=0.8) + else: + specax.plot(wave/1e4, flux, color=col1[icam], alpha=0.8) + specax.plot(wave/1e4, modelflux, color=col2[icam], lw=2, alpha=0.8) + + fullmodelspec = np.hstack(desimodelspec) + + if stackfit: + specax_twin = specax.twiny() + specax_twin.set_xlim(spec_wavelims[0]/(1+redshift)/1e4, spec_wavelims[1]/(1+redshift)/1e4) + specax_twin.xaxis.set_major_formatter(major_formatter) + restticks = np.array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1]) + restticks = restticks[(restticks >= spec_wavelims[0]/(1+redshift)/1e4) * (restticks <= spec_wavelims[1]/(1+redshift)/1e4)] + specax_twin.set_xticks(restticks) + else: + specax.spines[['top']].set_visible(False) + + specax.set_xlim(spec_wavelims[0]/1e4, spec_wavelims[1]/1e4) + specax.set_ylim(spec_ymin, spec_ymax) + specax.set_xlabel(r'Observed-frame Wavelength ($\mu$m)') + #specax.set_xlabel(r'Observed-frame Wavelength ($\AA$)') + specax.set_ylabel(r'$F_{\lambda}\ (10^{-17}~{\rm erg}~{\rm s}^{-1}~{\rm cm}^{-2}~\AA^{-1})$') + + # photometric SED + if not stackfit: + abmag_good = phot['abmag_ivar'] > 0 + abmag_goodlim = phot['abmag_limit'] > 0 + + if len(sedmodel) == 0: + log.warning('Best-fitting photometric continuum is all zeros or negative!') + if np.sum(abmag_good) > 0: + medmag = np.median(phot['abmag'][abmag_good]) + elif np.sum(abmag_goodlim) > 0: + medmag = np.median(phot['abmag_limit'][abmag_goodlim]) + else: + medmag = 0.0 + sedmodel_abmag = np.zeros_like(templatecache['templatewave']) + medmag + else: + factor = 10**(0.4 * 48.6) * sedwave**2 / (C_LIGHT * 1e13) / FLUXNORM / CTools.massnorm # [erg/s/cm2/A --> maggies] + sedmodel_abmag = -2.5*np.log10(sedmodel * factor) + sedax.plot(sedwave / 1e4, sedmodel_abmag, color='slategrey', alpha=0.9, zorder=1) + + sedax.scatter(sedphot['lambda_eff']/1e4, sedphot['abmag'], marker='s', + s=400, color='k', facecolor='none', linewidth=2, alpha=1.0, zorder=2) + + if not fastphot: + #factor = 10**(0.4 * 48.6) * fullwave**2 / (C_LIGHT * 1e13) / FLUXNORM # [erg/s/cm2/A --> maggies] + #good = fullmodelspec > 0 + #sedax.plot(fullwave[good]/1e4, -2.5*np.log10(fullmodelspec[good]*factor[good]), color='purple', alpha=0.8) + for icam in np.arange(len(data['cameras'])): + factor = 10**(0.4 * 48.6) * data['wave'][icam]**2 / (C_LIGHT * 1e13) / FLUXNORM # [erg/s/cm2/A --> maggies] + good = desimodelspec[icam] > 0 + _wave = data['wave'][icam][good]/1e4 + _flux = -2.5*np.log10(desimodelspec[icam][good]*factor[good]) + sedax.plot(_wave, _flux, color=col2[icam], alpha=0.8) + #pdb.set_trace() + + # we have to set the limits *before* we call errorbar, below! + dm = 1.5 + sed_ymin = np.nanmax(sedmodel_abmag) + dm + sed_ymax = np.nanmin(sedmodel_abmag) - dm + if np.sum(abmag_good) > 0 and np.sum(abmag_goodlim) > 0: + sed_ymin = np.max((np.nanmax(phot['abmag'][abmag_good]), np.nanmax(phot['abmag_limit'][abmag_goodlim]), np.nanmax(sedmodel_abmag))) + dm + sed_ymax = np.min((np.nanmin(phot['abmag'][abmag_good]), np.nanmin(phot['abmag_limit'][abmag_goodlim]), np.nanmin(sedmodel_abmag))) - dm + elif np.sum(abmag_good) > 0 and np.sum(abmag_goodlim) == 0: + sed_ymin = np.max((np.nanmax(phot['abmag'][abmag_good]), np.nanmax(sedmodel_abmag))) + dm + sed_ymax = np.min((np.nanmin(phot['abmag'][abmag_good]), np.nanmin(sedmodel_abmag))) - dm + elif np.sum(abmag_good) == 0 and np.sum(abmag_goodlim) > 0: + sed_ymin = np.max((np.nanmax(phot['abmag_limit'][abmag_goodlim]), np.nanmax(sedmodel_abmag))) + dm + sed_ymax = np.min((np.nanmin(phot['abmag_limit'][abmag_goodlim]), np.nanmin(sedmodel_abmag))) - dm + else: + abmag_good = phot['abmag'] > 0 + abmag_goodlim = phot['abmag_limit'] > 0 + if np.sum(abmag_good) > 0 and np.sum(abmag_goodlim) > 0: + sed_ymin = np.max((np.nanmax(phot['abmag'][abmag_good]), np.nanmax(phot['abmag_limit'][abmag_goodlim]))) + dm + sed_ymax = np.min((np.nanmin(phot['abmag'][abmag_good]), np.nanmin(phot['abmag_limit'][abmag_goodlim]))) - dm + elif np.sum(abmag_good) > 0 and np.sum(abmag_goodlim) == 0: + sed_ymin = np.nanmax(phot['abmag'][abmag_good]) + dm + sed_ymax = np.nanmin(phot['abmag'][abmag_good]) - dm + elif np.sum(abmag_good) == 0 and np.sum(abmag_goodlim) > 0: + sed_ymin = np.nanmax(phot['abmag_limit'][abmag_goodlim]) + dm + sed_ymax = np.nanmin(phot['abmag_limit'][abmag_goodlim]) - dm + #else: + # sed_ymin, sed_ymax = [30, 20] + + if sed_ymin > 30: + sed_ymin = 30 + if np.isnan(sed_ymin) or np.isnan(sed_ymax): + raise('Problem here!') + #print(sed_ymin, sed_ymax) + + #sedax.set_xlabel(r'Observed-frame Wavelength ($\mu$m)') + sedax.set_xlim(phot_wavelims[0], phot_wavelims[1]) + sedax.set_xscale('log') + sedax.set_ylabel('AB mag') + #sedax.set_ylabel(r'Apparent Brightness (AB mag)') + sedax.set_ylim(sed_ymin, sed_ymax) + + sedax.xaxis.set_major_formatter(major_formatter) + obsticks = np.array([0.1, 0.2, 0.5, 1.0, 1.5, 3.0, 5.0, 10.0, 20.0]) + obsticks = obsticks[(obsticks >= phot_wavelims[0]) * (obsticks <= phot_wavelims[1])] + sedax.set_xticks(obsticks) + + if fastphot: + sedax.set_xlabel(r'Observed-frame Wavelength ($\mu$m)') + + sedax_twin = sedax.twiny() + sedax_twin.set_xlim(phot_wavelims[0]/(1+redshift), phot_wavelims[1]/(1+redshift)) + sedax_twin.set_xscale('log') + #sedax_twin.set_xlabel(r'Rest-frame Wavelength ($\mu$m)') + + sedax_twin.xaxis.set_major_formatter(major_formatter) + restticks = np.array([0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 1.5, 3.0, 5.0, 10.0, 15.0, 20.0]) + restticks = restticks[(restticks >= phot_wavelims[0]/(1+redshift)) * (restticks <= phot_wavelims[1]/(1+redshift))] + sedax_twin.set_xticks(restticks) + + # integrated flux / photometry + abmag = np.squeeze(phot['abmag']) + abmag_limit = np.squeeze(phot['abmag_limit']) + abmag_fainterr = np.squeeze(phot['abmag_fainterr']) + abmag_brighterr = np.squeeze(phot['abmag_brighterr']) + yerr = np.squeeze([abmag_brighterr, abmag_fainterr]) + + markersize = 14 + + dofit = np.where(CTools.bands_to_fit)[0] + if len(dofit) > 0: + good = np.where((abmag[dofit] > 0) * (abmag_limit[dofit] == 0))[0] + upper = np.where(abmag_limit[dofit] > 0)[0] + if len(good) > 0: + sedax.errorbar(phot['lambda_eff'][dofit][good]/1e4, abmag[dofit][good], + yerr=yerr[:, dofit[good]], + fmt='o', markersize=markersize, markeredgewidth=1, markeredgecolor='k', + markerfacecolor=photcol1, elinewidth=3, ecolor=photcol1, capsize=4, + label=r'$grz\,W_{1}W_{2}W_{3}W_{4}$', zorder=2, alpha=1.0) + if len(upper) > 0: + sedax.errorbar(phot['lambda_eff'][dofit][upper]/1e4, abmag_limit[dofit][upper], + lolims=True, yerr=0.75, + fmt='o', markersize=markersize, markeredgewidth=3, markeredgecolor='k', + markerfacecolor=photcol1, elinewidth=3, ecolor=photcol1, capsize=4, alpha=0.7) + + ignorefit = np.where(CTools.bands_to_fit == False)[0] + if len(ignorefit) > 0: + good = np.where((abmag[ignorefit] > 0) * (abmag_limit[ignorefit] == 0))[0] + upper = np.where(abmag_limit[ignorefit] > 0)[0] + if len(good) > 0: + sedax.errorbar(phot['lambda_eff'][ignorefit][good]/1e4, abmag[ignorefit][good], + yerr=yerr[:, ignorefit[good]], + fmt='o', markersize=markersize, markeredgewidth=3, markeredgecolor='k', + markerfacecolor='none', elinewidth=3, ecolor=photcol1, capsize=4, alpha=0.7) + if len(upper) > 0: + sedax.errorbar(phot['lambda_eff'][ignorefit][upper]/1e4, abmag_limit[ignorefit][upper], + lolims=True, yerr=0.75, fmt='o', markersize=markersize, markeredgewidth=3, + markeredgecolor='k', markerfacecolor='none', elinewidth=3, + ecolor=photcol1, capsize=5, alpha=0.7) + + if fastphot: + txt = leg['rchi2_phot'] + else: + # Label the DESI wavelength range and the aperture correction. + sedax.plot([np.min(fullwave)/1e4, np.max(fullwave)/1e4], [sed_ymin-1, sed_ymin-1], + lw=2, ls='-', color='gray', marker='s')#, alpha=0.5) + sedax.text(((np.max(fullwave)-np.min(fullwave))/2+np.min(fullwave)*0.8)/1e4, sed_ymin-1.7, + 'DESI x {:.2f}'.format(apercorr), ha='center', va='center', fontsize=16, + color='k') + + txt = '\n'.join((leg['rchi2_cont'], leg['rchi2_phot'], leg['rchi2'])) + + sedax.text(0.02, 0.94, txt, ha='left', va='top', + transform=sedax.transAxes, fontsize=legfntsz)#, bbox=bbox) + + txt = '\n'.join(( + #r'{}'.format(leg['fagn']), + r'{}'.format(leg['zzsun']), + r'{}'.format(leg['AV']), + r'{}'.format(leg['sfr']), + r'{}'.format(leg['age']), + r'{}'.format(leg['mstar']), + )) + sedax.text(legxpos, legypos2, txt, ha='right', va='bottom', + transform=sedax.transAxes, fontsize=legfntsz1, bbox=bbox) + + if not fastphot: + # draw lines connecting the SED and spectral plots + sedax.add_artist(ConnectionPatch(xyA=(spec_wavelims[0]/1e4, sed_ymin), + xyB=(spec_wavelims[0]/1e4, spec_ymax), + coordsA='data', coordsB='data', + axesA=sedax, axesB=specax, color='k')) + sedax.add_artist(ConnectionPatch(xyA=(spec_wavelims[1]/1e4, sed_ymin), + xyB=(spec_wavelims[1]/1e4, spec_ymax), + coordsA='data', coordsB='data', + axesA=sedax, axesB=specax, color='k')) + + # zoom in on individual emission lines - use linetable! + if not fastphot: + linetable = EMFit.linetable + inrange = (linetable['restwave'] * (1+redshift) > np.min(fullwave)) * (linetable['restwave'] * (1+redshift) < np.max(fullwave)) + linetable = linetable[inrange] + + nline = len(set(linetable['plotgroup'])) + + plotsig_default = 200.0 # [km/s] + plotsig_default_balmer = 500.0 # [km/s] + plotsig_default_broad = 2000.0 # [km/s] + + minwaves, maxwaves, meanwaves, deltawaves, sigmas, linenames = [], [], [], [], [], [] + for plotgroup in set(linetable['plotgroup']): + I = np.where(plotgroup == linetable['plotgroup'])[0] + linename = linetable['nicename'][I[0]].replace('-', ' ') + linenames.append(linename) + minwaves.append(np.min(linetable['restwave'][I])) + maxwaves.append(np.max(linetable['restwave'][I])) + meanwaves.append(np.mean(linetable['restwave'][I])) + deltawaves.append((np.max(linetable['restwave'][I]) - np.min(linetable['restwave'][I])) / 2) + + sigmas1 = np.array([fastspec['{}_SIGMA'.format(line.upper())] for line in linetable[I]['name']]) + sigmas1 = sigmas1[sigmas1 > 0] + if len(sigmas1) > 0: + plotsig = 1.5*np.mean(sigmas1) + if plotsig < 50: + plotsig = 50.0 + else: + if np.any(linetable['isbroad'][I]): + if np.any(linetable['isbalmer'][I]): + plotsig = fastspec['BROAD_SIGMA'] + if plotsig < 50: + plotsig = fastspec['NARROW_SIGMA'] + if plotsig < 50: + plotsig = plotsig_default + #plotsig = plotsig_default_broad + else: + plotsig = fastspec['UV_SIGMA'] + if plotsig < 50: + plotsig = plotsig_default_broad + else: + plotsig = fastspec['NARROW_SIGMA'] + if plotsig < 50: + plotsig = plotsig_default + sigmas.append(plotsig) + + if len(linetable) == 0: + ax = [] + else: + srt = np.argsort(meanwaves) + minwaves = np.hstack(minwaves)[srt] + maxwaves = np.hstack(maxwaves)[srt] + meanwaves = np.hstack(meanwaves)[srt] + deltawaves = np.hstack(deltawaves)[srt] + sigmas = np.hstack(sigmas)[srt] + linenames = np.hstack(linenames)[srt] + + # Add the linenames to the spectrum plot. + for meanwave, linename in zip(meanwaves*(1+redshift), linenames): + #print(meanwave, ymax_spec) + if meanwave > spec_wavelims[0] and meanwave < spec_wavelims[1]: + if 'SiIII' in linename: + thislinename = '\n'+linename.replace('+', '+\n ') + elif '4363' in linename: + thislinename = linename+'\n' + else: + thislinename = linename + if stackfit: + specax.text(meanwave/1e4, spec_ymax*0.97, thislinename, ha='center', va='top', + rotation=270, fontsize=12, alpha=0.5) + else: + specax.text(meanwave/1e4, spec_ymax, thislinename, ha='center', va='top', + rotation=270, fontsize=12, alpha=0.5) + + removelabels = np.ones(nline, bool) + line_ymin, line_ymax = np.zeros(nline)+1e6, np.zeros(nline)-1e6 + + if stackfit: + ax, irow, colshift = [], 0, 5 + else: + ax, irow, colshift = [], 4, 5 # skip the gap row + + for iax, (minwave, maxwave, meanwave, deltawave, sig, linename) in enumerate( + zip(minwaves, maxwaves, meanwaves, deltawaves, sigmas, linenames)): + icol = iax % nlinecols + icol += colshift + if iax > 0 and iax % nlinecols == 0: + irow += 1 + #print(iax, irow, icol) + + xx = fig.add_subplot(gs[irow, icol]) + ax.append(xx) + + wmin = (minwave - deltawave) * (1+redshift) - 6 * sig * minwave * (1+redshift) / C_LIGHT + wmax = (maxwave + deltawave) * (1+redshift) + 6 * sig * maxwave * (1+redshift) / C_LIGHT + #wmin = (meanwave - deltawave) * (1+redshift) - 6 * sig * meanwave * (1+redshift) / C_LIGHT + #wmax = (meanwave + deltawave) * (1+redshift) + 6 * sig * meanwave * (1+redshift) / C_LIGHT + + # iterate over cameras + for icam in np.arange(len(data['cameras'])): # iterate over cameras + emlinewave = data['wave'][icam] + emlineflux = data['flux'][icam] - desicontinuum[icam] - desismoothcontinuum[icam] + emlinemodel = desiemlines[icam] + + emlinesigma, good = ivar2var(data['ivar'][icam], sigma=True, allmasked_ok=True, clip=0) + emlinewave = emlinewave[good] + emlineflux = emlineflux[good] + emlinesigma = emlinesigma[good] + emlinemodel = emlinemodel[good] + + emlinemodel_oneline = [] + for desiemlines_oneline1 in desiemlines_oneline: + emlinemodel_oneline.append(desiemlines_oneline1[icam][good]) + + indx = np.where((emlinewave > wmin) * (emlinewave < wmax))[0] + if len(indx) > 1: + removelabels[iax] = False + xx.plot(emlinewave[indx]/1e4, emlineflux[indx], color=col1[icam], alpha=0.5) + #if nsmoothspec > 1: + # xx.plot(emlinewave[indx]/1e4, gaussian_filter(emlineflux[indx], nsmoothspec), color=col1[icam], alpha=0.5) + #else: + # xx.plot(emlinewave[indx]/1e4, emlineflux[indx], color=col1[icam], alpha=0.5) + for emlinemodel_oneline1 in emlinemodel_oneline: + if np.sum(emlinemodel_oneline1[indx]) > 0: + xx.plot(emlinewave[indx]/1e4, emlinemodel_oneline1[indx], lw=1, alpha=0.8, color=col2[icam]) + #if nsmoothspec > 1: + # xx.plot(emlinewave[indx]/1e4, gaussian_filter(emlinemodel_oneline1[indx], nsmoothspec), lw=1, alpha=0.8, color=col2[icam]) + #else: + # xx.plot(emlinewave[indx]/1e4, emlinemodel_oneline1[indx], lw=1, alpha=0.8, color=col2[icam]) + xx.plot(emlinewave[indx]/1e4, emlinemodel[indx], color=col2[icam], lw=2, alpha=0.8) + #if nsmoothspec > 1: + # xx.plot(emlinewave[indx]/1e4, emlinemodel[indx], color=col2[icam], lw=2, alpha=0.8) + #else: + # xx.plot(emlinewave[indx]/1e4, gaussian_filter(emlinemodel[indx], nsmoothspec), color=col2[icam], lw=2, alpha=0.8) + + #xx.plot(emlinewave[indx], emlineflux[indx]-emlinemodel[indx], color='gray', alpha=0.3) + #xx.axhline(y=0, color='gray', ls='--') + + # get the robust range + sigflux = np.std(emlineflux[indx]) + filtflux = median_filter(emlineflux[indx], 3, mode='nearest') + + #_line_ymin, _line_ymax = -1.5 * sigflux, 4 * sigflux + #if np.max(emlinemodel[indx]) > _line_ymax: + # _line_ymax = np.max(emlinemodel[indx]) * 1.3 + _line_ymin, _line_ymax = -1.5 * sigflux, np.max(emlinemodel[indx]) * 1.4 + if 4 * sigflux > _line_ymax: + _line_ymax = 4 * sigflux + if np.max(filtflux) > _line_ymax: + _line_ymax = np.max(filtflux) + if np.min(emlinemodel[indx]) < _line_ymin: + _line_ymin = 0.8 * np.min(emlinemodel[indx]) + if _line_ymax > line_ymax[iax]: + line_ymax[iax] = _line_ymax + if _line_ymin < line_ymin[iax]: + line_ymin[iax] = _line_ymin + #print(linename, line_ymin[iax], line_ymax[iax]) + #if linename == '[OII] $\lambda\lambda$3726,29': + # pdb.set_trace() + + xx.set_xlim(wmin/1e4, wmax/1e4) + + xx.text(0.03, 0.89, linename, ha='left', va='center', + transform=xx.transAxes, fontsize=12) + xx.tick_params(axis='x', labelsize=16) + xx.tick_params(axis='y', labelsize=16) + + for iax, xx in enumerate(ax): + if removelabels[iax]: + xx.set_ylim(0, 1) + xx.set_xticklabels([]) + xx.set_yticklabels([]) + else: + xx.set_yticklabels([]) + xx.set_ylim(line_ymin[iax], line_ymax[iax]) + xx_twin = xx.twinx() + xx_twin.set_ylim(line_ymin[iax], line_ymax[iax]) + xlim = xx.get_xlim() + xx.xaxis.set_major_locator(ticker.MaxNLocator(2)) + + if fastphot: + plt.subplots_adjust(wspace=0.6, top=0.85, bottom=0.13, left=0.07, right=0.86) + + else: + plt.subplots_adjust(wspace=0.4, top=0.9, bottom=0.1, left=0.07, right=0.92, hspace=0.33) + + # common axis labels + if len(ax) > 0: + if len(ax) == 2: + xx = fig.add_subplot(gs[irow, icol+1]) + xx.axis('off') + ax.append(xx) + elif len(ax) == 1: + xx = fig.add_subplot(gs[irow, icol+1]) + xx.axis('off') + ax.append(xx) + xx = fig.add_subplot(gs[irow, icol+2]) + xx.axis('off') + ax.append(xx) + + ulpos = ax[0].get_position() + lpos = ax[nline-1].get_position() + urpos = ax[2].get_position() + xpos = (urpos.x1 - ulpos.x0) / 2 + ulpos.x0# + 0.03 + + ypos = lpos.y0 - 0.04 + fig.text(xpos, ypos, r'Observed-frame Wavelength ($\mu$m)', + ha='center', va='center', fontsize=fontsize2) + + xpos = urpos.x1 + 0.05 + ypos = (urpos.y1 - lpos.y0) / 2 + lpos.y0# + 0.03 + fig.text(xpos, ypos, + r'$F_{\lambda}\ (10^{-17}~{\rm erg}~{\rm s}^{-1}~{\rm cm}^{-2}~\AA^{-1})$', + ha='center', va='center', rotation=270, fontsize=fontsize2) + + legkeys = leg.keys() + legfntsz = 17 + + # rest wavelength plus targeting information + if fastphot: + ppos = sedax.get_position() + xpos = (ppos.x1 - ppos.x0) / 2 + ppos.x0 + ypos = ppos.y1 + 0.07 + + fig.text(xpos, ypos, r'Rest-frame Wavelength ($\mu$m)', + ha='center', va='bottom', fontsize=fontsize2) + fig.text(0.58, 0.87, '\n'.join(target), ha='left', va='bottom', + fontsize=fontsize2, linespacing=1.4) + + toppos, leftpos = 0.25, 0.59 + + txt = [ + r'{}'.format(leg['z']), + ] + if 'zredrock' in legkeys: + txt += [r'{}'.format(leg['zredrock'])] + txt += [ + #r'{}'.format(leg['zwarn']), + #'', + r'{}'.format(leg['vdisp']), + '', + r'{}'.format(leg['dn4000_model']), + ] + + fig.text(leftpos, toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, + bbox=bbox, linespacing=1.4) + + txt = [r'{}'.format(leg['absmag_r'])] + if 'absmag_gr' in legkeys: + txt += [r'{}'.format(leg['absmag_gr'])] + if 'absmag_rz' in legkeys: + txt += [r'{}'.format(leg['absmag_rz'])] + + fig.text(leftpos+0.18, toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, + bbox=bbox, linespacing=1.4) + + else: + if stackfit: + ppos = specax.get_position() + xpos = (ppos.x1 - ppos.x0) / 2 + ppos.x0 + ypos = ppos.y1 + 0.05 + + fig.text(xpos, ypos, r'Rest-frame Wavelength ($\mu$m)', + ha='center', va='bottom', fontsize=fontsize2) + fig.text(0.62, 0.91, '\n'.join(target), ha='left', va='bottom', + fontsize=fontsize2, linespacing=1.4) + else: + ppos = sedax.get_position() + xpos = (ppos.x1 - ppos.x0) / 2 + ppos.x0 + ypos = ppos.y1 + 0.03 + + fig.text(xpos, ypos, r'Rest-frame Wavelength ($\mu$m)', + ha='center', va='bottom', fontsize=fontsize2) + fig.text(0.647, 0.925, '\n'.join(target), ha='left', va='bottom', + fontsize=fontsize2, linespacing=1.4) + + # add some key results about the object at the bottom of the figure + + if stackfit: + #toppos, startpos, deltapos = 0.21, 0.04, 0.13 + toppos, leftpos, rightpos, adjust = 0.27, 0.05, 0.62, 0.01 + else: + toppos, leftpos, rightpos, adjust = 0.21, 0.03, 0.62, 0.01 + + nbox = 2 + 1*bool(leg_narrow) + bool(leg_broad) + boxpos = np.arange(nbox) * (rightpos - leftpos)/nbox + leftpos + + txt = [ + r'{}'.format(leg['z']), + ] + if 'zredrock' in legkeys: + txt += [r'{}'.format(leg['zredrock'])] + + txt += [ + #r'{}'.format(leg['zwarn']), + r'{}'.format(leg['vdisp']), + ] + if 'dv_narrow' in legkeys or 'dv_uv' in leg_uv.keys() or 'dv_broad' in leg_broad.keys(): + txt += [''] + + if 'dv_narrow' in legkeys: + txt += [ + r'{}'.format(leg['sigma_narrow']), + r'{}'.format(leg['dv_narrow']), + ] + if 'dv_uv' in leg_uv.keys(): + txt += [ + r'{}'.format(leg_uv['sigma_uv']), + r'{}'.format(leg_uv['dv_uv']), + ] + _ = leg_uv.pop('sigma_uv') + _ = leg_uv.pop('dv_uv') + if 'dv_broad' in leg_broad.keys(): + txt += [ + r'{}'.format(leg_broad['sigma_broad']), + r'{}'.format(leg_broad['dv_broad']), + ] + _ = leg_broad.pop('sigma_broad') + _ = leg_broad.pop('dv_broad') + + ibox = 0 + #fig.text(startpos, toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, + fig.text(boxpos[ibox], toppos, '\n'.join(txt), ha='left', va='top', fontsize=legfntsz, + bbox=bbox, linespacing=1.4) + ibox += 1 + + txt = [ + r'{}'.format(leg['absmag_r']), + r'{}'.format(leg['absmag_gr']), + r'{}'.format(leg['absmag_rz']), + '', + r'{}'.format(leg['dn4000_model']) + ] + if 'dn4000_spec' in legkeys: + txt += [r'{}'.format(leg['dn4000_spec'])] + + #fig.text(startpos+deltapos, toppos, '\n'.join(txt), ha='left', va='top', + fig.text(boxpos[ibox], toppos, '\n'.join(txt), ha='left', va='top', + fontsize=legfntsz, bbox=bbox, linespacing=1.4) + ibox += 1 + + #factor = 2 + if bool(leg_narrow): + txt = [] + for key in leg_narrow.keys(): + txt += [r'{}'.format(leg_narrow[key])] + #fig.text(startpos+deltapos*factor, toppos, '\n'.join(txt), ha='left', va='top', + fig.text(boxpos[ibox]-adjust*2, toppos, '\n'.join(txt), ha='left', va='top', + fontsize=legfntsz, bbox=bbox, linespacing=1.4) + ibox += 1 + #factor += 1.25 + + if bool(leg_broad): + txt = [] + for key in leg_broad.keys(): + txt += [r'{}'.format(leg_broad[key])] + + if bool(leg_uv): + if bool(leg_broad): + txt += [''] + else: + txt = [] + for key in leg_uv.keys(): + txt += [r'{}'.format(leg_uv[key])] + + if bool(leg_uv) or bool(leg_broad): + #fig.text(startpos+deltapos*factor, toppos, '\n'.join(txt), ha='left', va='top', + fig.text(boxpos[ibox]-adjust*1, toppos, '\n'.join(txt), ha='left', va='top', + fontsize=legfntsz, bbox=bbox, linespacing=1.4) + + log.info('Writing {}'.format(pngfile)) + fig.savefig(pngfile)#, dpi=150) + plt.close() + +def parse(options=None): + """Parse input arguments. + + """ + import sys, argparse + + parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter) + + parser.add_argument('--healpix', default=None, type=str, nargs='*', help="""Generate QA for all objects + with this healpixels (only defined for coadd-type 'healpix').""") + parser.add_argument('--tile', default=None, type=str, nargs='*', help='Generate QA for all objects on this tile.') + parser.add_argument('--night', default=None, type=str, nargs='*', help="""Generate QA for all objects observed on this + night (only defined for coadd-type 'pernight' and 'perexp').""") + parser.add_argument('--redux_dir', default=None, type=str, help='Top-level path to the reduced spectra.') + parser.add_argument('--redrockfiles', nargs='*', help='Optional full path to redrock file(s).') + parser.add_argument('--redrockfile-prefix', type=str, default='redrock-', help='Prefix of the input Redrock file name(s).') + parser.add_argument('--specfile-prefix', type=str, default='coadd-', help='Prefix of the spectral file(s).') + parser.add_argument('--qnfile-prefix', type=str, default='qso_qn-', help='Prefix of the QuasarNet afterburner file(s).') + parser.add_argument('--mapdir', type=str, default=None, help='Optional directory name for the dust maps.') + parser.add_argument('--fphotodir', type=str, default=None, help='Top-level location of the source photometry.') + parser.add_argument('--fphotoinfo', type=str, default=None, help='Photometric information file.') + + parser.add_argument('--emline_snrmin', type=float, default=0.0, help='Minimum emission-line S/N to be displayed.') + parser.add_argument('--nsmoothspec', type=int, default=0, help='Smoothing pixel value.') + + parser.add_argument('--minspecwave', type=float, default=3500., help='Minimum spectral wavelength (Angstrom).') + parser.add_argument('--maxspecwave', type=float, default=9900., help='Maximum spectral wavelength (Angstrom).') + parser.add_argument('--minphotwave', type=float, default=0.1, help='Minimum photometric wavelength (micron).') + parser.add_argument('--maxphotwave', type=float, default=35., help='Maximum photometric wavelength (micron).') + + parser.add_argument('--targetids', type=str, default=None, help='Comma-separated list of target IDs to process.') + parser.add_argument('-n', '--ntargets', type=int, help='Number of targets to process in each file.') + parser.add_argument('--firsttarget', type=int, default=0, help='Index of first object to to process in each file (0-indexed).') + parser.add_argument('--mp', type=int, default=1, help='Number of multiprocessing processes per MPI rank or node.') + parser.add_argument('--nophoto', action='store_true', help='Do not include the photometry in the model fitting.') + parser.add_argument('--overwrite', action='store_true', help='Overwrite existing files.') + + parser.add_argument('--imf', type=str, default='chabrier', help='Initial mass function.') + parser.add_argument('--templateversion', type=str, default='1.0.0', help='Template version number.') + parser.add_argument('--templates', type=str, default=None, help='Optional full path and filename to the templates.') + + parser.add_argument('--outprefix', default=None, type=str, help='Optional prefix for output filename.') + parser.add_argument('-o', '--outdir', default='.', type=str, help='Full path to desired output directory.') + + parser.add_argument('fastfitfile', nargs=1, help='Full path to fastspec or fastphot fitting results.') + + if options is None: + args = parser.parse_args() + log.info(' '.join(sys.argv)) + else: + args = parser.parse_args(options) + log.info('fastspecfit-qa {}'.format(' '.join(options))) + + return args + +def fastqa(args=None, comm=None): + """Main module. + + """ + import fitsio + from astropy.table import Table + from fastspecfit.io import (DESISpectra, get_templates_filename, get_qa_filename, + read_fastspecfit, DESI_ROOT_NERSC, select) + + if isinstance(args, (list, tuple, type(None))): + args = parse(args) + + if args.redux_dir is None: + args.redux_dir = os.path.join(os.environ.get('DESI_ROOT', DESI_ROOT_NERSC), 'spectro', 'redux') + if not os.path.isdir(args.redux_dir): + errmsg = 'Data reduction directory {} not found.'.format(args.redux_dir) + log.critical(errmsg) + raise IOError(errmsg) + + # Read the fitting results. + if not os.path.isfile(args.fastfitfile[0]): + log.warning('File {} not found.'.format(args.fastfitfile[0])) + return + + fastfit, metadata, coadd_type, fastphot = read_fastspecfit(args.fastfitfile[0]) + + # parse the targetids optional input + if args.targetids: + targetids = [int(x) for x in args.targetids.split(',')] + keep = np.where(np.isin(fastfit['TARGETID'], targetids))[0] + if len(keep) == 0: + log.warning('No matching targetids found!') + return + fastfit = fastfit[keep] + metadata = metadata[keep] + + if args.ntargets is not None: + keep = np.arange(args.ntargets) + args.firsttarget + log.info('Keeping {} targets.'.format(args.ntargets)) + fastfit = fastfit[keep] + metadata = metadata[keep] + + fastfit, metadata = select(fastfit, metadata, coadd_type, healpixels=args.healpix, + tiles=args.tile, nights=args.night) + + if coadd_type == 'custom' and args.redrockfiles is None: + errmsg = 'redrockfiles input is required if coadd_type==custom' + log.critical(errmsg) + raise IOError(errmsg) + + # check for the input redshifts header card + hdr = fitsio.read_header(args.fastfitfile[0]) + if 'INPUTZ' in hdr and hdr['INPUTZ']: + inputz = True + else: + inputz = False + + if 'NOSCORR' in hdr and hdr['NOSCORR']: + no_smooth_continuum = True + else: + no_smooth_continuum = False + + if args.outdir: + if not os.path.isdir(args.outdir): + os.makedirs(args.outdir, exist_ok=True) + + pngfile = get_qa_filename(metadata, coadd_type, outprefix=args.outprefix, + outdir=args.outdir, fastphot=fastphot, log=log) + + # are we overwriting? + if args.overwrite is False: + I = np.array([not os.path.isfile(_pngfile) for _pngfile in pngfile]) + J = ~I + if np.sum(J) > 0: + log.info('Skipping {} existing QA files.'.format(np.sum(J))) + fastfit = fastfit[I] + metadata = metadata[I] + + if len(metadata) == 0: + log.info('Done making all QA files!') + return + + log.info('Building QA for {} objects.'.format(len(metadata))) + + # Initialize the I/O class. + + Spec = DESISpectra(redux_dir=args.redux_dir, fphotodir=args.fphotodir, + fphotoinfo=args.fphotoinfo, mapdir=args.mapdir) + + templates = get_templates_filename(templateversion=args.templateversion, imf=args.imf) + + def _wrap_qa(redrockfile, indx=None, stackfit=False): + if indx is None: + indx = np.arange(len(fastfit)) + + if stackfit: + stackids = fastfit['STACKID'][indx] + data = Spec.read_stacked(redrockfile, stackids=stackids, mp=args.mp) + args.nophoto = True # force nophoto=True + + minspecwave = np.min(data[0]['coadd_wave']) - 20 + maxspecwave = np.max(data[0]['coadd_wave']) + 20 + else: + targetids = fastfit['TARGETID'][indx] + if inputz: + input_redshifts = fastfit['Z'][indx] + else: + input_redshifts = None + Spec.select(redrockfiles=redrockfile, targetids=targetids, + input_redshifts=input_redshifts, + redrockfile_prefix=args.redrockfile_prefix, + specfile_prefix=args.specfile_prefix, + qnfile_prefix=args.qnfile_prefix) + data = Spec.read_and_unpack(fastphot=fastphot, synthphot=True, mp=args.mp) + + minspecwave = args.minspecwave + maxspecwave = args.maxspecwave + + qaargs = [(data[igal], fastfit[indx[igal]], metadata[indx[igal]], + templates, coadd_type, Spec.fphoto, minspecwave, maxspecwave, + args.minphotwave, args.maxphotwave, args.emline_snrmin, + args.nsmoothspec, fastphot, args.nophoto, stackfit, inputz, + no_smooth_continuum, args.outdir, args.outprefix, log) + for igal in np.arange(len(indx))] + if args.mp > 1: + import multiprocessing + with multiprocessing.Pool(args.mp) as P: + P.map(_desiqa_one, qaargs) + else: + [desiqa_one(*_qaargs) for _qaargs in qaargs] + + t0 = time.time() + if coadd_type == 'healpix': + if args.redrockfiles is not None: + for redrockfile in args.redrockfiles: + _wrap_qa(redrockfile) + else: + allspecprods = metadata['SPECPROD'].data + allsurveys = metadata['SURVEY'].data + allprograms = metadata['PROGRAM'].data + allpixels = metadata['HEALPIX'].data + for specprod in set(allspecprods): + for survey in set(allsurveys): + for program in set(allprograms): + for pixel in set(allpixels): + indx = np.where((specprod == allspecprods) * (survey == allsurveys) * + (program == allprograms) * (pixel == allpixels))[0] + if len(indx) == 0: + #log.warning('No object found with specprod={}, survey={}, program={}, and healpixel={}!'.format( + # specprod, survey, program, pixel)) + continue + redrockfile = os.path.join(args.redux_dir, specprod, 'healpix', str(survey), str(program), str(pixel // 100), + str(pixel), 'redrock-{}-{}-{}.fits'.format(survey, program, pixel)) + _wrap_qa(redrockfile, indx) + elif coadd_type == 'custom': + for redrockfile in args.redrockfiles: + _wrap_qa(redrockfile) + elif coadd_type == 'stacked': + for redrockfile in args.redrockfiles: + _wrap_qa(redrockfile, stackfit=True) + else: + allspecprods = metadata['SPECPROD'].data + alltiles = metadata['TILEID'].astype(str).data + allnights = metadata['NIGHT'].astype(str).data + allpetals = metadata['FIBER'].data // 500 + if coadd_type == 'cumulative': + for specprod in set(allspecprods): + for tile in set(alltiles): + for petal in set(allpetals): + indx = np.where((specprod == allspecprods) * (tile == alltiles) * (petal == allpetals))[0] + if len(indx) == 0: + #log.warning('No object found with tileid={} and petal={}!'.format( + # tile, petal)) + continue + redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'cumulative', str(tile), allnights[indx[0]], + 'redrock-{}-{}-thru{}.fits'.format(petal, tile, allnights[indx[0]])) + _wrap_qa(redrockfile, indx) + elif coadd_type == 'pernight': + for specprod in set(allspecprods): + for night in set(allnights): + for tile in set(alltiles): + for petal in set(allpetals): + indx = np.where((specprod == allspecprods) * (night == allnights) * + (tile == alltiles) * (petal == allpetals))[0] + if len(indx) == 0: + continue + redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'pernight', str(tile), str(night), + 'redrock-{}-{}-{}.fits'.format(petal, tile, night)) + _wrap_qa(redrockfile, indx) + elif coadd_type == 'perexp': + allexpids = metadata['EXPID'].data + for specprod in set(allspecprods): + for night in set(allnights): + for expid in set(allexpids): + for tile in set(alltiles): + for petal in set(allpetals): + indx = np.where((specprod == allspecprods) * (night == allnights) * + (expid == allexpids) * (tile == alltiles) * + (petal == allpetals))[0] + if len(indx) == 0: + continue + redrockfile = os.path.join(args.redux_dir, specprod, 'tiles', 'perexp', str(tile), '{:08d}'.format(expid), + 'redrock-{}-{}-exp{:08d}.fits'.format(petal, tile, expid)) + _wrap_qa(redrockfile, indx) + + log.info('QA for everything took: {:.2f} sec'.format(time.time()-t0)) + diff --git a/py/fastspecfit/test/build_test_templates.py b/py/fastspecfit/test/build_test_templates.py deleted file mode 100644 index 6b1e7168..00000000 --- a/py/fastspecfit/test/build_test_templates.py +++ /dev/null @@ -1,49 +0,0 @@ -"""Build miniature template set. - -""" -import os, pdb -import numpy as np -import fitsio -from astropy.table import Table -from astropy.io import fits -from pkg_resources import resource_filename - -os.environ['FASTSPECFIT_TEMPLATES'] = os.path.join(os.environ.get('DESI_ROOT'), 'science', 'gqp', 'templates', 'SSP-CKC14z') - -sspversion = 'v1.0' -isochrone = 'Padova' -library = 'CKC14z' -imf = 'Kroupa' -metallicity = 'Z0.0190' - -templates_dir = os.environ.get('FASTSPECFIT_TEMPLATES', '.') -sspfile = os.path.join(templates_dir, sspversion, 'SSP_{}_{}_{}_{}.fits'.format( - isochrone, library, imf, metallicity)) -if not os.path.isfile(sspfile): - raise IOError('SSP templates file not found {}'.format(sspfile)) - -outfile = os.path.join(resource_filename('fastspecfit.test', 'data'), 'SSP_{}_{}_{}_{}.fits'.format( - isochrone, library, imf, metallicity)) - -print('Reading {}'.format(sspfile)) -wave, wavehdr = fitsio.read(sspfile, 'WAVE', header=True) -flux, fluxhdr = fitsio.read(sspfile, 'FLUX', header=True) -meta, metahdr = fitsio.read(sspfile, 'METADATA', header=True) -meta = Table(meta) - -# Trim the wavelengths and select the number/ages of the templates. -# https://www.sdss.org/dr14/spectro/galaxy_mpajhu -minwave, maxwave = 500.0, 30e4 -keep = np.where((wave >= minwave) * (wave <= maxwave))[0] -wave = wave[keep] - -myages = np.array([0.005, 0.01, 0.025, 0.05, 0.1, 0.15, 0.2, 0.4, 0.6, 0.9, 1.1, 1.4, 2.5, 5, 10.0, 13.3])*1e9 -iage = np.array([np.argmin(np.abs(meta['age']-myage)) for myage in myages]) -flux = flux[:, iage][keep, :] # flux[keep, ::5] -meta = meta[iage] - -print('Writing {}'.format(outfile)) -fitsio.write(outfile, flux, 'FLUX', header=fluxhdr, clobber=True) -fitsio.write(outfile, wave, 'WAVE', header=wavehdr) -fitsio.write(outfile, meta.as_array(), 'METADATA', header=metahdr) - diff --git a/py/fastspecfit/test/data/SSP_Padova_CKC14z_Kroupa_Z0.0190.fits b/py/fastspecfit/test/data/SSP_Padova_CKC14z_Kroupa_Z0.0190.fits deleted file mode 100644 index 4c2f9c29..00000000 Binary files a/py/fastspecfit/test/data/SSP_Padova_CKC14z_Kroupa_Z0.0190.fits and /dev/null differ diff --git a/py/fastspecfit/webapp/fastmodel/models.py b/py/fastspecfit/webapp/fastmodel/models.py index 9dbcd119..6b9a6533 100644 --- a/py/fastspecfit/webapp/fastmodel/models.py +++ b/py/fastspecfit/webapp/fastmodel/models.py @@ -154,27 +154,35 @@ class FastModel(Model): flux_synth_photmodel_w3 = FloatField(null=True) flux_synth_photmodel_w4 = FloatField(null=True) - kcorr_u = FloatField(null=True) - kcorr_b = FloatField(null=True) - kcorr_v = FloatField(null=True) - kcorr_w1 = FloatField(null=True) - kcorr_w2 = FloatField(null=True) - absmag_u = FloatField(null=True) - absmag_b = FloatField(null=True) - absmag_v = FloatField(null=True) - absmag_w1 = FloatField(null=True) - absmag_w2 = FloatField(null=True) - - kcorr_sdss_u = FloatField(null=True) - kcorr_sdss_g = FloatField(null=True) - kcorr_sdss_r = FloatField(null=True) - kcorr_sdss_i = FloatField(null=True) - kcorr_sdss_z = FloatField(null=True) - absmag_sdss_u = FloatField(null=True) - absmag_sdss_g = FloatField(null=True) - absmag_sdss_r = FloatField(null=True) - absmag_sdss_i = FloatField(null=True) - absmag_sdss_z = FloatField(null=True) + kcorr10_decam_g = FloatField(null=True) + kcorr10_decam_r = FloatField(null=True) + kcorr10_decam_z = FloatField(null=True) + absmag10_decam_g = FloatField(null=True) + absmag10_decam_r = FloatField(null=True) + absmag10_decam_z = FloatField(null=True) + + kcorr00_u = FloatField(null=True) + kcorr00_b = FloatField(null=True) + kcorr00_v = FloatField(null=True) + absmag00_u = FloatField(null=True) + absmag00_b = FloatField(null=True) + absmag00_v = FloatField(null=True) + + kcorr01_sdss_u = FloatField(null=True) + kcorr01_sdss_g = FloatField(null=True) + kcorr01_sdss_r = FloatField(null=True) + kcorr01_sdss_i = FloatField(null=True) + kcorr01_sdss_z = FloatField(null=True) + absmag01_sdss_u = FloatField(null=True) + absmag01_sdss_g = FloatField(null=True) + absmag01_sdss_r = FloatField(null=True) + absmag01_sdss_i = FloatField(null=True) + absmag01_sdss_z = FloatField(null=True) + + kcorr01_w1 = FloatField(null=True) + #kcorr_w2 = FloatField(null=True) + absmag01_w1 = FloatField(null=True) + #absmag01_w2 = FloatField(null=True) abmag_g = CharField(max_length=50, default='') abmag_r = CharField(max_length=50, default='') @@ -216,6 +224,9 @@ class FastModel(Model): loglnu_1500 = FloatField(null=True) loglnu_2800 = FloatField(null=True) + logl_1450 = FloatField(null=True) + logl_1700 = FloatField(null=True) + logl_3000 = FloatField(null=True) logl_5100 = FloatField(null=True) apercorr = FloatField(null=True) apercorr_g = FloatField(null=True) @@ -223,7 +234,8 @@ class FastModel(Model): apercorr_z = FloatField(null=True) rchi2_line = FloatField(null=True) - delta_linerchi2 = FloatField(null=True) + delta_linechi2 = FloatField(null=True) + delta_linendof = IntegerField(null=True) narrow_z = FloatField(null=True) broad_z = FloatField(null=True) diff --git a/py/fastspecfit/webapp/load.py b/py/fastspecfit/webapp/load.py index 7fb313f5..9a32bb25 100644 --- a/py/fastspecfit/webapp/load.py +++ b/py/fastspecfit/webapp/load.py @@ -144,45 +144,58 @@ def main(): 'FLUX_SYNTH_PHOTMODEL_W2', 'FLUX_SYNTH_PHOTMODEL_W3', 'FLUX_SYNTH_PHOTMODEL_W4', - 'KCORR_U', - 'ABSMAG_U', - 'ABSMAG_IVAR_U', - 'KCORR_B', - 'ABSMAG_B', - 'ABSMAG_IVAR_B', - 'KCORR_V', - 'ABSMAG_V', - 'ABSMAG_IVAR_V', - 'KCORR_SDSS_U', - 'ABSMAG_SDSS_U', - 'ABSMAG_IVAR_SDSS_U', - 'KCORR_SDSS_G', - 'ABSMAG_SDSS_G', - 'ABSMAG_IVAR_SDSS_G', - 'KCORR_SDSS_R', - 'ABSMAG_SDSS_R', - 'ABSMAG_IVAR_SDSS_R', - 'KCORR_SDSS_I', - 'ABSMAG_SDSS_I', - 'ABSMAG_IVAR_SDSS_I', - 'KCORR_SDSS_Z', - 'ABSMAG_SDSS_Z', - 'ABSMAG_IVAR_SDSS_Z', - 'KCORR_W1', - 'ABSMAG_W1', - 'ABSMAG_IVAR_W1', - 'KCORR_W2', - 'ABSMAG_W2', - 'ABSMAG_IVAR_W2', + 'ABSMAG10_DECAM_G', + 'ABSMAG10_IVAR_DECAM_G', + 'KCORR10_DECAM_G', + 'ABSMAG10_DECAM_R', + 'ABSMAG10_IVAR_DECAM_R', + 'KCORR10_DECAM_R', + 'ABSMAG10_DECAM_Z', + 'ABSMAG10_IVAR_DECAM_Z', + 'KCORR10_DECAM_Z', + 'ABSMAG00_U', + 'ABSMAG00_IVAR_U', + 'KCORR00_U', + 'ABSMAG00_B', + 'ABSMAG00_IVAR_B', + 'KCORR00_B', + 'ABSMAG00_V', + 'ABSMAG00_IVAR_V', + 'KCORR00_V', + 'ABSMAG01_SDSS_U', + 'ABSMAG01_IVAR_SDSS_U', + 'KCORR01_SDSS_U', + 'ABSMAG01_SDSS_G', + 'ABSMAG01_IVAR_SDSS_G', + 'KCORR01_SDSS_G', + 'ABSMAG01_SDSS_R', + 'ABSMAG01_IVAR_SDSS_R', + 'KCORR01_SDSS_R', + 'ABSMAG01_SDSS_I', + 'ABSMAG01_IVAR_SDSS_I', + 'KCORR01_SDSS_I', + 'ABSMAG01_SDSS_Z', + 'ABSMAG01_IVAR_SDSS_Z', + 'KCORR01_SDSS_Z', + 'ABSMAG01_W1', + 'ABSMAG01_IVAR_W1', + 'KCORR01_W1', + #'ABSMAG_W2', + #'ABSMAG_IVAR_W2', + #'KCORR_W2', 'LOGLNU_1500', 'LOGLNU_2800', + 'LOGL_1450', + 'LOGL_1700', + 'LOGL_3000', 'LOGL_5100', 'APERCORR', 'APERCORR_G', 'APERCORR_R', 'APERCORR_Z', 'RCHI2_LINE', - 'DELTA_LINERCHI2', + 'DELTA_LINECHI2', + 'DELTA_LINENDOF', 'NARROW_Z', 'BROAD_Z', 'UV_Z', diff --git a/py/fastspecfit/webapp/templates/target.html b/py/fastspecfit/webapp/templates/target.html index d5813322..f2028365 100644 --- a/py/fastspecfit/webapp/templates/target.html +++ b/py/fastspecfit/webapp/templates/target.html @@ -385,13 +385,25 @@
Physical Properties & Luminosities
Z/Zsun - AV
(mag) - SFR
(Msun/yr) + AV + SFR Age
(Gyr) log10(M/Msun) - log10(Lν, 1500 Å)
(10-28 erg s-1Hz-1) - log10(Lν, 2800 Å)
(10-28 erg s-1Hz-1) - log10(L5100 Å)
(1010 Lsun) + log10(Lν, 1500 Å) + log10(Lν, 2800 Å) + log10(L1450 Å) + log10(L1700 Å) + log10(L3000 Å) + log10(L5100 Å) + + +   + (mag) + (Msun/yr) + (Gyr) +   + (10-28 erg s-1Hz-1) + (1010 Lsun) @@ -404,6 +416,9 @@
Physical Properties & Luminosities
{{ target.logmstar|floatformat:3 }} {{ target.loglnu_1500|floatformat:3 }} {{ target.loglnu_2800|floatformat:3 }} + {{ target.logl_1450|floatformat:3 }} + {{ target.logl_1700|floatformat:3 }} + {{ target.logl_3000|floatformat:3 }} {{ target.logl_5100|floatformat:3 }} @@ -417,50 +432,71 @@
Rest-Frame Photometry & K-corrections

- SDSS + DECam + SDSS/WISE Johnson-Morgan +   + 1.0g + 1.0r + 1.0z 0.1u 0.1g 0.1r 0.1i 0.1z + 0.1W1 0.0U 0.0B 0.0V + Mabs - {{ target.absmag_sdss_u|floatformat:3 }} - {{ target.absmag_sdss_g|floatformat:3 }} - {{ target.absmag_sdss_r|floatformat:3 }} - {{ target.absmag_sdss_i|floatformat:3 }} - {{ target.absmag_sdss_z|floatformat:3 }} - {{ target.absmag_u|floatformat:3 }} - {{ target.absmag_b|floatformat:3 }} - {{ target.absmag_v|floatformat:3 }} + {{ target.absmag10_decam_g|floatformat:3 }} + {{ target.absmag10_decam_r|floatformat:3 }} + {{ target.absmag10_decam_z|floatformat:3 }} + {{ target.absmag01_sdss_u|floatformat:3 }} + {{ target.absmag01_sdss_g|floatformat:3 }} + {{ target.absmag01_sdss_r|floatformat:3 }} + {{ target.absmag01_sdss_i|floatformat:3 }} + {{ target.absmag01_sdss_z|floatformat:3 }} + {{ target.absmag01_w1|floatformat:3 }} + {{ target.absmag00_u|floatformat:3 }} + {{ target.absmag00_b|floatformat:3 }} + {{ target.absmag00_v|floatformat:3 }} + Kcorr - {{ target.kcorr_sdss_u|floatformat:3 }} - {{ target.kcorr_sdss_g|floatformat:3 }} - {{ target.kcorr_sdss_r|floatformat:3 }} - {{ target.kcorr_sdss_i|floatformat:3 }} - {{ target.kcorr_sdss_z|floatformat:3 }} - {{ target.kcorr_u|floatformat:3 }} - {{ target.kcorr_b|floatformat:3 }} - {{ target.kcorr_v|floatformat:3 }} + {{ target.kcorr10_decam_g|floatformat:3 }} + {{ target.kcorr10_decam_r|floatformat:3 }} + {{ target.kcorr10_decam_z|floatformat:3 }} + {{ target.kcorr01_sdss_u|floatformat:3 }} + {{ target.kcorr01_sdss_g|floatformat:3 }} + {{ target.kcorr01_sdss_r|floatformat:3 }} + {{ target.kcorr01_sdss_i|floatformat:3 }} + {{ target.kcorr01_sdss_z|floatformat:3 }} + {{ target.kcorr01_w1|floatformat:3 }} + {{ target.kcorr00_u|floatformat:3 }} + {{ target.kcorr00_b|floatformat:3 }} + {{ target.kcorr00_v|floatformat:3 }} + @@ -485,12 +521,13 @@
Spectral Fitting
- + - + + @@ -499,7 +536,8 @@
Spectral Fitting
- + +
   Doublet Ratios
χ2ν,line Δχ2ν,nobroadΔχ2nobroadΔνnobroad Mg II
λ2796/λ2803		 [OII]
λ3726/λ3729		 [SII]
λ6731/λ6716
{{ target.rchi2_line|floatformat:3 }}{{ target.delta_linerchi2|floatformat:3 }}{{ target.delta_linechi2|floatformat:3 }}{{ target.delta_linendof }} {{ target.mgii_doublet_ratio|floatformat:3 }} {{ target.oii_doublet_ratio|floatformat:3 }} {{ target.sii_doublet_ratio|floatformat:3 }}