"""
fastspecfit.continuum
=====================
Methods and tools for continuum-fitting.
"""
import pdb # for debugging
import os, time
import numpy as np
from astropy.table import Table
from fastspecfit.io import FLUXNORM
from fastspecfit.util import C_LIGHT
# SPS template constants (used by build-fsps-templates)
PIXKMS_BLU = 25. # [km/s]
PIXKMS_RED = 100. # [km/s]
PIXKMS_WAVESPLIT = 1e4 # [Angstrom]
[docs]
def _smooth_continuum(wave, flux, ivar, redshift, camerapix=None, medbin=175,
smooth_window=75, smooth_step=25, maskkms_uv=3000.0,
maskkms_balmer=1000.0, maskkms_narrow=200.0,
linetable=None, emlinesfile=None, linemask=None, png=None,
log=None, verbose=False):
"""Build a smooth, nonparametric continuum spectrum.
Parameters
----------
wave : :class:`numpy.ndarray` [npix]
Observed-frame wavelength array.
flux : :class:`numpy.ndarray` [npix]
Spectrum corresponding to `wave`.
ivar : :class:`numpy.ndarray` [npix]
Inverse variance spectrum corresponding to `flux`.
redshift : :class:`float`
Object redshift.
medbin : :class:`int`, optional, defaults to 150 pixels
Width of the median-smoothing kernel in pixels; a magic number.
smooth_window : :class:`int`, optional, defaults to 50 pixels
Width of the sliding window used to compute the iteratively clipped
statistics (mean, median, sigma); a magic number. Note: the nominal
extraction width (0.8 A) and observed-frame wavelength range
(3600-9800 A) corresponds to pixels that are 66-24 km/s. So
`smooth_window` of 50 corresponds to 3300-1200 km/s, which is
conservative for all but the broadest lines. A magic number.
smooth_step : :class:`int`, optional, defaults to 10 pixels
Width of the step size when computing smoothed statistics; a magic
number.
maskkms_uv : :class:`float`, optional, defaults to 3000 km/s
Masking width for UV emission lines. Pixels within +/-3*maskkms_uv
are masked before median-smoothing.
maskkms_balmer : :class:`float`, optional, defaults to 3000 km/s
Like `maskkms_uv` but for Balmer lines.
maskkms_narrow : :class:`float`, optional, defaults to 300 km/s
Like `maskkms_uv` but for narrow, forbidden lines.
linemask : :class:`numpy.ndarray` of type :class:`bool`, optional, defaults to `None`
Boolean mask with the same number of pixels as `wave` where `True`
means a pixel is (possibly) affected by an emission line
(specifically a strong line which likely cannot be median-smoothed).
png : :class:`str`, optional, defaults to `None`
Generate a simple QA plot and write it out to this filename.
Returns
-------
smooth :class:`numpy.ndarray` [npix]
Smooth continuum spectrum which can be subtracted from `flux` in
order to create a pure emission-line spectrum.
smoothsigma :class:`numpy.ndarray` [npix]
Smooth one-sigma uncertainty spectrum.
"""
from numpy.lib.stride_tricks import sliding_window_view
from scipy.ndimage import median_filter
from scipy.stats import sigmaclip
#from astropy.stats import sigma_clip
from scipy.interpolate import make_interp_spline
if log is None:
from desiutil.log import get_logger, DEBUG
if verbose:
log = get_logger(DEBUG)
else:
log = get_logger()
if linetable is None:
from fastspecfit.emlines import read_emlines
linetable = read_emlines(emlinesfile=emlinesfile)
npix = len(wave)
# If we're not given a linemask, make a conservative one.
if linemask is None:
linemask = np.zeros(npix, bool) # True = (possibly) affected by emission line
nsig = 3
# select just strong lines
zlinewaves = linetable['restwave'] * (1. + redshift)
inrange = (zlinewaves > np.min(wave)) * (zlinewaves < np.max(wave))
if np.sum(inrange) > 0:
linetable = linetable[inrange]
linetable = linetable[linetable['amp'] >= 1]
if len(linetable) > 0:
for oneline in linetable:
zlinewave = oneline['restwave'] * (1. + redshift)
if oneline['isbroad']:
if oneline['isbalmer']:
sigma = maskkms_balmer
else:
sigma = maskkms_uv
else:
sigma = maskkms_narrow
sigma *= zlinewave / C_LIGHT # [km/s --> Angstrom]
I = (wave >= (zlinewave - nsig*sigma)) * (wave <= (zlinewave + nsig*sigma))
if len(I) > 0:
linemask[I] = True
# Special: mask Ly-a (1215 A)
zlinewave = 1215. * (1. + redshift)
if (zlinewave > np.min(wave)) * (zlinewave < np.max(wave)):
sigma = maskkms_uv * zlinewave / C_LIGHT # [km/s --> Angstrom]
I = (wave >= (zlinewave - nsig*sigma)) * (wave <= (zlinewave + nsig*sigma))
if len(I) > 0:
linemask[I] = True
if len(linemask) != npix:
errmsg = 'Linemask must have the same number of pixels as the input spectrum.'
log.critical(errmsg)
raise ValueError(errmsg)
# If camerapix is given, mask the 10 (presumably noisy) pixels from the edge
# of each per-camera spectrum.
if camerapix is not None:
for campix in camerapix:
ivar[campix[0]:campix[0]+10] = 0.
ivar[campix[1]-10:campix[1]] = 0.
# Build the smooth (line-free) continuum by computing statistics in a
# sliding window, accounting for masked pixels and trying to be smart
# about broad lines. See:
# https://stackoverflow.com/questions/41851044/python-median-filter-for-1d-numpy-array
# https://numpy.org/devdocs/reference/generated/numpy.lib.stride_tricks.sliding_window_view.html
wave_win = sliding_window_view(wave, window_shape=smooth_window)
flux_win = sliding_window_view(flux, window_shape=smooth_window)
ivar_win = sliding_window_view(ivar, window_shape=smooth_window)
noline_win = sliding_window_view(np.logical_not(linemask), window_shape=smooth_window)
nminpix = 15
smooth_wave, smooth_flux, smooth_sigma, smooth_mask = [], [], [], []
for swave, sflux, sivar, noline in zip(wave_win[::smooth_step],
flux_win[::smooth_step],
ivar_win[::smooth_step],
noline_win[::smooth_step]):
# if there are fewer than XX good pixels after accounting for the
# line-mask, skip this window.
sflux = sflux[noline]
if len(sflux) < nminpix:
smooth_mask.append(True)
continue
swave = swave[noline]
sivar = sivar[noline]
cflux, _, _ = sigmaclip(sflux, low=2.0, high=2.0)
if len(cflux) < nminpix:
smooth_mask.append(True)
continue
# Toss out regions with too little good data.
if np.sum(sivar > 0) < nminpix:
smooth_mask.append(True)
continue
I = np.isin(sflux, cflux) # fragile?
sig = np.std(cflux) # simple median and sigma
mn = np.median(cflux)
# One more check for crummy spectral regions.
if mn == 0.0:
smooth_mask.append(True)
continue
smooth_wave.append(np.mean(swave[I]))
smooth_mask.append(False)
## astropy is too slow!!
#cflux = sigma_clip(sflux, sigma=2.0, cenfunc='median', stdfunc='std', masked=False, grow=1.5)
#if np.sum(np.isfinite(cflux)) < 10:
# smooth_mask.append(True)
# continue
#I = np.isfinite(cflux) # should never be fully masked!
#smooth_wave.append(np.mean(swave[I]))
#smooth_mask.append(False)
#sig = np.std(cflux[I])
#mn = np.median(cflux[I])
## inverse-variance weighted mean and sigma
#norm = np.sum(sivar[I])
#mn = np.sum(sivar[I] * cflux[I]) / norm # weighted mean
#sig = np.sqrt(np.sum(sivar[I] * (cflux[I] - mn)**2) / norm) # weighted sigma
smooth_sigma.append(sig)
smooth_flux.append(mn)
smooth_mask = np.array(smooth_mask)
smooth_wave = np.array(smooth_wave)
smooth_sigma = np.array(smooth_sigma)
smooth_flux = np.array(smooth_flux)
# For debugging.
if png:
_smooth_wave = smooth_wave.copy()
_smooth_flux = smooth_flux.copy()
# corner case for very wacky spectra
if len(smooth_flux) == 0:
smooth_flux = flux
#smooth_sigma = flux * 0 + np.std(flux)
else:
#smooth_flux = np.interp(wave, smooth_wave, smooth_flux)
#smooth_sigma = np.interp(wave, smooth_wave, smooth_sigma)
srt = np.argsort(smooth_wave)
bspl_flux = make_interp_spline(smooth_wave[srt], smooth_flux[srt], k=1)
smooth_flux = bspl_flux(wave)
# check for extrapolation
blu = np.where(wave < np.min(smooth_wave[srt]))[0]
red = np.where(wave > np.max(smooth_wave[srt]))[0]
if len(blu) > 0:
smooth_flux[:blu[-1]] = smooth_flux[blu[-1]]
if len(red) > 0:
smooth_flux[red[0]:] = smooth_flux[red[0]]
smooth = median_filter(smooth_flux, medbin, mode='nearest')
#smoothsigma = median_filter(smooth_sigma, medbin, mode='nearest')
#smoothsigma = median_filter(smooth_sigma, medbin, mode='nearest')
I = ivar > 0
if np.sum(I) > 0:
smoothsigma = np.zeros_like(wave) + np.median(1. / np.sqrt(ivar[I]))
smoothsigma[I] = 1. / np.sqrt(ivar[I])
# very important!
Z = (flux == 0.0) * (ivar == 0.0)
if np.sum(Z) > 0:
smooth[Z] = 0.0
# Optional QA.
if png:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(2, 1, figsize=(8, 10), sharex=True)
ax[0].plot(wave, flux)
ax[0].scatter(wave[linemask], flux[linemask], s=10, marker='s', color='k', zorder=2)
ax[0].plot(wave, smooth, color='red')
ax[0].plot(_smooth_wave, _smooth_flux, color='orange')
#ax[0].scatter(_smooth_wave, _smooth_flux, color='orange', marker='s', ls='-', s=20)
ax[1].plot(wave, flux - smooth)
ax[1].axhline(y=0, color='k')
#for xx in ax:
# #xx.set_xlim(3800, 4300)
# #xx.set_xlim(5200, 6050)
# #xx.set_xlim(7000, 9000)
# xx.set_xlim(8250, 8400)
#for xx in ax:
# xx.set_ylim(-0.2, 1.5)
zlinewaves = linetable['restwave'] * (1. + redshift)
linenames = linetable['name']
inrange = np.where((zlinewaves > np.min(wave)) * (zlinewaves < np.max(wave)))[0]
if len(inrange) > 0:
for linename, zlinewave in zip(linenames[inrange], zlinewaves[inrange]):
for xx in ax:
xx.axvline(x=zlinewave, color='gray')
fig.savefig(png)
return smooth, smoothsigma
[docs]
def _estimate_linesigmas(wave, flux, ivar, redshift=0.0, png=None,
refit=True, log=None, verbose=False):
"""Estimate the velocity width from potentially strong, isolated lines.
"""
import warnings
if log is None:
from desiutil.log import get_logger, DEBUG
if verbose:
log = get_logger(DEBUG)
else:
log = get_logger()
def _get_linesigma(zlinewaves, init_linesigma, label='Line', ax=None):
from scipy.optimize import curve_fit
linesigma, linesigma_snr = 0.0, 0.0
if ax:
_empty = True
inrange = (zlinewaves > np.min(wave)) * (zlinewaves < np.max(wave))
if np.sum(inrange) > 0:
stackdvel, stackflux, stackivar, contflux = [], [], [], []
for zlinewave in zlinewaves[inrange]:
I = ((wave >= (zlinewave - 5*init_linesigma * zlinewave / C_LIGHT)) *
(wave <= (zlinewave + 5*init_linesigma * zlinewave / C_LIGHT)) *
(ivar > 0))
J = np.logical_or(
(wave > (zlinewave - 8*init_linesigma * zlinewave / C_LIGHT)) *
(wave < (zlinewave - 5*init_linesigma * zlinewave / C_LIGHT)) *
(ivar > 0),
(wave < (zlinewave + 8*init_linesigma * zlinewave / C_LIGHT)) *
(wave > (zlinewave + 5*init_linesigma * zlinewave / C_LIGHT)) *
(ivar > 0))
if (np.sum(I) > 3) and np.max(flux[I]*ivar[I]) > 1:
stackdvel.append((wave[I] - zlinewave) / zlinewave * C_LIGHT)
norm = np.percentile(flux[I], 99)
if norm <= 0:
norm = 1.0
stackflux.append(flux[I] / norm)
stackivar.append(ivar[I] * norm**2)
if np.sum(J) > 3:
contflux.append(flux[J] / norm) # continuum pixels
#contflux.append(np.std(flux[J]) / norm) # error in the mean
#contflux.append(np.std(flux[J]) / np.sqrt(np.sum(J)) / norm) # error in the mean
else:
contflux.append(flux[I] / norm) # shouldn't happen...
#contflux.append(np.std(flux[I]) / norm) # shouldn't happen...
#contflux.append(np.std(flux[I]) / np.sqrt(np.sum(I)) / norm) # shouldn't happen...
if len(stackflux) > 0:
stackdvel = np.hstack(stackdvel)
stackflux = np.hstack(stackflux)
stackivar = np.hstack(stackivar)
contflux = np.hstack(contflux)
if len(stackflux) > 10: # require at least 10 pixels
#onegauss = lambda x, amp, sigma: amp * np.exp(-0.5 * x**2 / sigma**2) # no pedestal
onegauss = lambda x, amp, sigma, const: amp * np.exp(-0.5 * x**2 / sigma**2) + const
#onegauss = lambda x, amp, sigma, const, slope: amp * np.exp(-0.5 * x**2 / sigma**2) + const + slope*x
stacksigma = 1 / np.sqrt(stackivar)
try:
with warnings.catch_warnings():
warnings.simplefilter('ignore')
popt, _ = curve_fit(onegauss, xdata=stackdvel, ydata=stackflux,
sigma=stacksigma, p0=[1.0, init_linesigma, 0.0])
popt[1] = np.abs(popt[1])
if popt[0] > 0 and popt[1] > 0:
linesigma = popt[1]
robust_std = np.diff(np.percentile(contflux, [25, 75]))[0] / 1.349 # robust sigma
#robust_std = np.std(contflux)
if robust_std > 0:
linesigma_snr = popt[0] / robust_std
else:
linesigma_snr = 0.0
else:
popt = None
except RuntimeError:
popt = None
if ax:
_label = r'{} $\sigma$={:.0f} km/s S/N={:.1f}'.format(label, linesigma, linesigma_snr)
ax.scatter(stackdvel, stackflux, s=10, label=_label)
if popt is not None:
srt = np.argsort(stackdvel)
linemodel = onegauss(stackdvel[srt], *popt)
ax.plot(stackdvel[srt], linemodel, color='k')#, label='Gaussian Model')
else:
linemodel = stackflux * 0
#_min, _max = np.percentile(stackflux, [5, 95])
_max = np.max([np.max(linemodel), 1.05*np.percentile(stackflux, 99)])
ax.set_ylim(-2*np.median(contflux), _max)
if linesigma > 0:
if linesigma < np.max(stackdvel):
ax.set_xlim(-5*linesigma, +5*linesigma)
ax.set_xlabel(r'$\Delta v$ (km/s)')
ax.set_ylabel('Relative Flux')
ax.legend(loc='upper left', fontsize=8, frameon=False)
_empty = False
log.debug('{} masking sigma={:.0f} km/s and S/N={:.0f}'.format(label, linesigma, linesigma_snr))
if ax and _empty:
ax.plot([0, 0], [0, 0], label='{}-No Data'.format(label))
ax.axes.xaxis.set_visible(False)
ax.axes.yaxis.set_visible(False)
ax.legend(loc='upper left', fontsize=10)
return linesigma, linesigma_snr
linesigma_snr_min = 1.5
init_linesigma_balmer = 1000.0
init_linesigma_narrow = 200.0
init_linesigma_uv = 2000.0
if png:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 3, figsize=(12, 4))
else:
ax = [None] * 3
# [OII] doublet, [OIII] 4959,5007
zlinewaves = np.array([3728.483, 4960.295, 5008.239]) * (1. + redshift)
linesigma_narrow, linesigma_narrow_snr = _get_linesigma(zlinewaves, init_linesigma_narrow,
label='Forbidden', #label='[OII]+[OIII]',
ax=ax[0])
# refit with the new value
if refit and linesigma_narrow_snr > 0:
if (linesigma_narrow > init_linesigma_narrow) and (linesigma_narrow < 5*init_linesigma_narrow) and (linesigma_narrow_snr > linesigma_snr_min):
if ax[0] is not None:
ax[0].clear()
linesigma_narrow, linesigma_narrow_snr = _get_linesigma(
zlinewaves, linesigma_narrow, label='Forbidden', ax=ax[0])
if (linesigma_narrow < 50) or (linesigma_narrow > 5*init_linesigma_narrow) or (linesigma_narrow_snr < linesigma_snr_min):
linesigma_narrow_snr = 0.0
linesigma_narrow = init_linesigma_narrow
# Hbeta, Halpha
zlinewaves = np.array([4862.683, 6564.613]) * (1. + redshift)
linesigma_balmer, linesigma_balmer_snr = _get_linesigma(zlinewaves, init_linesigma_balmer,
label='Balmer', #label=r'H$\alpha$+H$\beta$',
ax=ax[1])
# refit with the new value
if refit and linesigma_balmer_snr > 0:
if (linesigma_balmer > init_linesigma_balmer) and (linesigma_balmer < 5*init_linesigma_balmer) and (linesigma_balmer_snr > linesigma_snr_min):
if ax[1] is not None:
ax[1].clear()
linesigma_balmer, linesigma_balmer_snr = _get_linesigma(zlinewaves, linesigma_balmer,
label='Balmer', #label=r'H$\alpha$+H$\beta$',
ax=ax[1])
# if no good fit, should we use narrow or Balmer??
if (linesigma_balmer < 50) or (linesigma_balmer > 5*init_linesigma_balmer) or (linesigma_balmer_snr < linesigma_snr_min):
linesigma_balmer_snr = 0.0
linesigma_balmer = init_linesigma_balmer
#linesigma_balmer = init_linesigma_narrow
# Lya, SiIV doublet, CIV doublet, CIII], MgII doublet
zlinewaves = np.array([1215.670, 1549.4795, 2799.942]) * (1. + redshift)
#zlinewaves = np.array([1215.670, 1398.2625, 1549.4795, 1908.734, 2799.942]) * (1. + redshift)
linesigma_uv, linesigma_uv_snr = _get_linesigma(zlinewaves, init_linesigma_uv,
label='UV/Broad', ax=ax[2])
# refit with the new value
if refit and linesigma_uv_snr > 0:
if (linesigma_uv > init_linesigma_uv) and (linesigma_uv < 5*init_linesigma_uv) and (linesigma_uv_snr > linesigma_snr_min):
if ax[2] is not None:
ax[2].clear()
linesigma_uv, linesigma_uv_snr = _get_linesigma(zlinewaves, linesigma_uv,
label='UV/Broad', ax=ax[2])
if (linesigma_uv < 300) or (linesigma_uv > 5*init_linesigma_uv) or (linesigma_uv_snr < linesigma_snr_min):
linesigma_uv_snr = 0.0
linesigma_uv = init_linesigma_uv
if png:
fig.subplots_adjust(left=0.1, bottom=0.15, wspace=0.2, right=0.95, top=0.95)
fig.savefig(png)
return (linesigma_narrow, linesigma_balmer, linesigma_uv,
linesigma_narrow_snr, linesigma_balmer_snr, linesigma_uv_snr)
[docs]
def restframe_photometry(redshift, zmodelflux, zmodelwave, maggies, ivarmaggies,
filters_in, absmag_filters, band_shift=None, snrmin=2.,
dmod=None, cosmo=None, log=None):
"""Compute K-corrections and rest-frame photometry for a single object.
Parameters
----------
redshift : :class:`float`
Galaxy or QSO redshift.
zmodelwave : `numpy.ndarray`
Observed-frame (redshifted) model wavelength array.
zmodelflux : `numpy.ndarray`
Observed-frame model spectrum.
maggies : `numpy.ndarray`
Input photometric fluxes in the `filters_in` bandpasses.
ivarmaggies : `numpy.ndarray`
Inverse variance photometry corresponding to `maggies`.
filters_in : `speclite.filters.FilterSequence`
Input filter curves.
absmag_filters : `speclite.filters.FilterSequence`
Filter curves corresponding to desired bandpasses.
band_shift : `numpy.ndarray` or `None`
Band-shift each bandpass in `absmag_filters` by this amount.
snrmin : :class:`float`, defaults to 2.
Minimum signal-to-noise ratio in the input photometry (`maggies`) in
order for that bandpass to be used to compute a K-correction.
dmod : :class:`float` or `None`
Distance modulus corresponding to `redshift`. Not needed if `cosmo` is
provided.
cosmo : `fastspecfit.util.TabulatedDESI` or `None`
Cosmological model class needed to compute the distance modulus.
log : `desiutil.log`
Logging object.
Returns
-------
kcorr : `numpy.ndarray`
K-corrections for each bandpass in `absmag_filters`.
absmag : `numpy.ndarray`
Absolute magnitudes, band-shifted according to `band_shift` (if
provided) for each bandpass in `absmag_filters`.
ivarabsmag : `numpy.ndarray`
Inverse variance corresponding to `absmag`.
synth_absmag : `numpy.ndarray`
Like `absmag`, but entirely based on synthesized photometry.
synth_maggies_in : `numpy.ndarray`
Synthesized input photometry (should closely match `maggies` if the
model fit is good).
Notes
-----
By default, the K-correction is computed by finding the observed-frame
bandpass closest in wavelength (and with a minimum signal-to-noise ratio) to
the desired band-shifted absolute magnitude bandpass. In other words, by
default we endeavor to minimize the K-correction. The inverse variance,
`ivarabsmag`, is derived from the inverse variance of the K-corrected
photometry. If no bandpass is available then `ivarabsmag` is set to zero and
`absmag` is derived from the synthesized rest-frame photometry.
"""
if log is None:
from desiutil.log import get_logger
log = get_logger()
nabs = len(absmag_filters)
kcorr = np.zeros(nabs, dtype='f4')
absmag = np.zeros(nabs, dtype='f4')
ivarabsmag = np.zeros(nabs, dtype='f4')
synth_absmag = np.zeros(nabs, dtype='f4')
synth_maggies_in = np.zeros(len(maggies))
if redshift <= 0.0:
errmsg = 'Input redshift not defined, zero, or negative!'
log.warning(errmsg)
return kcorr, absmag, ivarabsmag, synth_absmag, synth_maggies_in
if cosmo is None:
from fastspecfit.util import TabulatedDESI
cosmo = TabulatedDESI()
if dmod is None:
dmod = cosmo.distance_modulus(redshift)
modelwave = zmodelwave / (1. + redshift)
lambda_in = filters_in.effective_wavelengths.value
if band_shift is None:
band_shift = np.zeros_like(lambda_in)
# input bandpasses, observed frame; maggies and synth_maggies_in should be
# very close.
try:
synth_maggies_in = filters_in.get_ab_maggies(zmodelflux / FLUXNORM, zmodelwave)
except:
# pad in the case of an object at very high redshift (z>5.5).
log.warning('Padding model spectrum due to insufficient wavelength coverage to synthesize photometry.')
padflux, padwave = filters_in.pad_spectrum(zmodelflux / FLUXNORM, zmodelwave, axis=0, method='edge')
synth_maggies_in = filters_in.get_ab_maggies(padflux, padwave)
synth_maggies_in = np.array(synth_maggies_in.as_array().tolist()[0])
def _kcorr_and_absmag(filters_out, band_shift):
"""Little internal method to handle a single value of band_shift."""
nout = len(filters_out)
# note the factor of 1+band_shift
lambda_out = filters_out.effective_wavelengths.value / (1. + band_shift)
# Multiply by (1+z) to convert the best-fitting model to the "rest
# frame" and then divide by 1+band_shift to shift it and the wavelength
# vector to the band-shifted redshift. Also need one more factor of
# 1+band_shift in order maintain the AB mag normalization.
synth_outmaggies_rest = filters_out.get_ab_maggies(
zmodelflux * (1. + redshift) / (1. + band_shift) /
FLUXNORM, modelwave * (1. + band_shift))
synth_outmaggies_rest = np.array(synth_outmaggies_rest.as_array().tolist()[0]) / (1. + band_shift)
# Output bandpasses, observed frame; pad in the case of an object at
# very high redshift. Note that min(modelwave)=450 A is set in
# fastspec_one.
try:
synth_outmaggies_obs = filters_out.get_ab_maggies(zmodelflux / FLUXNORM, zmodelwave)
except:
log.warning('Padding model spectrum due to insufficient wavelength coverage to synthesize photometry.')
padflux, padwave = filters_out.pad_spectrum(
zmodelflux / FLUXNORM, zmodelwave, method='edge')
synth_outmaggies_obs = filters_out.get_ab_maggies(padflux, padwave)
synth_outmaggies_obs = np.array(synth_outmaggies_obs.as_array().tolist()[0])
kcorr = np.zeros(nout, dtype='f4')
absmag = np.zeros(nout, dtype='f4')
ivarabsmag = np.zeros(nout, dtype='f4')
synth_absmag = np.zeros(nout, dtype='f4')
for jj in np.arange(nout):
lambdadist = np.abs(lambda_in / (1. + redshift) - lambda_out[jj])
# K-correct from the nearest "good" bandpass (to minimizes the K-correction)
#oband = np.argmin(lambdadist)
#oband = np.argmin(lambdadist + (ivarmaggies == 0)*1e10)
oband = np.argmin(lambdadist + (maggies*np.sqrt(ivarmaggies) < snrmin)*1e10)
kcorr[jj] = + 2.5 * np.log10(synth_outmaggies_rest[jj] / synth_maggies_in[oband])
synth_absmag[jj] = -2.5 * np.log10(synth_outmaggies_rest[jj]) - dmod
# m_R = M_Q + DM(z) + K_QR(z) or
# M_Q = m_R - DM(z) - K_QR(z)
if maggies[oband] * np.sqrt(ivarmaggies[oband]) > snrmin:
absmag[jj] = -2.5 * np.log10(maggies[oband]) - dmod - kcorr[jj]
ivarabsmag[jj] = maggies[oband]**2 * ivarmaggies[oband] * (0.4 * np.log(10.))**2
else:
# if we use synthesized photometry then ivarabsmag is zero
# (which should never happen?)
absmag[jj] = synth_absmag[jj]
return kcorr, absmag, ivarabsmag, synth_absmag
for _band_shift in set(band_shift):
I = np.where(_band_shift == np.array(band_shift))[0]
# Unfortunately, absmag_filters is a FilterSequence object, which is an
# immutable list, so we have to calculate K-corrections using all the
# filters, every time, sigh.
_kcorr, _absmag, _ivarabsmag, _synth_absmag = _kcorr_and_absmag(absmag_filters, band_shift=_band_shift)
kcorr[I] = _kcorr[I]
absmag[I] = _absmag[I]
ivarabsmag[I] = _ivarabsmag[I]
synth_absmag[I] = _synth_absmag[I]
return kcorr, absmag, ivarabsmag, synth_absmag, synth_maggies_in
[docs]
def _convolve_vdisp(templateflux, vdisp, pixsize_kms):
"""Convolve by the velocity dispersion.
Parameters
----------
templateflux
vdisp
Returns
-------
Notes
-----
"""
from scipy.ndimage import gaussian_filter1d
if vdisp <= 0.0:
return templateflux
sigma = vdisp / pixsize_kms # [pixels]
smoothflux = gaussian_filter1d(templateflux, sigma=sigma, axis=0)
return smoothflux
class Filters(object):
def __init__(self, ignore_photometry=False, fphoto=None, load_filters=True):
"""Class to load filters, dust, and filter- and dust-related methods.
"""
super(Filters, self).__init__()
from speclite import filters
if fphoto is not None:
keys = fphoto.keys()
self.uniqueid = fphoto['uniqueid']
self.photounits = fphoto['photounits']
if 'readcols' in keys:
self.readcols = np.array(fphoto['readcols'])
if 'outcols' in keys:
self.outcols = np.array(fphoto['outcols'])
self.bands = np.array(fphoto['bands'])
self.bands_to_fit = np.array(fphoto['bands_to_fit'])
self.fluxcols = np.array(fphoto['fluxcols'])
self.fluxivarcols = np.array(fphoto['fluxivarcols'])
self.min_uncertainty = np.array(fphoto['min_uncertainty'])
if 'legacysurveydr' in keys:
self.legacysurveydr = fphoto['legacysurveydr']
if 'viewer_layer' in keys:
self.viewer_layer = fphoto['viewer_layer']
if 'viewer_pixscale' in keys:
self.viewer_pixscale = fphoto['viewer_pixscale']
if 'synth_bands' in keys:
self.synth_bands = np.array(fphoto['synth_bands'])
if 'fiber_bands' in keys:
self.fiber_bands = np.array(fphoto['fiber_bands'])
self.absmag_bands = np.array(fphoto['absmag_bands'])
self.band_shift = np.array(fphoto['band_shift'])
if load_filters:
# If fphoto['filters'] is a dictionary, then assume that there
# are N/S filters (as indicated by photsys).
self.filters = {}
for key in fphoto['filters'].keys():
self.filters[key] = filters.FilterSequence([filters.load_filter(filtname)
for filtname in fphoto['filters'][key]])
self.synth_filters = {}
for key in fphoto['synth_filters'].keys():
self.synth_filters[key] = filters.FilterSequence([filters.load_filter(filtname)
for filtname in fphoto['synth_filters'][key]])
if hasattr(self, 'fiber_bands'):
self.fiber_filters = {}
for key in fphoto['fiber_filters'].keys():
self.fiber_filters[key] = filters.FilterSequence([filters.load_filter(filtname)
for filtname in fphoto['fiber_filters'][key]])
# Simple list of filters.
self.absmag_filters = filters.FilterSequence([filters.load_filter(filtname) for filtname in fphoto['absmag_filters']])
else:
# This should never happen in production because we read the default
# fphoto file in io.DESISpectra.__init__.
self.uniqueid = 'TARGETID'
self.photounits = 'nanomaggies'
self.readcols = np.array(['TARGETID', 'RA', 'DEC', 'RELEASE', 'LS_ID',
'FIBERFLUX_G', 'FIBERFLUX_R', 'FIBERFLUX_Z',
'FIBERTOTFLUX_G', 'FIBERTOTFLUX_R', 'FIBERTOTFLUX_Z'])
self.outcols = np.array(['PHOTSYS', 'LS_ID'])
self.bands = np.array(['g', 'r', 'z', 'W1', 'W2', 'W3', 'W4'])
self.bands_to_fit = np.ones(len(self.bands), bool)
self.fluxcols = np.array(['FLUX_G', 'FLUX_R', 'FLUX_Z',
'FLUX_W1', 'FLUX_W2', 'FLUX_W3', 'FLUX_W4'])
self.fluxivarcols = np.array(['FLUX_IVAR_G', 'FLUX_IVAR_R', 'FLUX_IVAR_Z',
'FLUX_IVAR_W1', 'FLUX_IVAR_W2', 'FLUX_IVAR_W3', 'FLUX_IVAR_W4'])
self.min_uncertainty = np.array([0.02, 0.02, 0.02, 0.05, 0.05, 0.05, 0.05]) # mag
self.legacysurveydr = 'dr9'
self.viewer_layer = 'ls-dr9'
self.viewer_pixscale = 0.262
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 = ['decam_g', 'decam_r', 'decam_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',
'wise2010-W1', 'wise2010-W2', 'wise2010-W3', 'wise2010-W4'),
'S': filters.load_filters('decam2014-g', 'decam2014-r', 'decam2014-z',
'wise2010-W1', 'wise2010-W2', 'wise2010-W3', 'wise2010-W4')}
self.synth_filters = {'N': filters.load_filters('BASS-g', 'BASS-r', 'MzLS-z'),
'S': filters.load_filters('decam2014-g', 'decam2014-r', 'decam2014-z')}
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')))
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)
raise ValueError(errmsg)
# Do not fit the photometry.
if ignore_photometry:
self.bands_to_fit *= [False]
@staticmethod
def parse_photometry(bands, maggies, lambda_eff, ivarmaggies=None,
nanomaggies=True, nsigma=2.0, min_uncertainty=None,
log=None, verbose=False):
"""Parse input (nano)maggies to various outputs and pack into a table.
Parameters
----------
flam - 10-17 erg/s/cm2/A
fnu - 10-17 erg/s/cm2/Hz
abmag - AB mag
nanomaggies - input maggies are actually 1e-9 maggies
nsigma - magnitude limit
Returns
-------
phot - photometric table
Notes
-----
"""
from astropy.table import Column
if log is None:
from desiutil.log import get_logger, DEBUG
if verbose:
log = get_logger(DEBUG)
else:
log = get_logger()
shp = maggies.shape
if maggies.ndim == 1:
nband, ngal = shp[0], 1
else:
nband, ngal = shp[0], shp[1]
phot = Table()
phot.add_column(Column(name='band', data=bands))
phot.add_column(Column(name='lambda_eff', length=nband, dtype='f4'))
phot.add_column(Column(name='nanomaggies', length=nband, shape=(ngal, ), dtype='f4'))
phot.add_column(Column(name='nanomaggies_ivar', length=nband, shape=(ngal, ), dtype='f4'))
phot.add_column(Column(name='flam', length=nband, shape=(ngal, ), dtype='f8')) # note f8!
phot.add_column(Column(name='flam_ivar', length=nband, shape=(ngal, ), dtype='f8'))
phot.add_column(Column(name='abmag', length=nband, shape=(ngal, ), dtype='f4'))
phot.add_column(Column(name='abmag_ivar', length=nband, shape=(ngal, ), dtype='f4'))
#phot.add_column(Column(name='abmag_err', length=nband, shape=(ngal, ), dtype='f4'))
phot.add_column(Column(name='abmag_brighterr', length=nband, shape=(ngal, ), dtype='f4'))
phot.add_column(Column(name='abmag_fainterr', length=nband, shape=(ngal, ), dtype='f4'))
phot.add_column(Column(name='abmag_limit', length=nband, shape=(ngal, ), dtype='f4'))
if ivarmaggies is None:
ivarmaggies = np.zeros_like(maggies)
# Gaia-only targets can sometimes have grz=-99.
if np.any(ivarmaggies < 0) or np.any(maggies == -99.0):
errmsg = 'All ivarmaggies must be zero or positive!'
log.critical(errmsg)
raise ValueError(errmsg)
phot['lambda_eff'] = lambda_eff#.astype('f4')
if nanomaggies:
phot['nanomaggies'] = maggies#.astype('f4')
phot['nanomaggies_ivar'] = ivarmaggies#.astype('f4')
else:
phot['nanomaggies'] = (maggies * 1e9)#.astype('f4')
phot['nanomaggies_ivar'] = (ivarmaggies * 1e-18)#.astype('f4')
if nanomaggies:
nanofactor = 1e-9 # [nanomaggies-->maggies]
else:
nanofactor = 1.0
#print('Hack!!')
#if debug:
# maggies[3:5, 4:7] = 0.0
dims = maggies.shape
flatmaggies = maggies.flatten()
flativarmaggies = ivarmaggies.flatten()
abmag = np.zeros_like(flatmaggies)
abmag_limit = np.zeros_like(flatmaggies)
abmag_brighterr = np.zeros_like(flatmaggies)
abmag_fainterr = np.zeros_like(flatmaggies)
abmag_ivar = np.zeros_like(flatmaggies)
# deal with measurements
good = np.where(flatmaggies > 0)[0]
if len(good) > 0:
abmag[good] = -2.5 * np.log10(nanofactor * flatmaggies[good])
# deal with upper limits
snr = flatmaggies * np.sqrt(flativarmaggies)
upper = np.where((flativarmaggies > 0) * (snr <= nsigma))[0]
if len(upper) > 0:
abmag_limit[upper] = - 2.5 * np.log10(nanofactor * nsigma / np.sqrt(flativarmaggies[upper]))
# significant detections
good = np.where(snr > nsigma)[0]
if len(good) > 0:
errmaggies = 1 / np.sqrt(flativarmaggies[good])
abmag_brighterr[good] = errmaggies / (0.4 * np.log(10) * (flatmaggies[good]+errmaggies))#.astype('f4') # bright end (flux upper limit)
abmag_fainterr[good] = errmaggies / (0.4 * np.log(10) * (flatmaggies[good]-errmaggies))#.astype('f4') # faint end (flux lower limit)
abmag_ivar[good] = (flativarmaggies[good] * (flatmaggies[good] * 0.4 * np.log(10))**2)#.astype('f4')
phot['abmag'] = abmag.reshape(dims)
phot['abmag_limit'] = abmag_limit.reshape(dims)
phot['abmag_brighterr'] = abmag_brighterr.reshape(dims)
phot['abmag_fainterr'] = abmag_fainterr.reshape(dims)
phot['abmag_ivar'] = abmag_ivar.reshape(dims)
# Add a minimum uncertainty in quadrature **but only for flam**, which
# is used in the fitting.
if min_uncertainty is not None:
log.debug('Propagating minimum photometric uncertainties (mag): [{}]'.format(
' '.join(min_uncertainty.astype(str))))
good = np.where((maggies != 0) * (ivarmaggies > 0))[0]
if len(good) > 0:
factor = 2.5 / np.log(10.)
magerr = factor / (np.sqrt(ivarmaggies[good]) * maggies[good])
magerr2 = magerr**2 + min_uncertainty[good]**2
ivarmaggies[good] = factor**2 / (maggies[good]**2 * magerr2)
factor = nanofactor * 10**(-0.4 * 48.6) * C_LIGHT * 1e13 / lambda_eff**2 # [maggies-->erg/s/cm2/A]
if ngal > 1:
factor = factor[:, None] # broadcast for the models
phot['flam'] = (maggies * factor)
phot['flam_ivar'] = (ivarmaggies / factor**2)
return phot
@staticmethod
def get_dn4000(wave, flam, flam_ivar=None, redshift=None, rest=True, log=None):
"""Compute DN(4000) and, optionally, the inverse variance.
Parameters
----------
wave
flam
flam_ivar
redshift
rest
Returns
-------
Notes
-----
If `rest`=``False`` then `redshift` input is required.
Require full wavelength coverage over the definition of the index.
See eq. 11 in Bruzual 1983
(https://articles.adsabs.harvard.edu/pdf/1983ApJ...273..105B) but with
the "narrow" definition of Balogh et al. 1999.
"""
from fastspecfit.util import ivar2var
if log is None:
from desiutil.log import get_logger
log = get_logger()
dn4000, dn4000_ivar = 0.0, 0.0
if rest is False or redshift is not None:
restwave = wave / (1. + redshift) # [Angstrom]
flam2fnu = (1. + redshift) * restwave**2 / (C_LIGHT * 1e5) # [erg/s/cm2/A-->erg/s/cm2/Hz, rest]
else:
restwave = wave
flam2fnu = restwave**2 / (C_LIGHT * 1e5) # [erg/s/cm2/A-->erg/s/cm2/Hz, rest]
# Require a 2-Angstrom pad around the break definition.
wpad = 2.0
if np.min(restwave) > (3850-wpad) or np.max(restwave) < (4100+wpad):
log.warning('Too little wavelength coverage to compute Dn(4000)')
return dn4000, dn4000_ivar
fnu = flam * flam2fnu # [erg/s/cm2/Hz]
if flam_ivar is not None:
fnu_ivar = flam_ivar / flam2fnu**2
else:
fnu_ivar = np.ones_like(flam) # uniform weights
def _integrate(wave, flux, ivar, w1, w2):
from scipy import integrate, interpolate
# trim for speed
I = (wave > (w1-wpad)) * (wave < (w2+wpad))
J = np.logical_and(I, ivar > 0)
# Require no more than 20% of pixels are masked.
if np.sum(J) / np.sum(I) < 0.8:
log.warning('More than 20% of pixels in Dn(4000) definition are masked.')
return 0.0
wave = wave[J]
flux = flux[J]
ivar = ivar[J]
# should never have to extrapolate
f = interpolate.interp1d(wave, flux, assume_sorted=False, bounds_error=True)
f1 = f(w1)
f2 = f(w2)
i = interpolate.interp1d(wave, ivar, assume_sorted=False, bounds_error=True)
i1 = i(w1)
i2 = i(w2)
# insert the boundary wavelengths then integrate
I = np.where((wave > w1) * (wave < w2))[0]
wave = np.insert(wave[I], [0, len(I)], [w1, w2])
flux = np.insert(flux[I], [0, len(I)], [f1, f2])
ivar = np.insert(ivar[I], [0, len(I)], [i1, i2])
weight = integrate.simps(x=wave, y=ivar)
index = integrate.simps(x=wave, y=flux*ivar) / weight
index_var = 1 / weight
return index, index_var
blufactor = 3950.0 - 3850.0
redfactor = 4100.0 - 4000.0
try:
# yes, blue wavelength go with red integral bounds
numer, numer_var = _integrate(restwave, fnu, fnu_ivar, 4000, 4100)
denom, denom_var = _integrate(restwave, fnu, fnu_ivar, 3850, 3950)
except:
log.warning('Integration failed when computing DN(4000).')
return dn4000, dn4000_ivar
if denom == 0.0 or numer == 0.0:
log.warning('DN(4000) is ill-defined or could not be computed.')
return dn4000, dn4000_ivar
dn4000 = (blufactor / redfactor) * numer / denom
if flam_ivar is not None:
dn4000_ivar = (1.0 / (dn4000**2)) / (denom_var / (denom**2) + numer_var / (numer**2))
return dn4000, dn4000_ivar
class Inoue14(object):
def __init__(self, scale_tau=1.):
"""
IGM absorption from Inoue et al. (2014)
Parameters
----------
scale_tau : float
Parameter multiplied to the IGM :math:`\tau` values (exponential
in the linear absorption fraction).
I.e., :math:`f_\mathrm{igm} = e^{-\mathrm{scale\_tau} \tau}`.
Copyright (c) 2016-2022 Gabriel Brammer
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""
super(Inoue14, self).__init__()
self._load_data()
self.scale_tau = scale_tau
@staticmethod
def _pow(a, b):
"""C-like power, a**b
"""
return a**b
def _load_data(self):
from importlib import resources
LAF_file = resources.files('fastspecfit').joinpath('data/LAFcoeff.txt')
DLA_file = resources.files('fastspecfit').joinpath('data/DLAcoeff.txt')
data = np.loadtxt(LAF_file, unpack=True)
ix, lam, ALAF1, ALAF2, ALAF3 = data
self.lam = lam[:,np.newaxis]
self.ALAF1 = ALAF1[:,np.newaxis]
self.ALAF2 = ALAF2[:,np.newaxis]
self.ALAF3 = ALAF3[:,np.newaxis]
data = np.loadtxt(DLA_file, unpack=True)
ix, lam, ADLA1, ADLA2 = data
self.ADLA1 = ADLA1[:,np.newaxis]
self.ADLA2 = ADLA2[:,np.newaxis]
return True
@property
def NA(self):
"""
Number of Lyman-series lines
"""
return self.lam.shape[0]
def tLSLAF(self, zS, lobs):
"""
Lyman series, Lyman-alpha forest
"""
z1LAF = 1.2
z2LAF = 4.7
l2 = self.lam #[:, np.newaxis]
tLSLAF_value = np.zeros_like(lobs*l2).T
x0 = (lobs < l2*(1+zS))
x1 = x0 & (lobs < l2*(1+z1LAF))
x2 = x0 & ((lobs >= l2*(1+z1LAF)) & (lobs < l2*(1+z2LAF)))
x3 = x0 & (lobs >= l2*(1+z2LAF))
tLSLAF_value = np.zeros_like(lobs*l2)
tLSLAF_value[x1] += ((self.ALAF1/l2**1.2)*lobs**1.2)[x1]
tLSLAF_value[x2] += ((self.ALAF2/l2**3.7)*lobs**3.7)[x2]
tLSLAF_value[x3] += ((self.ALAF3/l2**5.5)*lobs**5.5)[x3]
return tLSLAF_value.sum(axis=0)
def tLSDLA(self, zS, lobs):
"""
Lyman Series, DLA
"""
z1DLA = 2.0
l2 = self.lam #[:, np.newaxis]
tLSDLA_value = np.zeros_like(lobs*l2)
x0 = (lobs < l2*(1+zS)) & (lobs < l2*(1.+z1DLA))
x1 = (lobs < l2*(1+zS)) & ~(lobs < l2*(1.+z1DLA))
tLSDLA_value[x0] += ((self.ADLA1/l2**2)*lobs**2)[x0]
tLSDLA_value[x1] += ((self.ADLA2/l2**3)*lobs**3)[x1]
return tLSDLA_value.sum(axis=0)
def tLCDLA(self, zS, lobs):
"""
Lyman continuum, DLA
"""
z1DLA = 2.0
lamL = 911.8
tLCDLA_value = np.zeros_like(lobs)
x0 = lobs < lamL*(1.+zS)
if zS < z1DLA:
tLCDLA_value[x0] = 0.2113 * self._pow(1.0+zS, 2) - 0.07661 * self._pow(1.0+zS, 2.3) * self._pow(lobs[x0]/lamL, (-3e-1)) - 0.1347 * self._pow(lobs[x0]/lamL, 2)
else:
x1 = lobs >= lamL*(1.+z1DLA)
tLCDLA_value[x0 & x1] = 0.04696 * self._pow(1.0+zS, 3) - 0.01779 * self._pow(1.0+zS, 3.3) * self._pow(lobs[x0 & x1]/lamL, (-3e-1)) - 0.02916 * self._pow(lobs[x0 & x1]/lamL, 3)
tLCDLA_value[x0 & ~x1] =0.6340 + 0.04696 * self._pow(1.0+zS, 3) - 0.01779 * self._pow(1.0+zS, 3.3) * self._pow(lobs[x0 & ~x1]/lamL, (-3e-1)) - 0.1347 * self._pow(lobs[x0 & ~x1]/lamL, 2) - 0.2905 * self._pow(lobs[x0 & ~x1]/lamL, (-3e-1))
return tLCDLA_value
def tLCLAF(self, zS, lobs):
"""
Lyman continuum, LAF
"""
z1LAF = 1.2
z2LAF = 4.7
lamL = 911.8
tLCLAF_value = np.zeros_like(lobs)
x0 = lobs < lamL*(1.+zS)
if zS < z1LAF:
tLCLAF_value[x0] = 0.3248 * (self._pow(lobs[x0]/lamL, 1.2) - self._pow(1.0+zS, -9e-1) * self._pow(lobs[x0]/lamL, 2.1))
elif zS < z2LAF:
x1 = lobs >= lamL*(1+z1LAF)
tLCLAF_value[x0 & x1] = 2.545e-2 * (self._pow(1.0+zS, 1.6) * self._pow(lobs[x0 & x1]/lamL, 2.1) - self._pow(lobs[x0 & x1]/lamL, 3.7))
tLCLAF_value[x0 & ~x1] = 2.545e-2 * self._pow(1.0+zS, 1.6) * self._pow(lobs[x0 & ~x1]/lamL, 2.1) + 0.3248 * self._pow(lobs[x0 & ~x1]/lamL, 1.2) - 0.2496 * self._pow(lobs[x0 & ~x1]/lamL, 2.1)
else:
x1 = lobs > lamL*(1.+z2LAF)
x2 = (lobs >= lamL*(1.+z1LAF)) & (lobs < lamL*(1.+z2LAF))
x3 = lobs < lamL*(1.+z1LAF)
tLCLAF_value[x0 & x1] = 5.221e-4 * (self._pow(1.0+zS, 3.4) * self._pow(lobs[x0 & x1]/lamL, 2.1) - self._pow(lobs[x0 & x1]/lamL, 5.5))
tLCLAF_value[x0 & x2] = 5.221e-4 * self._pow(1.0+zS, 3.4) * self._pow(lobs[x0 & x2]/lamL, 2.1) + 0.2182 * self._pow(lobs[x0 & x2]/lamL, 2.1) - 2.545e-2 * self._pow(lobs[x0 & x2]/lamL, 3.7)
tLCLAF_value[x0 & x3] = 5.221e-4 * self._pow(1.0+zS, 3.4) * self._pow(lobs[x0 & x3]/lamL, 2.1) + 0.3248 * self._pow(lobs[x0 & x3]/lamL, 1.2) - 3.140e-2 * self._pow(lobs[x0 & x3]/lamL, 2.1)
return tLCLAF_value
def full_IGM(self, z, lobs):
"""Get full Inoue IGM absorption
Parameters
----------
z : float
Redshift to evaluate IGM absorption
lobs : array
Observed-frame wavelength(s) in Angstroms.
Returns
-------
abs : array
IGM absorption
"""
tau_LS = self.tLSLAF(z, lobs) + self.tLSDLA(z, lobs)
tau_LC = self.tLCLAF(z, lobs) + self.tLCDLA(z, lobs)
### Upturn at short wavelengths, low-z
#k = 1./100
#l0 = 600-6/k
#clip = lobs/(1+z) < 600.
#tau_clip = 100*(1-1./(1+np.exp(-k*(lobs/(1+z)-l0))))
tau_clip = 0.
return np.exp(-self.scale_tau*(tau_LC + tau_LS + tau_clip))
def build_grid(self, zgrid, lrest):
"""Build a spline interpolation object for fast IGM models
Returns: self.interpolate
"""
from scipy.interpolate import CubicSpline
igm_grid = np.zeros((len(zgrid), len(lrest)))
for iz in range(len(zgrid)):
igm_grid[iz,:] = self.full_IGM(zgrid[iz], lrest*(1+zgrid[iz]))
self.interpolate = CubicSpline(zgrid, igm_grid)
[docs]
def continuum_specfit(data, result, templatecache, fphoto=None, emlinesfile=None,
constrain_age=False, no_smooth_continuum=False, ignore_photometry=False,
fastphot=False, log=None, debug_plots=False, verbose=False):
"""Fit the non-negative stellar continuum of a single spectrum.
Parameters
----------
data : :class:`dict`
Dictionary of input spectroscopy (plus ancillary data) populated by
:func:`fastspecfit.io.DESISpectra.read_and_unpack`.
Returns
-------
:class:`astropy.table.Table`
Table with all the continuum-fitting results with columns documented
in :func:`init_output`.
Notes
-----
- Consider using cross-correlation to update the redrock redshift.
- We solve for velocity dispersion if ...
"""
if log is None:
from desiutil.log import get_logger, DEBUG
if verbose:
log = get_logger(DEBUG)
else:
log = get_logger()
tall = time.time()
CTools = ContinuumTools(fphoto=fphoto, emlinesfile=emlinesfile,
ignore_photometry=ignore_photometry)
redshift = result['Z']
objflam = data['phot']['flam'].data * FLUXNORM
objflamivar = (data['phot']['flam_ivar'].data / FLUXNORM**2) * CTools.bands_to_fit
assert(np.all(objflamivar >= 0))
if np.any(CTools.bands_to_fit):
# Require at least one photometric optical band; do not just fit the IR
# because we will not be able to compute the aperture correction, below.
opt = (data['phot']['lambda_eff'] > 3e3) * (data['phot']['lambda_eff'] < 1e4)
if np.all(objflamivar[opt] == 0.0):
log.warning('All optical bands are masked; masking all photometry.')
objflamivar *= 0.0
# Optionally ignore templates which are older than the age of the
# universe at the redshift of the object.
def _younger_than_universe(age, tuniv, agepad=0.5):
"""Return the indices of the templates younger than the age of the universe
(plus an agepadding amount) at the given redshift. age in yr, agepad and
tuniv in Gyr
"""
return np.where(age <= 1e9 * (agepad + tuniv))[0]
nsed = len(templatecache['templateinfo'])
if constrain_age:
agekeep = _younger_than_universe(templatecache['templateinfo']['age'], data['tuniv'])
else:
agekeep = np.arange(nsed)
nage = len(agekeep)
ztemplatewave = templatecache['templatewave'] * (1. + redshift)
# Photometry-only fitting.
vdisp_nominal = templatecache['vdisp_nominal']
if fastphot:
vdispbest, vdispivar = vdisp_nominal, 0.0
log.info('Adopting nominal vdisp={:.0f} km/s.'.format(vdisp_nominal))
if np.all(objflamivar == 0):
log.info('All photometry is masked.')
coeff = np.zeros(nage, 'f4') # nage not nsed
rchi2_cont, rchi2_phot = 0.0, 0.0
dn4000_model = 0.0
sedmodel = np.zeros(len(templatecache['templatewave']))
else:
# Get the coefficients and chi2 at the nominal velocity dispersion.
t0 = time.time()
sedtemplates, sedphot = CTools.templates2data(
templatecache['templateflux_nomvdisp'][:, agekeep],
templatecache['templatewave'],
redshift=redshift, dluminosity=data['dluminosity'],
vdisp=None, synthphot=True, photsys=data['photsys'])
sedflam = sedphot['flam'].data * CTools.massnorm * FLUXNORM
coeff, rchi2_phot = CTools.call_nnls(sedflam, objflam, objflamivar)
rchi2_phot /= np.sum(objflamivar > 0) # dof???
log.info('Fitting {} models took {:.2f} seconds.'.format(
nage, time.time()-t0))
if np.all(coeff == 0):
log.warning('Continuum coefficients are all zero.')
sedmodel = np.zeros(len(templatecache['templatewave']))
dn4000_model = 0.0
else:
sedmodel = sedtemplates.dot(coeff)
# Measure Dn(4000) from the line-free model.
sedtemplates_nolines, _ = CTools.templates2data(
templatecache['templateflux_nolines_nomvdisp'][:, agekeep],
templatecache['templatewave'],
redshift=redshift, dluminosity=data['dluminosity'],
vdisp=None, synthphot=False)
sedmodel_nolines = sedtemplates_nolines.dot(coeff)
dn4000_model, _ = CTools.get_dn4000(templatecache['templatewave'],
sedmodel_nolines, rest=True, log=log)
log.info('Model Dn(4000)={:.3f}.'.format(dn4000_model))
else:
# Combine all three cameras; we will unpack them to build the
# best-fitting model (per-camera) below.
specwave = np.hstack(data['wave'])
specflux = np.hstack(data['flux'])
flamivar = np.hstack(data['ivar'])
specivar = flamivar * np.logical_not(np.hstack(data['linemask'])) # mask emission lines
if np.all(specivar == 0) or np.any(specivar < 0):
specivar = flamivar # not great...
if np.all(specivar == 0) or np.any(specivar < 0):
errmsg = 'All pixels are masked or some inverse variances are negative!'
log.critical(errmsg)
raise ValueError(errmsg)
npix = len(specwave)
# We'll need the filters for the aperture correction, below.
filters_in = CTools.synth_filters[data['photsys']]
# Prepare the spectral and photometric models at the galaxy
# redshift. And if the wavelength coverage is sufficient, also solve for
# the velocity dispersion.
#compute_vdisp = ((result['SNR_B'] > 1) and (result['SNR_R'] > 1)) and (redshift < 1.0)
restwave = specwave / (1+redshift)
Ivdisp = np.where((specivar > 0) * (restwave > 3500.0) * (restwave < 5500.0))[0]
#Ivdisp = np.where((specivar > 0) * (specsmooth != 0.0) * (restwave > 3500.0) * (restwave < 5500.0))[0]
compute_vdisp = (len(Ivdisp) > 0) and (np.ptp(restwave[Ivdisp]) > 500.0)
# stacked spectra do not have all three cameras
if 'SNR_B' in result.columns and 'SNR_R' in result.columns and 'SNR_Z' in result.columns:
log.info('S/N_b={:.2f}, S/N_r={:.2f}, S/N_z={:.2f}, rest wavelength coverage={:.0f}-{:.0f} A.'.format(
result['SNR_B'], result['SNR_R'], result['SNR_Z'], restwave[0], restwave[-1]))
if compute_vdisp:
t0 = time.time()
ztemplateflux_vdisp, _ = CTools.templates2data(
templatecache['vdispflux'], templatecache['vdispwave'], # [npix,vdispnsed,nvdisp]
redshift=redshift, dluminosity=data['dluminosity'],
specwave=data['wave'], specres=data['res'],
cameras=data['cameras'], synthphot=False, stack_cameras=True)
vdispchi2min, vdispbest, vdispivar, _ = CTools.call_nnls(
ztemplateflux_vdisp[Ivdisp, :, :],
specflux[Ivdisp], specivar[Ivdisp],
xparam=templatecache['vdisp'], xlabel=r'$\sigma$ (km/s)', log=log,
debug=debug_plots, png='deltachi2-vdisp.png')
log.info('Fitting for the velocity dispersion took {:.2f} seconds.'.format(time.time()-t0))
if vdispivar > 0:
# Require vdisp to be measured with S/N>1, which protects
# against tiny ivar becomming infinite in the output table.
vdispsnr = vdispbest * np.sqrt(vdispivar)
if vdispsnr < 1:
log.warning('vdisp signal-to-noise {:.2f} is less than one; adopting vdisp={:.0f} km/s.'.format(
vdispsnr, vdisp_nominal))
vdispbest, vdispivar = vdisp_nominal, 0.0
else:
log.info('Best-fitting vdisp={:.1f}+/-{:.1f} km/s.'.format(
vdispbest, 1/np.sqrt(vdispivar)))
else:
vdispbest = vdisp_nominal
log.info('Finding vdisp failed; adopting vdisp={:.0f} km/s.'.format(vdisp_nominal))
else:
vdispbest, vdispivar = vdisp_nominal, 0.0
log.info('Insufficient wavelength covereage to compute vdisp; adopting nominal vdisp={:.2f} km/s'.format(vdispbest))
# Derive the aperture correction.
t0 = time.time()
# First, do a quick fit of the DESI spectrum (including
# line-emission templates) so we can synthesize photometry from a
# noiseless model.
if vdispbest == vdisp_nominal:
# Use the cached templates.
use_vdisp = None
input_templateflux = templatecache['templateflux_nomvdisp'][:, agekeep]
input_templateflux_nolines = templatecache['templateflux_nolines_nomvdisp'][:, agekeep]
else:
use_vdisp = vdispbest
input_templateflux = templatecache['templateflux'][:, agekeep]
input_templateflux_nolines = templatecache['templateflux_nolines'][:, agekeep]
desitemplates, desitemplatephot = CTools.templates2data(
input_templateflux, templatecache['templatewave'],
redshift=redshift, dluminosity=data['dluminosity'],
specwave=data['wave'], specres=data['res'], specmask=data['mask'],
vdisp=use_vdisp, cameras=data['cameras'], stack_cameras=True,
synthphot=True, photsys=data['photsys'])
desitemplateflam = desitemplatephot['flam'].data * CTools.massnorm * FLUXNORM
apercorrs, apercorr = np.zeros(len(CTools.synth_bands), 'f4'), 0.0
sedtemplates, _ = CTools.templates2data(input_templateflux, templatecache['templatewave'],
vdisp=use_vdisp, redshift=redshift,
dluminosity=data['dluminosity'], synthphot=False)
if not np.any(CTools.bands_to_fit):
log.info('Skipping aperture correction since no bands were fit.')
apercorrs, apercorr = np.ones(len(CTools.synth_bands), 'f4'), 1.0
else:
# Fit using the templates with line-emission.
quickcoeff, _ = CTools.call_nnls(desitemplates, specflux, specivar)
if np.all(quickcoeff == 0):
log.warning('Quick continuum coefficients are all zero.')
else:
# Synthesize grz photometry from the full-wavelength SED to make
# sure we get the z-band correct.
quicksedflux = sedtemplates.dot(quickcoeff)
try:
quickmaggies = np.array(filters_in.get_ab_maggies(quicksedflux / FLUXNORM, ztemplatewave).as_array().tolist()[0])
except:
# pad in the case of an object at very high redshift (z>5.5).
log.warning('Padding model spectrum due to insufficient wavelength coverage to synthesize photometry.')
padflux, padwave = filters_in.pad_spectrum(quicksedflux / FLUXNORM, ztemplatewave, axis=0, method='edge')
quickmaggies = np.array(filters_in.get_ab_maggies(padflux, padwave, axis=0).as_array().tolist()[0])
quickphot = CTools.parse_photometry(CTools.synth_bands, quickmaggies,
filters_in.effective_wavelengths.value,
nanomaggies=False, log=log)
numer = np.hstack([data['phot']['nanomaggies'][data['phot']['band'] == band]
for band in CTools.synth_bands])
denom = quickphot['nanomaggies'].data
I = (numer > 0.0) * (denom > 0.0)
if np.any(I):
apercorrs[I] = numer[I] / denom[I]
I = apercorrs > 0
if np.any(I):
apercorr = np.median(apercorrs[I])
log.info('Median aperture correction = {:.3f} [{:.3f}-{:.3f}].'.format(
apercorr, np.min(apercorrs), np.max(apercorrs)))
if apercorr <= 0:
log.warning('Aperture correction not well-defined; adopting 1.0.')
apercorr = 1.0
#apercorr_g, apercorr_r, apercorr_z = apercorrs
data['apercorr'] = apercorr # needed for the line-fitting
# Perform the final fit using the line-free templates in the spectrum
# (since we mask those pixels) but the photometry synthesized from the
# templates with lines.
desitemplates_nolines, _ = CTools.templates2data(
input_templateflux_nolines, templatecache['templatewave'], redshift=redshift,
dluminosity=data['dluminosity'],
specwave=data['wave'], specres=data['res'], specmask=data['mask'],
vdisp=use_vdisp, cameras=data['cameras'], stack_cameras=True,
synthphot=False)
coeff, rchi2_cont = CTools.call_nnls(np.vstack((desitemplateflam, desitemplates_nolines)),
np.hstack((objflam, specflux * apercorr)),
np.hstack((objflamivar, specivar / apercorr**2)))
# full-continuum fitting rchi2
rchi2_cont /= (np.sum(objflamivar > 0) + np.sum(specivar > 0)) # dof???
log.info('Final fitting with {} models took {:.2f} seconds.'.format(
nage, time.time()-t0))
# rchi2 fitting just to the photometry, for analysis purposes
rchi2_phot = np.sum(objflamivar * (objflam - desitemplateflam.dot(coeff))**2)
if np.any(objflamivar > 0):
rchi2_phot /= np.sum(objflamivar > 0)
# Compute the full-wavelength best-fitting model.
if np.all(coeff == 0):
log.warning('Continuum coefficients are all zero.')
sedmodel = np.zeros(len(templatecache['templatewave']), 'f4')
desimodel = np.zeros_like(specflux)
desimodel_nolines = np.zeros_like(specflux)
dn4000_model = 0.0
else:
sedmodel = sedtemplates.dot(coeff)
desimodel = desitemplates.dot(coeff)
desimodel_nolines = desitemplates_nolines.dot(coeff)
# Measure Dn(4000) from the line-free model.
sedtemplates_nolines, _ = CTools.templates2data(
input_templateflux_nolines, templatecache['templatewave'],
vdisp=use_vdisp, redshift=redshift, dluminosity=data['dluminosity'],
synthphot=False)
sedmodel_nolines = sedtemplates_nolines.dot(coeff)
dn4000_model, _ = CTools.get_dn4000(templatecache['templatewave'], sedmodel_nolines, rest=True, log=log)
# Get DN(4000). Specivar is line-masked so we can't use it!
dn4000, dn4000_ivar = CTools.get_dn4000(specwave, specflux, flam_ivar=flamivar,
redshift=redshift, rest=False, log=log)
if dn4000_ivar > 0:
log.info('Spectroscopic DN(4000)={:.3f}+/-{:.3f}, Model Dn(4000)={:.3f}'.format(
dn4000, 1/np.sqrt(dn4000_ivar), dn4000_model))
else:
log.info('Spectroscopic DN(4000)={:.3f}, Model Dn(4000)={:.3f}'.format(
dn4000, dn4000_model))
png = None
#png = '/global/cfs/cdirs/desi/users/ioannis/tmp/junk.png'
linemask = np.hstack(data['linemask'])
if np.all(coeff == 0):
log.warning('Continuum coefficients are all zero.')
_smoothcontinuum = np.zeros_like(specwave)
else:
# Need to be careful we don't pass a large negative residual
# where there are gaps in the data.
residuals = apercorr*specflux - desimodel_nolines
I = (specflux == 0.0) * (specivar == 0.0)
if np.any(I):
residuals[I] = 0.0
_smoothcontinuum, _ = CTools.smooth_continuum(
specwave, residuals, specivar / apercorr**2,
redshift, camerapix=data['camerapix'], emlinesfile=emlinesfile,
linemask=linemask, log=log, png=png)
if no_smooth_continuum:
log.info('Zeroing out the smooth continuum correction.')
_smoothcontinuum *= 0
# Unpack the continuum into individual cameras.
continuummodel, smoothcontinuum = [], []
for camerapix in data['camerapix']:
continuummodel.append(desimodel_nolines[camerapix[0]:camerapix[1]])
smoothcontinuum.append(_smoothcontinuum[camerapix[0]:camerapix[1]])
for icam, cam in enumerate(data['cameras']):
nonzero = continuummodel[icam] != 0
if np.sum(nonzero) > 0:
corr = np.median(smoothcontinuum[icam][nonzero] / continuummodel[icam][nonzero])
result['SMOOTHCORR_{}'.format(cam.upper())] = corr * 100 # [%]
if 'SMOOTHCORR_B' in result.columns and 'SMOOTHCORR_R' in result.columns and 'SMOOTHCORR_Z' in result.columns:
log.info('Smooth continuum correction: b={:.3f}%, r={:.3f}%, z={:.3f}%'.format(
result['SMOOTHCORR_B'], result['SMOOTHCORR_R'], result['SMOOTHCORR_Z']))
# Compute K-corrections, rest-frame quantities, and physical properties.
if np.all(coeff == 0):
kcorr = np.zeros(len(CTools.absmag_bands))
absmag = np.zeros(len(CTools.absmag_bands))#-99.0
ivarabsmag = np.zeros(len(CTools.absmag_bands))
synth_bestmaggies = np.zeros(len(CTools.bands))
lums, cfluxes = {}, {}
AV, age, zzsun, logmstar, sfr = 0.0, 0.0, 0.0, 0.0, 0.0
else:
kcorr, absmag, ivarabsmag, _, synth_bestmaggies = CTools.kcorr_and_absmag(
data, templatecache['templatewave'], sedmodel, log=log)
lums, cfluxes = CTools.continuum_fluxes(data, templatecache['templatewave'], sedmodel, log=log)
# get the SPS properties
mstars = templatecache['templateinfo']['mstar'][agekeep] # [current mass in stars, Msun]
masstot = coeff.dot(mstars)
coefftot = np.sum(coeff)
logmstar = np.log10(CTools.massnorm * masstot)
zzsun = np.log10(coeff.dot(mstars * 10.**templatecache['templateinfo']['zzsun'][agekeep]) / masstot) # mass-weighted
AV = coeff.dot(templatecache['templateinfo']['av'][agekeep]) / coefftot # luminosity-weighted [mag]
age = coeff.dot(templatecache['templateinfo']['age'][agekeep]) / coefftot / 1e9 # luminosity-weighted [Gyr]
#age = coeff.dot(mstars * templatecache['templateinfo']['age'][agekeep]) / masstot / 1e9 # mass-weighted [Gyr]
sfr = CTools.massnorm * coeff.dot(templatecache['templateinfo']['sfr'][agekeep]) # [Msun/yr]
rindx = np.argmin(np.abs(CTools.absmag_filters.effective_wavelengths.value / (1.+CTools.band_shift) - 5600))
log.info(f'log(M/Msun)={logmstar:.2f}, M{CTools.absmag_bands[rindx]}={absmag[rindx]:.2f} mag, A(V)={AV:.3f}, Age={age:.3f} Gyr, SFR={sfr:.3f} Msun/yr, Z/Zsun={zzsun:.3f}')
#import matplotlib.pyplot as plt
#W = np.where(coeff > 0)[0]
#info = templatecache['templateinfo']
#sedwave = templatecache['templatewave'] * (1. + redshift)
#I = np.where((sedwave > 3000) * (sedwave < 40000))[0]
#plt.plot(sedwave[I], sedmodel[I])
#for ww in W:
# plt.plot(sedwave[I], coeff[ww] * sedtemplates[I, ww],
# label=f'{info["age"][ww]/1e9:.2f} Gyr, A(V)={info["av"][ww]:.2f} mag')
#plt.legend(fontsize=10)
#plt.scatter(sedphot['lambda_eff'][:4], objflam[:4], color='k', s=100, zorder=10)
#plt.scatter(sedphot['lambda_eff'][:4], sedflam.dot(coeff)[:4], color='gray', s=100, marker='x', zorder=11)
#plt.xscale('log')
#plt.yscale('log')
#plt.savefig('/global/cfs/cdirs/desi/users/ioannis/tmp/foo.png')
#
#print(info[W])
#print(coeff[W])
# Pack it in and return.
result['COEFF'][agekeep] = coeff
result['RCHI2_PHOT'] = rchi2_phot
result['VDISP'] = vdispbest # * u.kilometer/u.second
result['AV'] = AV # * u.mag
result['AGE'] = age # * u.Gyr
result['ZZSUN'] = zzsun
result['LOGMSTAR'] = logmstar
result['SFR'] = sfr
result['DN4000_MODEL'] = dn4000_model
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):
for lum in lums.keys():
result[lum] = lums[lum]
if bool(cfluxes):
for cflux in cfluxes.keys():
result[cflux] = cfluxes[cflux]
if not fastphot:
result['RCHI2_CONT'] = rchi2_cont
result['VDISP_IVAR'] = vdispivar # * (u.second/u.kilometer)**2
result['APERCORR'] = apercorr
for iband, band in enumerate(CTools.synth_bands):
result['APERCORR_{}'.format(band.upper())] = apercorrs[iband]
result['DN4000_OBS'] = dn4000
result['DN4000_IVAR'] = dn4000_ivar
log.info('Continuum-fitting took {:.2f} seconds.'.format(time.time()-tall))
if fastphot:
return sedmodel, None
else:
# divide out the aperture correction
continuummodel = [_continuummodel / apercorr for _continuummodel in continuummodel]
smoothcontinuum = [_smoothcontinuum / apercorr for _smoothcontinuum in smoothcontinuum]
return continuummodel, smoothcontinuum