"""
desisim.lya_spectra
===================
Function to simulate a QSO spectrum including Lyman-alpha absorption.
"""
from __future__ import division, print_function
import numpy as np
from desisim.dla import insert_dlas
from desiutil.log import get_logger
lambda_RF_LYA = 1215.67
absorber_IGM = {
'MgI(2853)' : { 'LRF':2852.96, 'COEF':1.e-4 },
'MgII(2804)' : { 'LRF':2803.5324, 'COEF':5.e-4 },
'MgII(2796)' : { 'LRF':2796.3511, 'COEF':9.e-4 },
'FeII(2600)' : { 'LRF':2600.1724835, 'COEF':1.e-4 },
'FeII(2587)' : { 'LRF':2586.6495659, 'COEF':1.e-4 },
'MnII(2577)' : { 'LRF':2576.877, 'COEF':1.e-4 },
'FeII(2383)' : { 'LRF':2382.7641781, 'COEF':1.e-4 },
'FeII(2374)' : { 'LRF':2374.4603294, 'COEF':1.e-4 },
'FeII(2344)' : { 'LRF':2344.2129601, 'COEF':1.e-4 },
'AlIII(1863)' : { 'LRF':1862.79113, 'COEF':1.e-4 },
'AlIII(1855)' : { 'LRF':1854.71829, 'COEF':1.e-4 },
'AlII(1671)' : { 'LRF':1670.7886, 'COEF':1.e-4 },
'FeII(1608)' : { 'LRF':1608.4511, 'COEF':1.e-4 },
'CIV(1551)' : { 'LRF':1550.77845, 'COEF':5.435e-4 },
'CIV(1548)' : { 'LRF':1548.2049, 'COEF':1.487e-3 },
'SiII(1527)' : { 'LRF':1526.70698, 'COEF':1.e-4 },
'SiIV(1403)' : { 'LRF':1402.77291, 'COEF':5.e-4 },
'SiIV(1394)' : { 'LRF':1393.76018, 'COEF':9.e-4 },
'CII(1335)' : { 'LRF':1334.5323, 'COEF':1.e-4 },
'SiII(1304)' : { 'LRF':1304.3702, 'COEF':1.e-4 },
'OI(1302)' : { 'LRF':1302.1685, 'COEF':1.e-4 },
'SiII(1260)' : { 'LRF':1260.4221, 'COEF':3.542e-4 },
'NV(1243)' : { 'LRF':1242.804, 'COEF':5.e-4 },
'NV(1239)' : { 'LRF':1238.821, 'COEF':5.e-4 },
'SiIII(1207)' : { 'LRF':1206.500, 'COEF':1.8919e-3 },
'NI(1200)' : { 'LRF':1200., 'COEF':1.e-3 },
'SiII(1193)' : { 'LRF':1193.2897, 'COEF':9.0776e-4 },
'SiII(1190)' : { 'LRF':1190.4158, 'COEF':6.4239e-4 },
'OI(1039)' : { 'LRF':1039.230, 'COEF':1.e-3 },
'OVI(1038)' : { 'LRF':1037.613, 'COEF':3.382-3 },
'OVI(1032)' : { 'LRF':1031.912, 'COEF':5.358e-3 },
'LYB' : { 'LRF':1025.72, 'COEF':0.1901 },
'CIII(977)' : { 'LRF':977.020, 'COEF':5.e-3 },
'OI(989)' : { 'LRF':988.7, 'COEF':1.e-3 },
'SiII(990)' : { 'LRF':989.8731, 'COEF':1.e-3 },
'LY3' : { 'LRF':972.537, 'COEF':0.0697 },
'LY4' : { 'LRF':949.7431, 'COEF':0.0335 },
'LY5' : { 'LRF':937.8035, 'COEF':0.0187 },
}
[docs]def read_lya_skewers(lyafile,indices=None,read_dlas=False,add_metals=False,add_lyb=False) :
'''
Reads Lyman alpha transmission skewers (from CoLoRe, format v2.x.y)
Args:
lyafile: full path to input FITS filename
Options:
indices: indices of input file to sub-select
read_dlas: try read DLA HDU from file
add_metals: try to read metals HDU and multiply transmission
Returns:
wave[nwave]
transmission[nlya, nwave]
metadata[nlya]
dlas[ndla] (if read_dlas=True, otherwise None)
Input file must have WAVELENGTH, TRANSMISSION, and METADATA HDUs
'''
log = get_logger()
import fitsio
h = fitsio.FITS(lyafile)
if "WAVELENGTH" in h :
wave = h["WAVELENGTH"].read()
else :
log.warning("I assume WAVELENGTH is HDU 2")
wave = h[2].read()
if "F_LYA" in h :
trans = h["F_LYA"].read()
elif "TRANSMISSION" in h:
trans=h["TRANSMISSION"].read()
else :
log.warning("I assume TRANSMISSION is HDU 3")
trans = h[3].read()
if trans.shape[1] != wave.size :
if trans.shape[0] == wave.size :
trans = trans.T # now shape is (nqso,nwave)
else :
log.error("shape of wavelength={} and transmission={} don't match".format(wave.shape,trans.shape))
raise ValueError("shape of wavelength={} and transmission={} don't match".format(wave.shape,trans.shape))
if "METADATA" in h :
meta = h["METADATA"].read()
else :
log.warning("I assume METADATA is HDU 1")
meta = h[1].read()
if (add_lyb):
if ("F_LYB" in h) :
lyb = h["F_LYB"].read()
trans*=lyb
log.info("Added LYB from transmission file")
else:
nolyb="No HDU with EXTNAME='F_LYB' in transmission file {}".format(lyafile)
log.error(nolyb)
raise KeyError(nolyb)
if add_metals:
if add_metals=='all':
#For format london v>7.3
if "F_METALS" in h:
metals = h["F_METALS"].read()
trans *= metals
log.info("Added F_Metals from transmision file")
#For format london v<7.3
elif "METALS" in h :
metals = h["METALS"].read()
trans *= metals
log.info('Added Metals from file')
else:
nom="No HDU with EXTNAME='METALS' or EXTNAME='F_METALS' in transmission file {}".format(lyafile)
log.error(nom)
raise KeyError(nom)
else:
if add_metals=='all-dev':
metal_list=['F_SI1260','F_SI1207','F_SI1193','F_SI1190']
else:
metal_list=['F_'+m for m in add_metals.split(',')]
log.info("add {} metals from transmision file".format(metal_list))
for metal in metal_list:
if (metal in h):
metals = h[metal].read()
trans *= metals
else:
nom="No HDU with EXTNAME={} in transmission file {} ".format(metal,lyafile)
log.error(nom)
raise KeyError(nom)
if indices is not None :
trans = trans[indices]
meta=meta[:][indices]
if (read_dlas):
if "DLA" in h:
dlas=h["DLA"].read()
log.info("Read DLAs from transmision file")
else:
mess="No HDU with EXTNAME='DLA' in transmission file {}".format(lyafile)
log.error(mess)
raise KeyError(mess)
else:
dlas=None
return wave,trans,meta,dlas
[docs]def apply_lya_transmission(qso_wave,qso_flux,trans_wave,trans) :
'''
Apply transmission to input flux, interpolating if needed. Note that the
transmission might include Lyman-beta and metal absorption, so we should
probably change the name of this function.
Args:
qso_wave: 1D[nwave] array of QSO wavelengths
qso_flux: 2D[nqso, nwave] array of fluxes
trans_wave: 1D[ntranswave ] array of transmission wavelength samples
trans: 2D[nqso, ntranswave] transmissions [0-1]
Returns:
output_flux[nqso, nwave]
This routine simply apply the transmission
the only thing besides multiplication is a wavelength interpolation
of transmission to the QSO wavelength grid
'''
if qso_flux.shape[0] != trans.shape[0] :
raise(ValueError("not same number of qso {} {}".format(qso_flux.shape[0],trans.shape[0])))
output_flux = qso_flux.copy()
if qso_wave.ndim == 2: # desisim.QSO(resample=True) returns a 2D wavelength array
for q in range(qso_flux.shape[0]) :
output_flux[q, :] *= np.interp(qso_wave[q, :],trans_wave,trans[q, :],left=0,right=1)
else:
for q in range(qso_flux.shape[0]) :
output_flux[q, :] *= np.interp(qso_wave,trans_wave,trans[q, :],left=0,right=1)
return output_flux
[docs]def get_spectra(lyafile, nqso=None, wave=None, templateid=None, normfilter='sdss2010-g',
seed=None, rand=None, qso=None, add_dlas=False, debug=False, nocolorcuts=True):
"""Generate a QSO spectrum which includes Lyman-alpha absorption.
Args:
lyafile (str): name of the Lyman-alpha spectrum file to read.
nqso (int, optional): number of spectra to generate (starting from the
first spectrum; if more flexibility is needed use TEMPLATEID).
wave (numpy.ndarray, optional): desired output wavelength vector.
templateid (int numpy.ndarray, optional): indices of the spectra
(0-indexed) to read from LYAFILE (default is to read everything). If
provided together with NQSO, TEMPLATEID wins.
normfilter (str, optional): normalization filter
seed (int, optional): Seed for random number generator.
rand (numpy.RandomState, optional): RandomState object used for the
random number generation. If provided together with SEED, this
optional input superseeds the numpy.RandomState object instantiated by
SEED.
qso (desisim.templates.QSO, optional): object with which to generate
individual spectra/templates.
add_dlas (bool): Inject damped Lya systems into the Lya forest
These are done according to the current best estimates for the incidence dN/dz
(Prochaska et al. 2008, ApJ, 675, 1002)
Set in calc_lz
These are *not* inserted according to overdensity along the sightline
nocolorcuts (bool, optional): Do not apply the fiducial rzW1W2 color-cuts
cuts (default True).
Returns (flux, wave, meta, dla_meta) where:
* flux (numpy.ndarray): Array [nmodel, npix] of observed-frame spectra
(erg/s/cm2/A).
* wave (numpy.ndarray): Observed-frame [npix] wavelength array (Angstrom).
* meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum
with columns defined in desisim.io.empty_metatable *plus* RA, DEC.
* objmeta (astropy.Table): Table of additional object-specific meta-data
[nmodel] for each output spectrum with columns defined in
desisim.io.empty_metatable.
* dla_meta (astropy.Table): Table of meta-data [ndla] for the DLAs injected
into the spectra. Only returned if add_dlas=True
Note: `dla_meta` is only included if add_dlas=True.
"""
from scipy.interpolate import interp1d
import fitsio
from speclite.filters import load_filters
from desisim.templates import QSO
from desisim.io import empty_metatable
h = fitsio.FITS(lyafile)
if templateid is None:
if nqso is None:
nqso = len(h)-1
templateid = np.arange(nqso)
else:
templateid = np.asarray(templateid)
nqso = len(np.atleast_1d(templateid))
if rand is None:
rand = np.random.RandomState(seed)
templateseed = rand.randint(2**32, size=nqso)
#heads = [head.read_header() for head in h[templateid + 1]]
heads = []
for indx in templateid:
heads.append(h[indx + 1].read_header())
zqso = np.array([head['ZQSO'] for head in heads])
ra = np.array([head['RA'] for head in heads])
dec = np.array([head['DEC'] for head in heads])
mag_g = np.array([head['MAG_G'] for head in heads])
# Hard-coded filtername! Should match MAG_G!
normfilt = load_filters(normfilter)
if qso is None:
qso = QSO(normfilter_south=normfilter, wave=wave)
wave = qso.wave
flux = np.zeros([nqso, len(wave)], dtype='f4')
meta, objmeta = empty_metatable(objtype='QSO', nmodel=nqso)
meta['TEMPLATEID'][:] = templateid
meta['REDSHIFT'][:] = zqso
meta['MAG'][:] = mag_g
meta['MAGFILTER'][:] = normfilter
meta['SEED'][:] = templateseed
meta['RA'] = ra
meta['DEC'] = dec
# Lists for DLA meta data
if add_dlas:
dla_NHI, dla_z, dla_id = [], [], []
# Loop on quasars
for ii, indx in enumerate(templateid):
flux1, _, meta1, objmeta1 = qso.make_templates(nmodel=1, redshift=np.atleast_1d(zqso[ii]),
mag=np.atleast_1d(mag_g[ii]), seed=templateseed[ii],
nocolorcuts=nocolorcuts, lyaforest=False)
flux1 *= 1e-17
for col in meta1.colnames:
meta[col][ii] = meta1[col][0]
for col in objmeta1.colnames:
objmeta[col][ii] = objmeta1[col][0]
# read lambda and forest transmission
data = h[indx + 1].read()
la = data['LAMBDA'][:]
tr = data['FLUX'][:]
if len(tr):
# Interpolate the transmission at the spectral wavelengths, if
# outside the forest, the transmission is 1.
itr = interp1d(la, tr, bounds_error=False, fill_value=1.0)
flux1 *= itr(wave)
# Inject a DLA?
if add_dlas:
if np.min(wave/lambda_RF_LYA - 1) < zqso[ii]: # Any forest?
dlas, dla_model = insert_dlas(wave, zqso[ii], seed=templateseed[ii])
ndla = len(dlas)
if ndla > 0:
flux1 *= dla_model
# Meta
dla_z += [idla['z'] for idla in dlas]
dla_NHI += [idla['N'] for idla in dlas]
dla_id += [indx]*ndla
padflux, padwave = normfilt.pad_spectrum(flux1, wave, method='edge')
normmaggies = np.array(normfilt.get_ab_maggies(padflux, padwave,
mask_invalid=True)[normfilter])
factor = 10**(-0.4 * mag_g[ii]) / normmaggies
flux1 *= factor
for key in ('FLUX_G', 'FLUX_R', 'FLUX_Z', 'FLUX_W1', 'FLUX_W2'):
meta[key][ii] *= factor
flux[ii, :] = flux1[:]
h.close()
# Finish
if add_dlas:
ndla = len(dla_id)
if ndla > 0:
from astropy.table import Table
dla_meta = Table()
dla_meta['NHI'] = dla_NHI # log NHI values
dla_meta['z'] = dla_z
dla_meta['ID'] = dla_id
else:
dla_meta = None
return flux*1e17, wave, meta, objmeta, dla_meta
else:
return flux*1e17, wave, meta, objmeta