PyQSOFit is a flexible tool to decompose quasar spectra. It was originally translated from Yue Shen's IDL qsofit code, but not exactly the same. With this program, people can get the information of quasar continuum and emission lines, e.g., continuum slope, line FWHM, line dispersion, EW, peak, etc. The main flows of this code are:
For the details, you can read Shen et al. (2015), Guo et al. (2018) and Wang et al.(2018, in prepration)
To install this code, just download and install it with pip
:
git clone https://github.com/legolason/PyQSOFit
cd PyQSOFit
python -m pip install .
This code was tested Python 3.9. These packages are required, but should have been installed automatically during the step above:
These packages are optional. They are not installed automatically but must be installed for parameter uncertainty estimates:
If using pPXF to get the host component (use_ppxf=True
),
Eigenspectra of galaxy and quasar from Yip et al. (2004a), (2004b) and Optical/UV Fe II templates are included in this package. We suggest that install python through anaconda, which incorporates most packages you need.
wave_range
parameterlmfit.params
). The best way to calculate the error is using MCMC. In previous versions of the code, we used a Monte Carlo resampling method, which perturbs the flux based on the input error. Computing parameter uncertainties will take longer, but in practice we found that MCMC is not significantly slower than the resampling method and MCMC probably gives more reliable parameter uncertainties and covariances.Hengxiao Guo, hengxiaoguo AT gmail.com (SHAO)
Yue Shen, shenyue AT illinois.edu (UIUC)
Shu Wang, wangshukiaa AT pku.edu.cn (Seoul National University)
Colin J. Burke, colinjb2 AT illinois.edu (UIUC)
Wenke Ren, rwk AT mail.ustc.edu.cn (USTC)
If our code makes your life easier, it would be appreciated to cite us:
%matplotlib inline
import glob, os, sys, timeit
import matplotlib
import numpy as np
sys.path.append('../')
from pyqsofit.PyQSOFit import QSOFit
from astropy.io import fits
from astropy.table import Table
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")
QSOFit.set_mpl_style()
# Show the versions so we know what works
import astropy
import lmfit
import pyqsofit
print(astropy.__version__)
print(lmfit.__version__)
print(pyqsofit.__version__)
import emcee # optional, for MCMC
print(emcee.__version__)
5.3.1 1.2.1 2.0.0 3.1.4
print(pyqsofit.__path__)
['/Users/wenke/data/20230710CodeRefine/20230710PyQSOFit/pyqsofit']
Step 1: Set up the model input parameters
Firstly, run the script below to produce the line list file, qsopar.fits, containing lines and their constraints, which will be needed in the following fitting program. From this file, you can change some specific parameters to suit your requirements, e.g., fitting range, line width, tie line center, tie line sigma, etc. If you want to fit extra lines, please append it to corresponding complex. Note that our line wavelength and sigma in the list are in Ln scale, like Lnlambda, Lnsigma.
path_ex = '.' #os.path.join(pyqsofit.__path__[0], '..', 'example')
# create a header
hdr0 = fits.Header()
hdr0['Author'] = 'Hengxiao Guo'
primary_hdu = fits.PrimaryHDU(header=hdr0)
"""
In this table, we specify the priors / initial conditions and boundaries for the line fitting parameters.
"""
line_priors = np.rec.array([
(6564.61, 'Ha', 6400, 6800, 'Ha_br', 2, 0.0, 0.0, 1e10, 5e-3, 0.004, 0.05, 0.015, 0, 0, 0, 0.05, 1),
(6564.61, 'Ha', 6400, 6800, 'Ha_na', 1, 0.0, 0.0, 1e10, 1e-3, 5e-4, 0.00169, 0.01, 1, 1, 0, 0.002, 1),
(6549.85, 'Ha', 6400, 6800, 'NII6549', 1, 0.0, 0.0, 1e10, 1e-3, 2.3e-4, 0.00169, 5e-3, 1, 1, 1, 0.001, 1),
(6585.28, 'Ha', 6400, 6800, 'NII6585', 1, 0.0, 0.0, 1e10, 1e-3, 2.3e-4, 0.00169, 5e-3, 1, 1, 1, 0.003, 1),
(6718.29, 'Ha', 6400, 6800, 'SII6718', 1, 0.0, 0.0, 1e10, 1e-3, 2.3e-4, 0.00169, 5e-3, 1, 1, 2, 0.001, 1),
(6732.67, 'Ha', 6400, 6800, 'SII6732', 1, 0.0, 0.0, 1e10, 1e-3, 2.3e-4, 0.00169, 5e-3, 1, 1, 2, 0.001, 1),
(4862.68, 'Hb', 4640, 5100, 'Hb_br', 2, 0.0, 0.0, 1e10, 5e-3, 0.004, 0.05, 0.01, 0, 0, 0, 0.01, 1),
(4862.68, 'Hb', 4640, 5100, 'Hb_na', 1, 0.0, 0.0, 1e10, 1e-3, 2.3e-4, 0.00169, 0.01, 1, 1, 0, 0.002, 1),
(4960.30, 'Hb', 4640, 5100, 'OIII4959c', 1, 0.0, 0.0, 1e10, 1e-3, 2.3e-4, 0.00169, 0.01, 1, 1, 0, 0.002, 1),
(5008.24, 'Hb', 4640, 5100, 'OIII5007c', 1, 0.0, 0.0, 1e10, 1e-3, 2.3e-4, 0.00169, 0.01, 1, 1, 0, 0.004, 1),
(4960.30, 'Hb', 4640, 5100, 'OIII4959w', 1, 0.0, 0.0, 1e10, 3e-3, 2.3e-4, 0.004, 0.01, 2, 2, 0, 0.001, 1),
(5008.24, 'Hb', 4640, 5100, 'OIII5007w', 1, 0.0, 0.0, 1e10, 3e-3, 2.3e-4, 0.004, 0.01, 2, 2, 0, 0.002, 1),
#(4687.02, 'Hb', 4640, 5100, 'HeII4687_br', 1, 0.0, 0.0, 1e10, 5e-3, 0.004, 0.05, 0.005, 0, 0, 0, 0.001, 1),
#(4687.02, 'Hb', 4640, 5100, 'HeII4687_na', 1, 0.0, 0.0, 1e10, 1e-3, 2.3e-4, 0.00169, 0.005, 1, 1, 0, 0.001, 1),
#(3934.78, 'CaII', 3900, 3960, 'CaII3934' , 2, 0.0, 0.0, 1e10, 1e-3, 3.333e-4, 0.00169, 0.01, 99, 0, 0, -0.001, 1),
#(3728.48, 'OII', 3650, 3800, 'OII3728', 1, 0.0, 0.0, 1e10, 1e-3, 3.333e-4, 0.00169, 0.01, 1, 1, 0, 0.001, 1),
#(3426.84, 'NeV', 3380, 3480, 'NeV3426', 1, 0.0, 0.0, 1e10, 1e-3, 3.333e-4, 0.00169, 0.01, 0, 0, 0, 0.001, 1),
#(3426.84, 'NeV', 3380, 3480, 'NeV3426_br', 1, 0.0, 0.0, 1e10, 5e-3, 0.0025, 0.02, 0.01, 0, 0, 0, 0.001, 1),
(2798.75, 'MgII', 2700, 2900, 'MgII_br', 2, 0.0, 0.0, 1e10, 5e-3, 0.004, 0.05, 0.015, 0, 0, 0, 0.05, 1),
(2798.75, 'MgII', 2700, 2900, 'MgII_na', 1, 0.0, 0.0, 1e10, 1e-3, 5e-4, 0.00169, 0.01, 1, 1, 0, 0.002, 1),
(1908.73, 'CIII', 1700, 1970, 'CIII_br', 2, 0.0, 0.0, 1e10, 5e-3, 0.004, 0.05, 0.015, 99, 0, 0, 0.01, 1),
#(1908.73, 'CIII', 1700, 1970, 'CIII_na', 1, 0.0, 0.0, 1e10, 1e-3, 5e-4, 0.00169, 0.01, 1, 1, 0, 0.002, 1),
#(1892.03, 'CIII', 1700, 1970, 'SiIII1892', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.015, 0.003, 1, 1, 0, 0.005, 1),
#(1857.40, 'CIII', 1700, 1970, 'AlIII1857', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.015, 0.003, 1, 1, 0, 0.005, 1),
#(1816.98, 'CIII', 1700, 1970, 'SiII1816', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.015, 0.01, 1, 1, 0, 0.0002, 1),
#(1786.7, 'CIII', 1700, 1970, 'FeII1787', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.015, 0.01, 1, 1, 0, 0.0002, 1),
#(1750.26, 'CIII', 1700, 1970, 'NIII1750', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.015, 0.01, 1, 1, 0, 0.001, 1),
#(1718.55, 'CIII', 1700, 1900, 'NIV1718', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.015, 0.01, 1, 1, 0, 0.001, 1),
(1549.06, 'CIV', 1500, 1700, 'CIV_br', 2, 0.0, 0.0, 1e10, 5e-3, 0.004, 0.05, 0.015, 0, 0, 0, 0.05, 1),
# (1549.06, 'CIV', 1500, 1700, 'CIV_na', 1, 0.0, 0.0, 1e10, 1e-3, 5e-4, 0.00169, 0.01, 1, 1, 0, 0.002, 1),
#(1640.42, 'CIV', 1500, 1700, 'HeII1640', 1, 0.0, 0.0, 1e10, 1e-3, 5e-4, 0.00169, 0.008, 1, 1, 0, 0.002, 1),
#(1663.48, 'CIV', 1500, 1700, 'OIII1663', 1, 0.0, 0.0, 1e10, 1e-3, 5e-4, 0.00169, 0.008, 1, 1, 0, 0.002, 1),
#(1640.42, 'CIV', 1500, 1700, 'HeII1640_br', 1, 0.0, 0.0, 1e10, 5e-3, 0.0025, 0.02, 0.008, 1, 1, 0, 0.002, 1),
#(1663.48, 'CIV', 1500, 1700, 'OIII1663_br', 1, 0.0, 0.0, 1e10, 5e-3, 0.0025, 0.02, 0.008, 1, 1, 0, 0.002, 1),
#(1402.06, 'SiIV', 1290, 1450, 'SiIV_OIV1', 1, 0.0, 0.0, 1e10, 5e-3, 0.002, 0.05, 0.015, 1, 1, 0, 0.05, 1),
#(1396.76, 'SiIV', 1290, 1450, 'SiIV_OIV2', 1, 0.0, 0.0, 1e10, 5e-3, 0.002, 0.05, 0.015, 1, 1, 0, 0.05, 1),
#(1335.30, 'SiIV', 1290, 1450, 'CII1335', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.015, 0.01, 1, 1, 0, 0.001, 1),
#(1304.35, 'SiIV', 1290, 1450, 'OI1304', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.015, 0.01, 1, 1, 0, 0.001, 1),
(1215.67, 'Lya', 1150, 1290, 'Lya_br', 3, 0.0, 0.0, 1e10, 5e-3, 0.002, 0.05, 0.02, 0, 0, 0, 0.05, 1),
(1240.14, 'Lya', 1150, 1290, 'NV1240', 1, 0.0, 0.0, 1e10, 2e-3, 0.001, 0.01, 0.005, 0, 0, 0, 0.002, 1),
# (1215.67, 'Lya', 1150, 1290, 'Lya_na', 1, 0.0, 0.0, 1e10, 1e-3, 5e-4, 0.00169, 0.01, 0, 0, 0, 0.002, 1),
],
formats='float32, a20, float32, float32, a20, int32, float32, float32, float32, float32, float32, float32, float32, int32, int32, int32, float32, int32',
names=' lambda, compname, minwav, maxwav, linename, ngauss, inisca, minsca, maxsca, inisig, minsig, maxsig, voff, vindex, windex, findex, fvalue, vary')
# Header
hdr1 = fits.Header()
hdr1['lambda'] = 'Vacuum Wavelength in Ang'
hdr1['minwav'] = 'Lower complex fitting wavelength range'
hdr1['maxwav'] = 'Upper complex fitting wavelength range'
hdr1['ngauss'] = 'Number of Gaussians for the line'
# Can be set to negative for absorption lines if you want
hdr1['inisca'] = 'Initial guess of line scale [flux]'
hdr1['minsca'] = 'Lower range of line scale [flux]'
hdr1['maxsca'] = 'Upper range of line scale [flux]'
hdr1['inisig'] = 'Initial guess of linesigma [lnlambda]'
hdr1['minsig'] = 'Lower range of line sigma [lnlambda]'
hdr1['maxsig'] = 'Upper range of line sigma [lnlambda]'
hdr1['voff '] = 'Limits on velocity offset from the central wavelength [lnlambda]'
hdr1['vindex'] = 'Entries w/ same NONZERO vindex constrained to have same velocity'
hdr1['windex'] = 'Entries w/ same NONZERO windex constrained to have same width'
hdr1['findex'] = 'Entries w/ same NONZERO findex have constrained flux ratios'
hdr1['fvalue'] = 'Relative scale factor for entries w/ same findex'
hdr1['vary'] = 'Whether or not to vary the parameter (set to 0 to fix the line parameter to initial values)'
# Save line info
hdu1 = fits.BinTableHDU(data=line_priors, header=hdr1, name='line_priors')
"""
In this table, we specify the windows and priors / initial conditions and boundaries for the continuum fitting parameters.
"""
conti_windows = np.rec.array([
(1150., 1170.),
(1275., 1290.),
(1350., 1360.),
(1445., 1465.),
(1690., 1705.),
(1770., 1810.),
(1970., 2400.),
(2480., 2675.),
(2925., 3400.),
(3775., 3832.),
(4000., 4050.),
(4200., 4230.),
(4435., 4640.),
(5100., 5535.),
(6005., 6035.),
(6110., 6250.),
(6800., 7000.),
(7160., 7180.),
(7500., 7800.),
(8050., 8150.), # Continuum fitting windows (to avoid emission line, etc.) [AA]
],
formats = 'float32, float32',
names = 'min, max')
hdu2 = fits.BinTableHDU(data=conti_windows, name='conti_windows')
conti_priors = np.rec.array([
('Fe_uv_norm', 0.0, 0.0, 1e10, 1), # Normalization of the MgII Fe template [flux]
('Fe_uv_FWHM', 3000, 1200, 18000, 1), # FWHM of the MgII Fe template [AA]
('Fe_uv_shift', 0.0, -0.01, 0.01, 1), # Wavelength shift of the MgII Fe template [lnlambda]
('Fe_op_norm', 0.0, 0.0, 1e10, 1), # Normalization of the Hbeta/Halpha Fe template [flux]
('Fe_op_FWHM', 3000, 1200, 18000, 1), # FWHM of the Hbeta/Halpha Fe template [AA]
('Fe_op_shift', 0.0, -0.01, 0.01, 1), # Wavelength shift of the Hbeta/Halpha Fe template [lnlambda]
('PL_norm', 1.0, 0.0, 1e10, 1), # Normalization of the power-law (PL) continuum f_lambda = (lambda/3000)^-alpha
('PL_slope', -1.5, -5.0, 3.0, 1), # Slope of the power-law (PL) continuum
('Blamer_norm', 0.0, 0.0, 1e10, 1), # Normalization of the Balmer continuum at < 3646 AA [flux] (Dietrich et al. 2002)
('Balmer_Te', 15000, 10000, 50000, 1), # Te of the Balmer continuum at < 3646 AA [K?]
('Balmer_Tau', 0.5, 0.1, 2.0, 1), # Tau of the Balmer continuum at < 3646 AA
('conti_a_0', 0.0, None, None, 1), # 1st coefficient of the polynomial continuum
('conti_a_1', 0.0, None, None, 1), # 2nd coefficient of the polynomial continuum
('conti_a_2', 0.0, None, None, 1), # 3rd coefficient of the polynomial continuum
# Note: The min/max bounds on the conti_a_0 coefficients are ignored by the code,
# so they can be determined automatically for numerical stability.
],
formats = 'a20, float32, float32, float32, int32',
names = 'parname, initial, min, max, vary')
hdr3 = fits.Header()
hdr3['ini'] = 'Initial guess of line scale [flux]'
hdr3['min'] = 'FWHM of the MgII Fe template'
hdr3['max'] = 'Wavelength shift of the MgII Fe template'
hdr3['vary'] = 'Whether or not to vary the parameter (set to 0 to fix the continuum parameter to initial values)'
hdu3 = fits.BinTableHDU(data=conti_priors, header=hdr3, name='conti_priors')
"""
In this table, we allow user to customized some key parameters in our result measurements.
"""
measure_info = Table(
[
[[1350, 1450, 3000, 4200, 5100]],
[[
# [2240, 2650],
[4435, 4685],
]]
],
names=([
'cont_loc',
'Fe_flux_range'
]),
dtype=([
'float32',
'float32'
])
)
hdr4 = fits.Header()
hdr4['cont_loc'] = 'The wavelength of continuum luminosity in results'
hdr4['Fe_flux_range'] = 'Fe emission wavelength range calculated in results'
hdu4 = fits.BinTableHDU(data=measure_info, header=hdr4, name='measure_info')
hdu_list = fits.HDUList([primary_hdu, hdu1, hdu2, hdu3, hdu4])
hdu_list.writeto(os.path.join(path_ex, 'qsopar.fits'), overwrite=True)
Print the table:
Table(line_priors)
lambda | compname | minwav | maxwav | linename | ngauss | inisca | minsca | maxsca | inisig | minsig | maxsig | voff | vindex | windex | findex | fvalue | vary |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
float32 | bytes20 | float32 | float32 | bytes20 | int32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | int32 | int32 | int32 | float32 | int32 |
6564.61 | Ha | 6400.0 | 6800.0 | Ha_br | 2 | 0.0 | 0.0 | 10000000000.0 | 0.005 | 0.004 | 0.05 | 0.015 | 0 | 0 | 0 | 0.05 | 1 |
6564.61 | Ha | 6400.0 | 6800.0 | Ha_na | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.0005 | 0.00169 | 0.01 | 1 | 1 | 0 | 0.002 | 1 |
6549.85 | Ha | 6400.0 | 6800.0 | NII6549 | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.00023 | 0.00169 | 0.005 | 1 | 1 | 1 | 0.001 | 1 |
6585.28 | Ha | 6400.0 | 6800.0 | NII6585 | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.00023 | 0.00169 | 0.005 | 1 | 1 | 1 | 0.003 | 1 |
6718.29 | Ha | 6400.0 | 6800.0 | SII6718 | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.00023 | 0.00169 | 0.005 | 1 | 1 | 2 | 0.001 | 1 |
6732.67 | Ha | 6400.0 | 6800.0 | SII6732 | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.00023 | 0.00169 | 0.005 | 1 | 1 | 2 | 0.001 | 1 |
4862.68 | Hb | 4640.0 | 5100.0 | Hb_br | 2 | 0.0 | 0.0 | 10000000000.0 | 0.005 | 0.004 | 0.05 | 0.01 | 0 | 0 | 0 | 0.01 | 1 |
4862.68 | Hb | 4640.0 | 5100.0 | Hb_na | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.00023 | 0.00169 | 0.01 | 1 | 1 | 0 | 0.002 | 1 |
4960.3 | Hb | 4640.0 | 5100.0 | OIII4959c | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.00023 | 0.00169 | 0.01 | 1 | 1 | 0 | 0.002 | 1 |
5008.24 | Hb | 4640.0 | 5100.0 | OIII5007c | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.00023 | 0.00169 | 0.01 | 1 | 1 | 0 | 0.004 | 1 |
4960.3 | Hb | 4640.0 | 5100.0 | OIII4959w | 1 | 0.0 | 0.0 | 10000000000.0 | 0.003 | 0.00023 | 0.004 | 0.01 | 2 | 2 | 0 | 0.001 | 1 |
5008.24 | Hb | 4640.0 | 5100.0 | OIII5007w | 1 | 0.0 | 0.0 | 10000000000.0 | 0.003 | 0.00023 | 0.004 | 0.01 | 2 | 2 | 0 | 0.002 | 1 |
2798.75 | MgII | 2700.0 | 2900.0 | MgII_br | 2 | 0.0 | 0.0 | 10000000000.0 | 0.005 | 0.004 | 0.05 | 0.015 | 0 | 0 | 0 | 0.05 | 1 |
2798.75 | MgII | 2700.0 | 2900.0 | MgII_na | 1 | 0.0 | 0.0 | 10000000000.0 | 0.001 | 0.0005 | 0.00169 | 0.01 | 1 | 1 | 0 | 0.002 | 1 |
1908.73 | CIII | 1700.0 | 1970.0 | CIII_br | 2 | 0.0 | 0.0 | 10000000000.0 | 0.005 | 0.004 | 0.05 | 0.015 | 99 | 0 | 0 | 0.01 | 1 |
1549.06 | CIV | 1500.0 | 1700.0 | CIV_br | 2 | 0.0 | 0.0 | 10000000000.0 | 0.005 | 0.004 | 0.05 | 0.015 | 0 | 0 | 0 | 0.05 | 1 |
1215.67 | Lya | 1150.0 | 1290.0 | Lya_br | 3 | 0.0 | 0.0 | 10000000000.0 | 0.005 | 0.002 | 0.05 | 0.02 | 0 | 0 | 0 | 0.05 | 1 |
1240.14 | Lya | 1150.0 | 1290.0 | NV1240 | 1 | 0.0 | 0.0 | 10000000000.0 | 0.002 | 0.001 | 0.01 | 0.005 | 0 | 0 | 0 | 0.002 | 1 |
Note: if you want to tie the line properties in different complex, you can enlarge the complex range.
Step 2: Read the spectrum
Setup the paths and read in your spectrum. Our code is written under the frame of SDSS spectral data format. Other data is also available as long as they include wavelength, flux, error, and redshift, and make sure the wavelength resolution is the same as SDSS spectrum (For SDSS the pixel scale is 1.e-4 in log space).
path_out = os.path.join(pyqsofit.__path__[0], '../', 'example/data/')
# Requried
data = fits.open(os.path.join(path_ex, 'data/spec-0332-52367-0639.fits'))
lam = 10 ** data[1].data['loglam'] # OBS wavelength [A]
flux = data[1].data['flux'] # OBS flux [erg/s/cm^2/A]
err = 1 / np.sqrt(data[1].data['ivar']) # 1 sigma error
z = data[2].data['z'][0] # Redshift
# Optional
ra = data[0].header['plug_ra'] # RA
dec = data[0].header['plug_dec'] # DEC
plateid = data[0].header['plateid'] # SDSS plate ID
mjd = data[0].header['mjd'] # SDSS MJD
fiberid = data[0].header['fiberid'] # SDSS fiber ID
Step 3: Fit the spectrum
Use QSOFit to input the lam, flux, err, z, and other optinal parameters. Use function Fit to perform the fitting. Default settings cannot meet all needs. Please change settings for your own requirements. It depends on what science you need. The following example set dereddening, host decomposition to True.
The broad_fwhm
parameter can be adjusted depending on your definition (default is 1200 km s$^{-1}$).
# Prepare data
q_mle = QSOFit(lam, flux, err, z, ra=ra, dec=dec, plateid=plateid, mjd=mjd, fiberid=fiberid, path=path_ex)
# Double check the installation path with the PCA / Fe template files
# print('install path:', q_mle.install_path)
# Change it if you installed them somewhere else
#q_mle.install_path = '...'
start = timeit.default_timer()
# Do the fitting
q_mle.Fit(name=None, # customize the name of given targets. Default: plate-mjd-fiber
# prepocessing parameters
nsmooth=1, # do n-pixel smoothing to the raw input flux and err spectra
and_mask=False, # delete the and masked pixels
or_mask=False, # delete the or masked pixels
reject_badpix=True, # reject 10 most possible outliers by the test of pointDistGESD
deredden=True, # correct the Galactic extinction
wave_range=None, # trim input wavelength
wave_mask=None, # 2-D array, mask the given range(s)
# host decomposition parameters
decompose_host=True, # If True, the host galaxy-QSO decomposition will be applied
host_prior=True, # If True, the code will adopt prior-informed method to assist decomposition. Currently, only 'CZBIN1' and 'DZBIN1' model for QSO PCA are available. And the model for galaxy must be PCA too.
host_prior_scale=0.2, # scale of prior panelty. Usually, 0.2 works fine for SDSS spectra. Adjust it smaller if you find the prior affect the fitting results too much.
host_line_mask=True, # If True, the line region of galaxy will be masked when subtracted from original spectra.
decomp_na_mask=True, # If True, the narrow line region will be masked when perform decomposition
qso_type='CZBIN1', # PCA template name for quasar
npca_qso=10, # numebr of quasar templates
host_type='PCA', # template name for galaxy
npca_gal=5, # number of galaxy templates
# continuum model fit parameters
Fe_uv_op=True, # If True, fit continuum with UV and optical FeII template
poly=True, # If True, fit continuum with the polynomial component to account for the dust reddening
BC=False, # If True, fit continuum with Balmer continua from 1000 to 3646A
initial_guess=None, # Initial parameters for continuum model, read the annotation of this function for detail
rej_abs_conti=False, # If True, it will iterately reject 3 sigma outlier absorption pixels in the continuum
n_pix_min_conti=100, # Minimum number of negative pixels for host continuuum fit to be rejected.
# emission line fit parameters
linefit=True, # If True, the emission line will be fitted
rej_abs_line=False,
# If True, it will iterately reject 3 sigma outlier absorption pixels in the emission lines
# fitting method selection
MC=False,
# If True, do Monte Carlo resampling of the spectrum based on the input error array to produce the MC error array
MCMC=False,
# If True, do Markov Chain Monte Carlo sampling of the posterior probability densities to produce the error array
nsamp=200,
# The number of trials of the MC process (if MC=True) or number samples to run MCMC chain (if MCMC=True)
# advanced fitting parameters
param_file_name='qsopar.fits', # Name of the qso fitting parameter FITS file.
nburn=20, # The number of burn-in samples to run MCMC chain
nthin=10, # To set the MCMC chain returns every n samples
epsilon_jitter=0.,
# Initial jitter for every initial guass to avoid local minimum. (Under test, not recommanded to change)
# customize the results
save_result=False, # If True, all the fitting results will be saved to a fits file
save_fits_name=None, # The output name of the result fits
save_fits_path=path_out, # The output path of the result fits
plot_fig=True, # If True, the fitting results will be plotted
save_fig=False, # If True, the figure will be saved
plot_corner=True, # Whether or not to plot the corner plot results if MCMC=True
# debugging mode
verbose=True, # turn on (True) or off (False) debugging output
# sublevel parameters for figure plot and emcee
kwargs_plot={
'save_fig_path': '.', # The output path of the figure
'broad_fwhm' : 1200 # km/s, lower limit that code decide if a line component belongs to broad component
},
kwargs_conti_emcee={},
kwargs_line_emcee={})
end = timeit.default_timer()
print(f'Fitting finished in {np.round(end - start, 1)}s')
Name Value Min Max Stderr Vary Expr Brute_Step Balmer_Tau 0.5 0.1 2 None False None None Balmer_Te 1.5e+04 1e+04 5e+04 None False None None Blamer_norm 0 0 1e+10 None False None None Fe_op_FWHM 3000 1200 1.8e+04 None True None None Fe_op_norm 0 0 1e+10 None True None None Fe_op_shift 0 -0.01 0.01 None True None None Fe_uv_FWHM 3000 1200 1.8e+04 None False None None Fe_uv_norm 0 0 1e+10 None False None None Fe_uv_shift 0 -0.01 0.01 None False None None PL_norm 1 0 1e+10 None True None None PL_slope -1.5 -5 3 None True None None conti_a_0 0 -inf inf None True None None conti_a_1 0 -inf inf None True None None conti_a_2 0 -inf inf None True None None Fitting continuum Fit report [[Variables]] Fe_uv_norm: 0 (fixed) Fe_uv_FWHM: 3000 (fixed) Fe_uv_shift: 0 (fixed) Fe_op_norm: 3.87007482 (init = 0) Fe_op_FWHM: 18000.0000 (init = 3000) Fe_op_shift: 0.00999996 (init = 0) PL_norm: 35.8934976 (init = 1) PL_slope: -1.94745412 (init = -1.5) Blamer_norm: 0 (fixed) Balmer_Te: 15000 (fixed) Balmer_Tau: 0.5 (fixed) conti_a_0: -4123.34586 (init = 0) conti_a_1: 2.72341681 (init = 0) conti_a_2: -3.4858e-04 (init = 0) Name Value Min Max Stderr Vary Expr Brute_Step Hb_br_1_dwave 0 -0.01 0.01 None True None None Hb_br_1_scale 0 0 1e+10 None True None None Hb_br_1_sigma 0.005 0.004 0.05 None True None None Hb_br_2_dwave 0 -0.01 0.01 None True None None Hb_br_2_scale 0 0 1e+10 None True None None Hb_br_2_sigma 0.005 0.004 0.05 None True None None Hb_na_1_dwave 0 -0.01 0.01 None True None None Hb_na_1_scale 0 0 1e+10 None True None None Hb_na_1_sigma 0.001 0.00023 0.00169 None True None None OIII4959c_1_dwave 0 -0.01 0.01 None False Hb_na_1_dwave None OIII4959c_1_scale 0 0 1e+10 None True None None OIII4959c_1_sigma 0.001 0.00023 0.00169 None False Hb_na_1_sigma None OIII4959w_1_dwave 0 -0.01 0.01 None True None None OIII4959w_1_scale 0 0 1e+10 None True None None OIII4959w_1_sigma 0.003 0.00023 0.004 None True None None OIII5007c_1_dwave 0 -0.01 0.01 None False Hb_na_1_dwave None OIII5007c_1_scale 0 0 1e+10 None True None None OIII5007c_1_sigma 0.001 0.00023 0.00169 None False Hb_na_1_sigma None OIII5007w_1_dwave 0 -0.01 0.01 None False OIII4959w_1_dwave None OIII5007w_1_scale 0 0 1e+10 None True None None OIII5007w_1_sigma 0.003 0.00023 0.004 None False OIII4959w_1_sigma None Fitting complex Hb Fit report [[Variables]] Hb_br_1_scale: 7.32474748 +/- 0.61572363 (8.41%) (init = 0) Hb_br_1_dwave: -1.4564e-04 +/- 4.5392e-04 (311.66%) (init = 0) Hb_br_1_sigma: 0.00629871 +/- 6.7530e-04 (10.72%) (init = 0.005) Hb_br_2_scale: 4.00759370 +/- 0.56651017 (14.14%) (init = 0) Hb_br_2_dwave: 0.00682185 +/- 0.00185862 (27.25%) (init = 0) Hb_br_2_sigma: 0.02018566 +/- 0.00199827 (9.90%) (init = 0.005) Hb_na_1_scale: 8.83420725 +/- 1.23806683 (14.01%) (init = 0) Hb_na_1_dwave: -5.1273e-04 +/- 2.5470e-05 (4.97%) (init = 0) Hb_na_1_sigma: 4.6127e-04 +/- 1.9917e-05 (4.32%) (init = 0.001) OIII4959c_1_scale: 30.2192860 +/- 1.43132809 (4.74%) (init = 0) OIII4959c_1_dwave: -5.1273e-04 +/- 2.5470e-05 (4.97%) == 'Hb_na_1_dwave' OIII4959c_1_sigma: 4.6127e-04 +/- 1.9917e-05 (4.32%) == 'Hb_na_1_sigma' OIII5007c_1_scale: 95.2483925 +/- 3.04570984 (3.20%) (init = 0) OIII5007c_1_dwave: -5.1273e-04 +/- 2.5470e-05 (4.97%) == 'Hb_na_1_dwave' OIII5007c_1_sigma: 4.6127e-04 +/- 1.9917e-05 (4.32%) == 'Hb_na_1_sigma' OIII4959w_1_scale: 16.6302971 +/- 1.57002025 (9.44%) (init = 0) OIII4959w_1_dwave: 6.6893e-04 +/- 3.4721e-05 (5.19%) (init = 0) OIII4959w_1_sigma: 3.6380e-04 +/- 2.6231e-05 (7.21%) (init = 0.003) OIII5007w_1_scale: 54.6045026 +/- 2.68006899 (4.91%) (init = 0) OIII5007w_1_dwave: 6.6893e-04 +/- 3.4721e-05 (5.19%) == 'OIII4959w_1_dwave' OIII5007w_1_sigma: 3.6380e-04 +/- 2.6231e-05 (7.21%) == 'OIII4959w_1_sigma' [[Correlations]] (unreported correlations are < 0.100) C(Hb_br_2_scale, Hb_br_2_sigma) = -0.8144 C(Hb_na_1_dwave, OIII4959w_1_dwave) = +0.8069 C(Hb_br_1_sigma, Hb_br_2_scale) = -0.7858 C(Hb_na_1_sigma, OIII4959w_1_dwave) = +0.7724 C(Hb_br_1_scale, Hb_br_2_scale) = -0.7642 C(OIII4959w_1_dwave, OIII4959w_1_sigma) = -0.7597 C(Hb_na_1_dwave, Hb_na_1_sigma) = +0.7516 C(Hb_na_1_dwave, OIII4959w_1_sigma) = -0.7393 C(Hb_na_1_sigma, OIII4959w_1_sigma) = -0.6535 C(Hb_br_1_scale, Hb_br_2_sigma) = +0.6344 C(Hb_br_2_scale, Hb_br_2_dwave) = -0.5485 C(Hb_br_1_sigma, Hb_br_2_sigma) = +0.5457 C(Hb_na_1_sigma, OIII5007c_1_scale) = -0.5297 C(Hb_br_1_scale, Hb_br_2_dwave) = +0.5283 C(Hb_br_1_sigma, Hb_br_2_dwave) = +0.5028 C(Hb_br_2_dwave, Hb_br_2_sigma) = +0.4451 C(Hb_na_1_dwave, OIII5007c_1_scale) = -0.3908 C(Hb_br_1_scale, Hb_br_1_sigma) = +0.3739 C(Hb_br_1_dwave, Hb_br_2_sigma) = +0.3733 C(Hb_br_1_dwave, Hb_br_2_scale) = -0.3278 C(Hb_na_1_sigma, OIII5007w_1_scale) = -0.2764 C(Hb_br_1_sigma, Hb_na_1_scale) = +0.2604 C(OIII5007c_1_scale, OIII4959w_1_dwave) = -0.2394 C(OIII5007c_1_scale, OIII5007w_1_scale) = +0.2298 C(Hb_br_2_scale, OIII4959c_1_scale) = -0.2281 C(Hb_na_1_sigma, OIII4959w_1_scale) = -0.2220 C(Hb_na_1_dwave, OIII4959w_1_scale) = -0.2144 C(Hb_br_1_dwave, Hb_br_1_sigma) = +0.2022 C(OIII4959w_1_sigma, OIII5007w_1_scale) = -0.1975 C(OIII5007c_1_scale, OIII4959w_1_scale) = +0.1935 C(Hb_br_1_scale, Hb_br_1_dwave) = +0.1917 C(OIII4959w_1_scale, OIII5007w_1_scale) = +0.1856 C(Hb_na_1_dwave, OIII5007w_1_scale) = -0.1839 C(Hb_br_1_scale, Hb_na_1_scale) = -0.1802 C(Hb_br_2_scale, OIII4959w_1_scale) = -0.1781 C(Hb_br_1_scale, OIII4959c_1_scale) = +0.1778 C(Hb_br_2_sigma, OIII4959c_1_scale) = +0.1672 C(OIII4959w_1_scale, OIII4959w_1_dwave) = -0.1654 C(Hb_br_1_sigma, OIII4959c_1_scale) = +0.1540 C(OIII4959c_1_scale, OIII5007c_1_scale) = +0.1488 C(Hb_br_1_dwave, Hb_br_2_dwave) = -0.1403 C(Hb_br_1_scale, OIII4959w_1_scale) = +0.1391 C(OIII5007c_1_scale, OIII4959w_1_sigma) = +0.1385 C(Hb_na_1_sigma, OIII4959c_1_scale) = -0.1364 C(OIII4959c_1_scale, OIII5007w_1_scale) = +0.1247 C(Hb_br_2_sigma, OIII4959w_1_scale) = +0.1238 C(Hb_br_1_sigma, OIII4959w_1_scale) = +0.1235 C(Hb_br_2_scale, Hb_na_1_scale) = -0.1209 Name Value Min Max Stderr Vary Expr Brute_Step Ha_br_1_dwave 0 -0.015 0.015 None True None None Ha_br_1_scale 0 0 1e+10 None True None None Ha_br_1_sigma 0.005 0.004 0.05 None True None None Ha_br_2_dwave 0 -0.015 0.015 None True None None Ha_br_2_scale 0 0 1e+10 None True None None Ha_br_2_sigma 0.005 0.004 0.05 None True None None Ha_na_1_dwave 0 -0.01 0.01 None True None None Ha_na_1_scale 0 0 1e+10 None True None None Ha_na_1_sigma 0.001 0.0005 0.00169 None True None None NII6549_1_dwave 0 -0.005 0.005 None False Ha_na_1_dwave None NII6549_1_scale 0 0 1e+10 None True None None NII6549_1_sigma 0.001 0.00023 0.00169 None False Ha_na_1_sigma None NII6585_1_dwave 0 -0.005 0.005 None False Ha_na_1_dwave None NII6585_1_scale 0 0 1e+10 None False 3.0 * NII6549_1_scale None NII6585_1_sigma 0.001 0.00023 0.00169 None False Ha_na_1_sigma None SII6718_1_dwave 0 -0.005 0.005 None False Ha_na_1_dwave None SII6718_1_scale 0 0 1e+10 None True None None SII6718_1_sigma 0.001 0.00023 0.00169 None False Ha_na_1_sigma None SII6732_1_dwave 0 -0.005 0.005 None False Ha_na_1_dwave None SII6732_1_scale 0 0 1e+10 None False 1.0 * SII6718_1_scale None SII6732_1_sigma 0.001 0.00023 0.00169 None False Ha_na_1_sigma None Fitting complex Ha Fit report [[Variables]] Ha_br_1_scale: 17.5528986 +/- 2.26291140 (12.89%) (init = 0) Ha_br_1_dwave: -8.4732e-04 +/- 3.6979e-04 (43.64%) (init = 0) Ha_br_1_sigma: 0.00984311 +/- 4.1604e-04 (4.23%) (init = 0.005) Ha_br_2_scale: 16.9729131 +/- 2.11344416 (12.45%) (init = 0) Ha_br_2_dwave: 0.00150989 +/- 2.3588e-04 (15.62%) (init = 0) Ha_br_2_sigma: 0.00511150 +/- 4.4511e-04 (8.71%) (init = 0.005) Ha_na_1_scale: 22.0197149 +/- 1.38708356 (6.30%) (init = 0) Ha_na_1_dwave: -2.2102e-04 +/- 2.4586e-05 (11.12%) (init = 0) Ha_na_1_sigma: 6.6333e-04 +/- 2.5643e-05 (3.87%) (init = 0.001) NII6549_1_scale: 9.61714097 +/- 0.46326451 (4.82%) (init = 0) NII6549_1_dwave: -2.2102e-04 +/- 2.4586e-05 (11.12%) == 'Ha_na_1_dwave' NII6549_1_sigma: 6.6333e-04 +/- 2.5643e-05 (3.87%) == 'Ha_na_1_sigma' NII6585_1_scale: 28.8514229 +/- 1.38979352 (4.82%) == '3.0 * NII6549_1_scale' NII6585_1_dwave: -2.2102e-04 +/- 2.4586e-05 (11.12%) == 'Ha_na_1_dwave' NII6585_1_sigma: 6.6333e-04 +/- 2.5643e-05 (3.87%) == 'Ha_na_1_sigma' SII6718_1_scale: 11.4394444 +/- 0.51053274 (4.46%) (init = 0) SII6718_1_dwave: -2.2102e-04 +/- 2.4586e-05 (11.12%) == 'Ha_na_1_dwave' SII6718_1_sigma: 6.6333e-04 +/- 2.5643e-05 (3.87%) == 'Ha_na_1_sigma' SII6732_1_scale: 11.4394444 +/- 0.51053273 (4.46%) == '1.0 * SII6718_1_scale' SII6732_1_dwave: -2.2102e-04 +/- 2.4586e-05 (11.12%) == 'Ha_na_1_dwave' SII6732_1_sigma: 6.6333e-04 +/- 2.5643e-05 (3.87%) == 'Ha_na_1_sigma' [[Correlations]] (unreported correlations are < 0.100) C(Ha_br_1_scale, Ha_br_2_scale) = -0.9125 C(Ha_br_1_scale, Ha_br_1_sigma) = -0.9087 C(Ha_br_1_sigma, Ha_br_2_scale) = +0.8853 C(Ha_br_1_scale, Ha_br_2_sigma) = -0.8818 C(Ha_br_1_scale, Ha_br_1_dwave) = +0.7680 C(Ha_br_1_dwave, Ha_br_2_sigma) = -0.7625 C(Ha_br_1_sigma, Ha_br_2_sigma) = +0.6776 C(Ha_br_1_dwave, Ha_br_2_scale) = -0.6707 C(Ha_br_2_scale, Ha_br_2_sigma) = +0.6528 C(Ha_br_1_dwave, Ha_br_1_sigma) = -0.5767 C(Ha_br_1_sigma, Ha_br_2_dwave) = -0.5542 C(Ha_na_1_scale, NII6549_1_scale) = +0.4817 C(Ha_na_1_sigma, SII6718_1_scale) = -0.4796 C(Ha_br_2_scale, Ha_br_2_dwave) = -0.4712 C(Ha_br_1_scale, Ha_br_2_dwave) = +0.4004 C(Ha_br_2_sigma, Ha_na_1_sigma) = +0.3883 C(Ha_br_2_sigma, Ha_na_1_scale) = +0.3876 C(Ha_br_2_sigma, NII6549_1_scale) = +0.3774 C(Ha_br_1_dwave, Ha_na_1_scale) = -0.2612 C(Ha_br_1_dwave, Ha_na_1_sigma) = -0.2517 C(Ha_br_2_dwave, Ha_na_1_scale) = +0.2433 C(Ha_na_1_scale, Ha_na_1_sigma) = +0.2280 C(Ha_br_1_dwave, NII6549_1_scale) = -0.2104 C(Ha_br_1_scale, NII6549_1_scale) = -0.1943 C(Ha_br_2_dwave, Ha_na_1_sigma) = +0.1858 C(Ha_br_2_sigma, SII6718_1_scale) = -0.1845 C(Ha_br_2_dwave, Ha_br_2_sigma) = -0.1803 C(Ha_br_1_scale, Ha_na_1_scale) = -0.1765 C(Ha_br_1_scale, Ha_na_1_sigma) = -0.1734 C(Ha_br_2_dwave, Ha_na_1_dwave) = -0.1729 C(Ha_br_1_sigma, SII6718_1_scale) = -0.1705 C(Ha_na_1_sigma, NII6549_1_scale) = +0.1617 C(Ha_br_1_scale, SII6718_1_scale) = +0.1445 C(Ha_br_1_dwave, Ha_br_2_dwave) = -0.1333 C(Ha_br_1_sigma, NII6549_1_scale) = +0.1024 C(Ha_br_2_scale, Ha_na_1_sigma) = -0.1006 Fitting finished in 12.4s
The gray shaded bars at the top are the continuum windows used in the fitting.
Now you are already done with the QSO fitting part!
Step 3(optional): Compute parameter unertaintines using MC resampling:
For comparison, let's fit the same quasar spectrum, but this time using MC resampling of the spectrum based on the spectrum error. Obviously, this will take significantly longer depending on the number of samples you want to run nsamp
.
q_mc = QSOFit(lam, flux, err, z, ra=ra, dec=dec, plateid=plateid, mjd=mjd, fiberid=fiberid, path=path_ex)
start = timeit.default_timer()
# Do the fitting
q_mc.Fit(name=None, nsmooth=1, deredden=True, reject_badpix=False, wave_range=None, \
wave_mask=None, decompose_host=True, host_prior=True, decomp_na_mask=True, npca_gal=5, npca_qso=10, qso_type='CZBIN1',\
Fe_uv_op=True, poly=True, rej_abs_conti=False, rej_abs_line=True, MC=True, nsamp=200, linefit=True, \
save_result=True, kwargs_plot={'save_fig_path': '.'}, save_fits_name=None, verbose=True)
end = timeit.default_timer()
print(f'Fitting finished in {np.round(end - start, 1)}s')
Name Value Min Max Stderr Vary Expr Brute_Step Balmer_Tau 0.5 0.1 2 None False None None Balmer_Te 1.5e+04 1e+04 5e+04 None False None None Blamer_norm 0 0 1e+10 None False None None Fe_op_FWHM 3000 1200 1.8e+04 None True None None Fe_op_norm 0 0 1e+10 None True None None Fe_op_shift 0 -0.01 0.01 None True None None Fe_uv_FWHM 3000 1200 1.8e+04 None False None None Fe_uv_norm 0 0 1e+10 None False None None Fe_uv_shift 0 -0.01 0.01 None False None None PL_norm 1 0 1e+10 None True None None PL_slope -1.5 -5 3 None True None None conti_a_0 0 -inf inf None True None None conti_a_1 0 -inf inf None True None None conti_a_2 0 -inf inf None True None None Fitting continuum Fit report [[Variables]] Fe_uv_norm: 0 (fixed) Fe_uv_FWHM: 3000 (fixed) Fe_uv_shift: 0 (fixed) Fe_op_norm: 3.46537687 +/- 0.69362124 (20.02%) (init = 0) Fe_op_FWHM: 16514.2687 +/- 2414.29619 (14.62%) (init = 3000) Fe_op_shift: 0.00999970 +/- 0.00248865 (24.89%) (init = 0) PL_norm: 35.9833241 +/- 4.71896116 (13.11%) (init = 1) PL_slope: -2.17651038 +/- 1.46981162 (67.53%) (init = -1.5) Blamer_norm: 0 (fixed) Balmer_Te: 15000 (fixed) Balmer_Tau: 0.5 (fixed) conti_a_0: -1418.26808 +/- 7803.23636 (550.19%) (init = 0) conti_a_1: 1.54024167 +/- 2.73797425 (177.76%) (init = 0) conti_a_2: -2.0215e-04 +/- 2.9635e-04 (146.60%) (init = 0) [[Correlations]] (unreported correlations are < 0.100) C(conti_a_1, conti_a_2) = -0.9895 C(conti_a_0, conti_a_1) = -0.9844 C(PL_norm, PL_slope) = -0.9643 C(PL_slope, conti_a_0) = -0.9506 C(conti_a_0, conti_a_2) = +0.9488 C(Fe_op_norm, Fe_op_FWHM) = +0.9223 C(PL_slope, conti_a_1) = +0.8816 C(PL_norm, conti_a_0) = +0.8527 C(PL_slope, conti_a_2) = -0.8051 C(PL_norm, conti_a_1) = -0.7547 C(PL_norm, conti_a_2) = +0.6572 C(Fe_op_norm, PL_norm) = +0.4294 C(Fe_op_FWHM, PL_norm) = +0.4115 C(Fe_op_shift, conti_a_2) = +0.4047 C(Fe_op_shift, conti_a_1) = -0.3867 C(Fe_op_shift, conti_a_0) = +0.3472 C(Fe_op_norm, PL_slope) = -0.3196 C(Fe_op_FWHM, PL_slope) = -0.3191 C(Fe_op_norm, conti_a_2) = -0.2552 C(Fe_op_shift, PL_slope) = -0.2360 C(Fe_op_FWHM, conti_a_2) = -0.1926 C(Fe_op_norm, Fe_op_shift) = -0.1696 C(Fe_op_norm, conti_a_1) = +0.1279 C(Fe_op_shift, PL_norm) = +0.1221 Name Value Min Max Stderr Vary Expr Brute_Step Hb_br_1_dwave 0 -0.01 0.01 None True None None Hb_br_1_scale 0 0 1e+10 None True None None Hb_br_1_sigma 0.005 0.004 0.05 None True None None Hb_br_2_dwave 0 -0.01 0.01 None True None None Hb_br_2_scale 0 0 1e+10 None True None None Hb_br_2_sigma 0.005 0.004 0.05 None True None None Hb_na_1_dwave 0 -0.01 0.01 None True None None Hb_na_1_scale 0 0 1e+10 None True None None Hb_na_1_sigma 0.001 0.00023 0.00169 None True None None OIII4959c_1_dwave 0 -0.01 0.01 None False Hb_na_1_dwave None OIII4959c_1_scale 0 0 1e+10 None True None None OIII4959c_1_sigma 0.001 0.00023 0.00169 None False Hb_na_1_sigma None OIII4959w_1_dwave 0 -0.01 0.01 None True None None OIII4959w_1_scale 0 0 1e+10 None True None None OIII4959w_1_sigma 0.003 0.00023 0.004 None True None None OIII5007c_1_dwave 0 -0.01 0.01 None False Hb_na_1_dwave None OIII5007c_1_scale 0 0 1e+10 None True None None OIII5007c_1_sigma 0.001 0.00023 0.00169 None False Hb_na_1_sigma None OIII5007w_1_dwave 0 -0.01 0.01 None False OIII4959w_1_dwave None OIII5007w_1_scale 0 0 1e+10 None True None None OIII5007w_1_sigma 0.003 0.00023 0.004 None False OIII4959w_1_sigma None Fitting complex Hb Fit report [[Variables]] Hb_br_1_scale: 0.97885255 +/- 0.36645922 (37.44%) (init = 0) Hb_br_1_dwave: -0.01000000 +/- 0.01463778 (146.38%) (init = 0) Hb_br_1_sigma: 0.03895222 +/- 0.02502504 (64.25%) (init = 0.005) Hb_br_2_scale: 9.28924215 +/- 0.53274340 (5.74%) (init = 0) Hb_br_2_dwave: 4.0417e-04 +/- 3.4322e-04 (84.92%) (init = 0) Hb_br_2_sigma: 0.00838179 +/- 5.7216e-04 (6.83%) (init = 0.005) Hb_na_1_scale: 8.79697926 +/- 0.89487909 (10.17%) (init = 0) Hb_na_1_dwave: -1.3496e-04 +/- 2.0595e-05 (15.26%) (init = 0) Hb_na_1_sigma: 7.4126e-04 +/- 2.3662e-05 (3.19%) (init = 0.001) OIII4959c_1_scale: 26.7812861 +/- 1.52114381 (5.68%) (init = 0) OIII4959c_1_dwave: -1.3496e-04 +/- 2.0595e-05 (15.26%) == 'Hb_na_1_dwave' OIII4959c_1_sigma: 7.4126e-04 +/- 2.3662e-05 (3.19%) == 'Hb_na_1_sigma' OIII5007c_1_scale: 81.9650858 +/- 2.68002209 (3.27%) (init = 0) OIII5007c_1_dwave: -1.3496e-04 +/- 2.0595e-05 (15.26%) == 'Hb_na_1_dwave' OIII5007c_1_sigma: 7.4126e-04 +/- 2.3662e-05 (3.19%) == 'Hb_na_1_sigma' OIII4959w_1_scale: 2.79441081 +/- 1.31726704 (47.14%) (init = 0) OIII4959w_1_dwave: -0.00104684 +/- 3.7752e-04 (36.06%) (init = 0) OIII4959w_1_sigma: 0.00227867 +/- 3.1955e-04 (14.02%) (init = 0.003) OIII5007w_1_scale: 6.13927797 +/- 1.37851630 (22.45%) (init = 0) OIII5007w_1_dwave: -0.00104684 +/- 3.7752e-04 (36.06%) == 'OIII4959w_1_dwave' OIII5007w_1_sigma: 0.00227867 +/- 3.1955e-04 (14.02%) == 'OIII4959w_1_sigma' [[Correlations]] (unreported correlations are < 0.100) C(Hb_br_1_scale, Hb_br_1_sigma) = -0.8865 C(Hb_br_1_scale, Hb_br_1_dwave) = +0.8172 C(Hb_br_1_scale, Hb_br_2_scale) = -0.7835 C(Hb_br_1_scale, Hb_br_2_sigma) = -0.7634 C(Hb_br_1_dwave, Hb_br_1_sigma) = -0.7424 C(Hb_br_1_dwave, Hb_br_2_scale) = -0.7299 C(Hb_na_1_sigma, OIII5007w_1_scale) = -0.7030 C(Hb_br_1_sigma, Hb_br_2_scale) = +0.6883 C(Hb_br_1_sigma, Hb_br_2_sigma) = +0.6469 C(OIII4959c_1_scale, OIII4959w_1_scale) = -0.6248 C(OIII4959w_1_dwave, OIII5007w_1_scale) = +0.6247 C(OIII4959w_1_sigma, OIII5007w_1_scale) = -0.6189 C(Hb_br_1_dwave, Hb_br_2_sigma) = -0.6122 C(OIII4959w_1_scale, OIII5007w_1_scale) = +0.5669 C(Hb_na_1_sigma, OIII4959w_1_dwave) = -0.5428 C(Hb_na_1_sigma, OIII4959w_1_scale) = -0.5305 C(OIII4959w_1_scale, OIII4959w_1_dwave) = +0.4506 C(Hb_na_1_sigma, OIII4959w_1_sigma) = +0.4001 C(OIII4959w_1_scale, OIII4959w_1_sigma) = -0.3987 C(OIII5007c_1_scale, OIII5007w_1_scale) = -0.3832 C(Hb_br_1_dwave, OIII4959w_1_sigma) = -0.3544 C(Hb_na_1_dwave, OIII4959w_1_sigma) = -0.3500 C(Hb_br_2_scale, Hb_br_2_sigma) = +0.3420 C(OIII4959w_1_dwave, OIII4959w_1_sigma) = -0.3267 C(OIII4959c_1_scale, OIII4959w_1_dwave) = -0.3072 C(Hb_br_2_sigma, Hb_na_1_scale) = +0.2948 C(OIII5007c_1_scale, OIII4959w_1_dwave) = -0.2759 C(Hb_br_2_scale, OIII4959w_1_sigma) = +0.2753 C(OIII5007c_1_scale, OIII4959w_1_sigma) = +0.2747 C(OIII4959c_1_scale, OIII4959w_1_sigma) = +0.2711 C(Hb_na_1_dwave, Hb_na_1_sigma) = -0.2698 C(Hb_br_1_scale, OIII4959w_1_sigma) = -0.2683 C(Hb_na_1_dwave, OIII5007w_1_scale) = +0.2652 C(OIII4959c_1_scale, OIII5007w_1_scale) = -0.2543 C(Hb_na_1_dwave, OIII4959w_1_scale) = +0.2171 C(Hb_br_2_scale, Hb_na_1_scale) = -0.2120 C(Hb_br_1_sigma, Hb_br_2_dwave) = -0.1937 C(Hb_br_2_scale, OIII4959w_1_scale) = +0.1766 C(Hb_br_1_scale, Hb_br_2_dwave) = +0.1621 C(OIII4959c_1_scale, OIII5007c_1_scale) = +0.1528 C(Hb_na_1_sigma, OIII4959c_1_scale) = +0.1503 C(Hb_br_2_scale, Hb_br_2_dwave) = -0.1389 C(Hb_br_1_sigma, OIII4959w_1_sigma) = +0.1388 C(Hb_br_2_sigma, OIII4959w_1_sigma) = +0.1314 C(Hb_br_1_dwave, OIII4959w_1_scale) = -0.1302 C(OIII5007c_1_scale, OIII4959w_1_scale) = -0.1144 C(Hb_br_2_dwave, OIII4959w_1_dwave) = +0.1128 C(Hb_na_1_dwave, OIII5007c_1_scale) = -0.1113 C(Hb_br_1_scale, Hb_na_1_scale) = -0.1094 C(Hb_br_2_sigma, OIII5007w_1_scale) = +0.1070 C(Hb_na_1_scale, Hb_na_1_sigma) = -0.1010 Name Value Min Max Stderr Vary Expr Brute_Step Ha_br_1_dwave 0 -0.015 0.015 None True None None Ha_br_1_scale 0 0 1e+10 None True None None Ha_br_1_sigma 0.005 0.004 0.05 None True None None Ha_br_2_dwave 0 -0.015 0.015 None True None None Ha_br_2_scale 0 0 1e+10 None True None None Ha_br_2_sigma 0.005 0.004 0.05 None True None None Ha_na_1_dwave 0 -0.01 0.01 None True None None Ha_na_1_scale 0 0 1e+10 None True None None Ha_na_1_sigma 0.001 0.0005 0.00169 None True None None NII6549_1_dwave 0 -0.005 0.005 None False Ha_na_1_dwave None NII6549_1_scale 0 0 1e+10 None True None None NII6549_1_sigma 0.001 0.00023 0.00169 None False Ha_na_1_sigma None NII6585_1_dwave 0 -0.005 0.005 None False Ha_na_1_dwave None NII6585_1_scale 0 0 1e+10 None False 3.0 * NII6549_1_scale None NII6585_1_sigma 0.001 0.00023 0.00169 None False Ha_na_1_sigma None SII6718_1_dwave 0 -0.005 0.005 None False Ha_na_1_dwave None SII6718_1_scale 0 0 1e+10 None True None None SII6718_1_sigma 0.001 0.00023 0.00169 None False Ha_na_1_sigma None SII6732_1_dwave 0 -0.005 0.005 None False Ha_na_1_dwave None SII6732_1_scale 0 0 1e+10 None False 1.0 * SII6718_1_scale None SII6732_1_sigma 0.001 0.00023 0.00169 None False Ha_na_1_sigma None Fitting complex Ha Fit report [[Variables]] Ha_br_1_scale: 17.6485809 +/- 1.79398724 (10.17%) (init = 0) Ha_br_1_dwave: 0.00139980 +/- 1.9594e-04 (14.00%) (init = 0) Ha_br_1_sigma: 0.00499016 +/- 3.8904e-04 (7.80%) (init = 0.005) Ha_br_2_scale: 17.5080589 +/- 1.92454297 (10.99%) (init = 0) Ha_br_2_dwave: -7.5188e-04 +/- 3.1032e-04 (41.27%) (init = 0) Ha_br_2_sigma: 0.01001453 +/- 3.8666e-04 (3.86%) (init = 0.005) Ha_na_1_scale: 21.0241530 +/- 1.33829912 (6.37%) (init = 0) Ha_na_1_dwave: -2.3870e-04 +/- 2.4802e-05 (10.39%) (init = 0) Ha_na_1_sigma: 6.8370e-04 +/- 2.6017e-05 (3.81%) (init = 0.001) NII6549_1_scale: 9.36066780 +/- 0.44070666 (4.71%) (init = 0) NII6549_1_dwave: -2.3870e-04 +/- 2.4802e-05 (10.39%) == 'Ha_na_1_dwave' NII6549_1_sigma: 6.8370e-04 +/- 2.6017e-05 (3.81%) == 'Ha_na_1_sigma' NII6585_1_scale: 28.0820034 +/- 1.32212000 (4.71%) == '3.0 * NII6549_1_scale' NII6585_1_dwave: -2.3870e-04 +/- 2.4802e-05 (10.39%) == 'Ha_na_1_dwave' NII6585_1_sigma: 6.8370e-04 +/- 2.6017e-05 (3.81%) == 'Ha_na_1_sigma' SII6718_1_scale: 11.1556758 +/- 0.47862446 (4.29%) (init = 0) SII6718_1_dwave: -2.3870e-04 +/- 2.4802e-05 (10.39%) == 'Ha_na_1_dwave' SII6718_1_sigma: 6.8370e-04 +/- 2.6017e-05 (3.81%) == 'Ha_na_1_sigma' SII6732_1_scale: 11.1556758 +/- 0.47862446 (4.29%) == '1.0 * SII6718_1_scale' SII6732_1_dwave: -2.3870e-04 +/- 2.4802e-05 (10.39%) == 'Ha_na_1_dwave' SII6732_1_sigma: 6.8370e-04 +/- 2.6017e-05 (3.81%) == 'Ha_na_1_sigma' [[Correlations]] (unreported correlations are < 0.100) C(Ha_br_2_scale, Ha_br_2_sigma) = -0.9136 C(Ha_br_1_sigma, Ha_br_2_scale) = -0.8752 C(Ha_br_1_scale, Ha_br_2_scale) = -0.8739 C(Ha_br_1_scale, Ha_br_2_sigma) = +0.8645 C(Ha_br_2_scale, Ha_br_2_dwave) = +0.7274 C(Ha_br_1_sigma, Ha_br_2_dwave) = -0.7228 C(Ha_br_1_sigma, Ha_br_2_sigma) = +0.6808 C(Ha_br_1_scale, Ha_br_2_dwave) = -0.5955 C(Ha_br_1_scale, Ha_br_1_sigma) = +0.5753 C(Ha_br_2_dwave, Ha_br_2_sigma) = -0.5481 C(Ha_na_1_scale, NII6549_1_scale) = +0.5196 C(Ha_na_1_sigma, SII6718_1_scale) = -0.5106 C(Ha_br_1_dwave, Ha_br_2_sigma) = -0.4897 C(Ha_br_1_sigma, Ha_na_1_scale) = +0.4448 C(Ha_br_1_scale, Ha_br_1_dwave) = -0.4298 C(Ha_br_1_sigma, NII6549_1_scale) = +0.4135 C(Ha_br_1_sigma, Ha_na_1_sigma) = +0.3909 C(Ha_br_1_dwave, Ha_br_2_scale) = +0.3463 C(Ha_br_2_dwave, Ha_na_1_scale) = -0.2855 C(Ha_na_1_scale, Ha_na_1_sigma) = +0.2452 C(Ha_br_2_dwave, Ha_na_1_sigma) = -0.2427 C(Ha_br_1_dwave, Ha_na_1_scale) = +0.2383 C(Ha_br_2_dwave, NII6549_1_scale) = -0.2187 C(Ha_br_1_sigma, SII6718_1_scale) = -0.2177 C(Ha_br_2_scale, Ha_na_1_scale) = -0.2177 C(Ha_br_2_scale, NII6549_1_scale) = -0.2173 C(Ha_br_1_dwave, Ha_br_2_dwave) = -0.2146 C(Ha_br_2_sigma, SII6718_1_scale) = -0.2005 C(Ha_br_1_dwave, Ha_na_1_dwave) = -0.1913 C(Ha_br_2_scale, SII6718_1_scale) = +0.1763 C(Ha_br_2_scale, Ha_na_1_sigma) = -0.1758 C(Ha_na_1_sigma, NII6549_1_scale) = +0.1757 C(Ha_br_1_dwave, Ha_na_1_sigma) = +0.1723 C(Ha_br_1_scale, Ha_na_1_sigma) = -0.1471 C(Ha_br_1_dwave, Ha_br_1_sigma) = -0.1337 C(Ha_br_2_sigma, NII6549_1_scale) = +0.1256 C(Ha_br_1_scale, Ha_na_1_scale) = -0.1203 C(Ha_br_1_scale, NII6549_1_scale) = -0.1115 Saving figure as ./0332-52367-0639.pdf Fitting finished in 39.0s
Step 3(optional): Compute parameter unertaintines using MCMC:
For comparison, let's fit the same quasar spectrum, but this time using MCMC to get uncertainties on the parameters. Obviously, this will take significantly longer depending on the number of samples you want to run. You will need the pandas
and emcee
packages to do this. To visualize the posterior parameter distributions and their covariances, you will need the corner
package. It's also usually a good idea to perturb the parameter starting points by a small random number using the epsilon_jitter
argument, and exclude unnecessary line components, otherwise you may get an ''Initial state has a large condition number'' error.
q_mcmc = QSOFit(lam, flux, err, z, ra=ra, dec=dec, plateid=plateid, mjd=mjd, fiberid=fiberid, path=path_ex)
start = timeit.default_timer()
# Do the fitting
q_mcmc.Fit(name=None, nsmooth=1, deredden=True, reject_badpix=False, wave_range=None, \
wave_mask=None, decompose_host=True, host_prior=True, decomp_na_mask=True, npca_gal=5, npca_qso=10, qso_type='CZBIN1',
Fe_uv_op=True, poly=True, rej_abs_conti=False, \
MCMC=True, epsilon_jitter=0, nburn=10, nsamp=400, nthin=10, linefit=True, save_result=True, \
plot_fig=True, save_fig=False, plot_corner=True, kwargs_plot={'save_fig_path': '.'}, save_fits_name=None,
verbose=False)
end = timeit.default_timer()
print(f'Fitting finished in {np.round(end - start, 1)}s')
100%|██████████| 400/400 [00:17<00:00, 22.31it/s]
The chain is shorter than 50 times the integrated autocorrelation time for 8 parameter(s). Use this estimate with caution and run a longer chain! N/50 = 8; tau: [37.4819265 35.51595007 41.7142146 34.94352008 38.95248101 37.13491757 36.19169448 35.51427 ]
WARNING:root:Too few points to create valid contours 100%|██████████| 400/400 [00:02<00:00, 182.23it/s]
The chain is shorter than 50 times the integrated autocorrelation time for 15 parameter(s). Use this estimate with caution and run a longer chain! N/50 = 8; tau: [33.70053779 41.26235867 33.7287362 37.86308569 44.7735224 35.90809577 34.13139874 34.64956895 33.99660018 33.54001231 37.00544666 36.96485327 37.58347963 35.39112833 40.56997327]
WARNING:root:Too few points to create valid contours WARNING:root:Too few points to create valid contours 100%|██████████| 400/400 [00:02<00:00, 135.43it/s]
The chain is shorter than 50 times the integrated autocorrelation time for 11 parameter(s). Use this estimate with caution and run a longer chain! N/50 = 8; tau: [33.43639297 32.96977501 35.09388324 33.73922487 32.69749807 31.74444633 34.18940149 31.55221483 31.92628651 31.61301339 32.19395998] Fitting finished in 50.6s
When MCMC=True
or MC=True
(they cannot both be True
), the 1$\sigma$ uncertainties (from the 16th and 84th percentiles) are added to the results (the "_err" parameters). The initial MLE fit will again be computed and is used as the starting point for emcee.
You may need to turn off/on continuum fitting components or adjust the parameter bounds as neccessary to fit your data. From the above plots, it is clear that the parameter posteriors are reasonably well-constrained. For more details, such as to view the acceptance fraction, use verbose=True
.
View results: The results arrays will include the uncertainties on each parameter only if MCMC is used.
# Continuum fitting results
print(q_mcmc.conti_result_name)
print('')
print(q_mcmc.conti_result)
['ra' 'dec' 'plateid' 'MJD' 'fiberid' 'redshift' 'SN_ratio_conti' 'Fe_uv_norm' 'Fe_uv_norm_err' 'Fe_uv_FWHM' 'Fe_uv_FWHM_err' 'Fe_uv_shift' 'Fe_uv_shift_err' 'Fe_op_norm' 'Fe_op_norm_err' 'Fe_op_FWHM' 'Fe_op_FWHM_err' 'Fe_op_shift' 'Fe_op_shift_err' 'PL_norm' 'PL_norm_err' 'PL_slope' 'PL_slope_err' 'Blamer_norm' 'Blamer_norm_err' 'Balmer_Te' 'Balmer_Te_err' 'Balmer_Tau' 'Balmer_Tau_err' 'conti_a_0' 'conti_a_0_err' 'conti_a_1' 'conti_a_1_err' 'conti_a_2' 'conti_a_2_err' 'L1350' 'L1350_err' 'L1450' 'L1450_err' 'L3000' 'L3000_err' 'L4200' 'L4200_err' 'L5100' 'L5100_err' 'Fe_flux_4435_4685' 'Fe_flux_4435_4685_err' 'SN_host' 'rchi2_decomp' 'frac_host_4200' 'frac_host_5100' 'Dn4000' 'sigma' 'sigma_err' 'v_off' 'v_off_err' 'rchi2_ppxf' 'gal_par_0' 'gal_par_1' 'gal_par_2' 'gal_par_3' 'gal_par_4' 'qso_par_0' 'qso_par_1' 'qso_par_2' 'qso_par_3' 'qso_par_4' 'qso_par_5' 'qso_par_6' 'qso_par_7' 'qso_par_8' 'qso_par_9'] ['184.03066' '-2.2382793' '332' '52367' '639' '0.100609764' '27.49352773461182' '0.0' '0.0' '3000.0' '0.0' '0.0' '0.0' '3.4653768743453384' '0.4520118294501865' '16514.26873660736' '1994.2299200536136' '0.00999970337371383' '0.0006495549071946693' '35.98332409815441' '1.5073727190304567' '-2.1765103796164773' '0.22315560538563628' '0.0' '0.0' '15000.0' '0.0' '0.5' '0.0' '-1418.2680823659036' '1064.6747402732592' '1.540241670201466' '0.5204545035260251' '-0.00020214517831770686' '7.26597653010191e-05' '-1.0' '-1.0' '-1.0' '-1.0' '-1.0' '-1.0' '43.275094703694236' '0.013776195387460888' '43.24038929224117' '0.013259365194382866' '1189.579667283605' '136.88153778397674' '21.2176479440606' '7.193250913730492' '0.5976026679302436' '0.6896394465041954' '1.3852998823187828' '182.47746578157992' '14.066053270793734' '22.43197328418017' '12.233758125545412' '0.8298778955844628' '2246.1465595808727' '-134.09707711409297' '-224.52256594147457' '77.92240725419694' '126.81605381012216' '1141.1640555062716' '103.63580924406342' '40.9744563618422' '4.358677949595329' '13.622487431287192' '-16.983273438966318' '-8.250254529165726' '6.718401087449956' '-3.9941820690989633' '-1.2141964427794165']
# Gaussian fitting results
print(q_mcmc.gauss_result_name)
print('')
print(q_mcmc.gauss_result)
['Hb_br_1_scale' 'Hb_br_1_scale_err' 'Hb_br_1_centerwave' 'Hb_br_1_centerwave_err' 'Hb_br_1_sigma' 'Hb_br_1_sigma_err' 'Hb_br_2_scale' 'Hb_br_2_scale_err' 'Hb_br_2_centerwave' 'Hb_br_2_centerwave_err' 'Hb_br_2_sigma' 'Hb_br_2_sigma_err' 'Hb_na_1_scale' 'Hb_na_1_scale_err' 'Hb_na_1_centerwave' 'Hb_na_1_centerwave_err' 'Hb_na_1_sigma' 'Hb_na_1_sigma_err' 'OIII4959c_1_scale' 'OIII4959c_1_scale_err' 'OIII4959c_1_centerwave' 'OIII4959c_1_centerwave_err' 'OIII4959c_1_sigma' 'OIII4959c_1_sigma_err' 'OIII5007c_1_scale' 'OIII5007c_1_scale_err' 'OIII5007c_1_centerwave' 'OIII5007c_1_centerwave_err' 'OIII5007c_1_sigma' 'OIII5007c_1_sigma_err' 'OIII4959w_1_scale' 'OIII4959w_1_scale_err' 'OIII4959w_1_centerwave' 'OIII4959w_1_centerwave_err' 'OIII4959w_1_sigma' 'OIII4959w_1_sigma_err' 'OIII5007w_1_scale' 'OIII5007w_1_scale_err' 'OIII5007w_1_centerwave' 'OIII5007w_1_centerwave_err' 'OIII5007w_1_sigma' 'OIII5007w_1_sigma_err' 'Ha_br_1_scale' 'Ha_br_1_scale_err' 'Ha_br_1_centerwave' 'Ha_br_1_centerwave_err' 'Ha_br_1_sigma' 'Ha_br_1_sigma_err' 'Ha_br_2_scale' 'Ha_br_2_scale_err' 'Ha_br_2_centerwave' 'Ha_br_2_centerwave_err' 'Ha_br_2_sigma' 'Ha_br_2_sigma_err' 'Ha_na_1_scale' 'Ha_na_1_scale_err' 'Ha_na_1_centerwave' 'Ha_na_1_centerwave_err' 'Ha_na_1_sigma' 'Ha_na_1_sigma_err' 'NII6549_1_scale' 'NII6549_1_scale_err' 'NII6549_1_centerwave' 'NII6549_1_centerwave_err' 'NII6549_1_sigma' 'NII6549_1_sigma_err' 'NII6585_1_scale' 'NII6585_1_scale_err' 'NII6585_1_centerwave' 'NII6585_1_centerwave_err' 'NII6585_1_sigma' 'NII6585_1_sigma_err' 'SII6718_1_scale' 'SII6718_1_scale_err' 'SII6718_1_centerwave' 'SII6718_1_centerwave_err' 'SII6718_1_sigma' 'SII6718_1_sigma_err' 'SII6732_1_scale' 'SII6732_1_scale_err' 'SII6732_1_centerwave' 'SII6732_1_centerwave_err' 'SII6732_1_sigma' 'SII6732_1_sigma_err'] [9.24570753e+00 4.57458560e-01 8.48977663e+00 1.50289488e-04 8.38733724e-03 4.61882186e-04 1.00050968e+00 2.20800023e-01 8.47934460e+00 2.90558737e-04 3.56431243e-02 6.12564319e-03 8.70221728e+00 9.47224151e-01 8.48916819e+00 1.61297655e-05 7.68945372e-04 1.79167482e-05 2.58883515e+01 1.25915969e+00 8.50904467e+00 0.00000000e+00 7.68945372e-04 0.00000000e+00 7.57251267e+01 2.07416766e+00 8.51866343e+00 0.00000000e+00 7.68945372e-04 0.00000000e+00 2.42857234e+00 4.94428636e-01 8.50831894e+00 2.96583598e-04 2.52713678e-03 3.02580673e-04 4.93992180e+00 8.12884070e-01 8.51793770e+00 0.00000000e+00 2.52713678e-03 0.00000000e+00 1.75345372e+01 2.20059849e+00 8.79090418e+00 2.50947656e-04 5.17661357e-03 3.56777786e-04 1.70226511e+01 2.27295287e+00 8.78855166e+00 2.37341465e-04 1.00053187e-02 5.01190245e-04 2.19960189e+01 1.35685644e+00 8.78922856e+00 2.19726550e-05 6.63311069e-04 2.21560669e-05 9.62130375e+00 4.39754413e-01 8.78697693e+00 0.00000000e+00 6.63311069e-04 0.00000000e+00 2.88639113e+01 0.00000000e+00 8.79237187e+00 0.00000000e+00 6.63311069e-04 0.00000000e+00 1.14277110e+01 5.39845869e-01 8.81236851e+00 0.00000000e+00 6.63311069e-04 0.00000000e+00 1.14277110e+01 0.00000000e+00 8.81450665e+00 0.00000000e+00 6.63311069e-04 0.00000000e+00]
To get the results of each MC or MCMC sample, use gauss_result_all
. You can then use the line_prop
function to get the line properties for each sample, and take median values of each parameter's samples. The parameters in gauss_result
are used for the MLE solution for MCMC (or just one MC sample), and may not be the same as the median of the posterior distributions!
# Gaussian fitting results
print(q_mcmc.gauss_result_name[::2])
print('')
print(q_mcmc.gauss_result_all)
print(np.shape(q_mcmc.gauss_result_all))
['Hb_br_1_scale' 'Hb_br_1_centerwave' 'Hb_br_1_sigma' 'Hb_br_2_scale' 'Hb_br_2_centerwave' 'Hb_br_2_sigma' 'Hb_na_1_scale' 'Hb_na_1_centerwave' 'Hb_na_1_sigma' 'OIII4959c_1_scale' 'OIII4959c_1_centerwave' 'OIII4959c_1_sigma' 'OIII5007c_1_scale' 'OIII5007c_1_centerwave' 'OIII5007c_1_sigma' 'OIII4959w_1_scale' 'OIII4959w_1_centerwave' 'OIII4959w_1_sigma' 'OIII5007w_1_scale' 'OIII5007w_1_centerwave' 'OIII5007w_1_sigma' 'Ha_br_1_scale' 'Ha_br_1_centerwave' 'Ha_br_1_sigma' 'Ha_br_2_scale' 'Ha_br_2_centerwave' 'Ha_br_2_sigma' 'Ha_na_1_scale' 'Ha_na_1_centerwave' 'Ha_na_1_sigma' 'NII6549_1_scale' 'NII6549_1_centerwave' 'NII6549_1_sigma' 'NII6585_1_scale' 'NII6585_1_centerwave' 'NII6585_1_sigma' 'SII6718_1_scale' 'SII6718_1_centerwave' 'SII6718_1_sigma' 'SII6732_1_scale' 'SII6732_1_centerwave' 'SII6732_1_sigma'] [[9.23674880e+00 8.48977688e+00 8.36799256e-03 ... 1.14277110e+01 8.81450665e+00 6.63311069e-04] [9.24108953e+00 8.48977529e+00 8.36693646e-03 ... 1.14277110e+01 8.81450665e+00 6.63311069e-04] [9.28367498e+00 8.48977365e+00 8.34374060e-03 ... 1.14277110e+01 8.81450665e+00 6.63311069e-04] ... [9.17438998e+00 8.49018783e+00 9.37855183e-03 ... 1.14277110e+01 8.81450665e+00 6.63311069e-04] [9.65746530e+00 8.48946934e+00 8.20387795e-03 ... 1.14277110e+01 8.81450665e+00 6.63311069e-04] [9.19351983e+00 8.49020467e+00 8.84218664e-03 ... 1.14277110e+01 8.81450665e+00 6.63311069e-04]] (3900, 42)
# line fitting results
print(q_mcmc.line_result_name)
print('')
print(q_mcmc.line_result)
['1_complex_name' '1_line_status' '1_line_min_chi2' '1_line_bic' '1_line_red_chi2' '1_niter' '1_ndof' '2_complex_name' '2_line_status' '2_line_min_chi2' '2_line_bic' '2_line_red_chi2' '2_niter' '2_ndof' 'Hb_br_1_scale' 'Hb_br_1_scale_err' 'Hb_br_1_centerwave' 'Hb_br_1_centerwave_err' 'Hb_br_1_sigma' 'Hb_br_1_sigma_err' 'Hb_br_2_scale' 'Hb_br_2_scale_err' 'Hb_br_2_centerwave' 'Hb_br_2_centerwave_err' 'Hb_br_2_sigma' 'Hb_br_2_sigma_err' 'Hb_na_1_scale' 'Hb_na_1_scale_err' 'Hb_na_1_centerwave' 'Hb_na_1_centerwave_err' 'Hb_na_1_sigma' 'Hb_na_1_sigma_err' 'OIII4959c_1_scale' 'OIII4959c_1_scale_err' 'OIII4959c_1_centerwave' 'OIII4959c_1_centerwave_err' 'OIII4959c_1_sigma' 'OIII4959c_1_sigma_err' 'OIII5007c_1_scale' 'OIII5007c_1_scale_err' 'OIII5007c_1_centerwave' 'OIII5007c_1_centerwave_err' 'OIII5007c_1_sigma' 'OIII5007c_1_sigma_err' 'OIII4959w_1_scale' 'OIII4959w_1_scale_err' 'OIII4959w_1_centerwave' 'OIII4959w_1_centerwave_err' 'OIII4959w_1_sigma' 'OIII4959w_1_sigma_err' 'OIII5007w_1_scale' 'OIII5007w_1_scale_err' 'OIII5007w_1_centerwave' 'OIII5007w_1_centerwave_err' 'OIII5007w_1_sigma' 'OIII5007w_1_sigma_err' 'Ha_br_1_scale' 'Ha_br_1_scale_err' 'Ha_br_1_centerwave' 'Ha_br_1_centerwave_err' 'Ha_br_1_sigma' 'Ha_br_1_sigma_err' 'Ha_br_2_scale' 'Ha_br_2_scale_err' 'Ha_br_2_centerwave' 'Ha_br_2_centerwave_err' 'Ha_br_2_sigma' 'Ha_br_2_sigma_err' 'Ha_na_1_scale' 'Ha_na_1_scale_err' 'Ha_na_1_centerwave' 'Ha_na_1_centerwave_err' 'Ha_na_1_sigma' 'Ha_na_1_sigma_err' 'NII6549_1_scale' 'NII6549_1_scale_err' 'NII6549_1_centerwave' 'NII6549_1_centerwave_err' 'NII6549_1_sigma' 'NII6549_1_sigma_err' 'NII6585_1_scale' 'NII6585_1_scale_err' 'NII6585_1_centerwave' 'NII6585_1_centerwave_err' 'NII6585_1_sigma' 'NII6585_1_sigma_err' 'SII6718_1_scale' 'SII6718_1_scale_err' 'SII6718_1_centerwave' 'SII6718_1_centerwave_err' 'SII6718_1_sigma' 'SII6718_1_sigma_err' 'SII6732_1_scale' 'SII6732_1_scale_err' 'SII6732_1_centerwave' 'SII6732_1_centerwave_err' 'SII6732_1_sigma' 'SII6732_1_sigma_err' 'Hb_whole_br_fwhm' 'Hb_whole_br_fwhm_err' 'Hb_whole_br_sigma' 'Hb_whole_br_sigma_err' 'Hb_whole_br_ew' 'Hb_whole_br_ew_err' 'Hb_whole_br_peak' 'Hb_whole_br_peak_err' 'Hb_whole_br_area' 'Hb_whole_br_area_err' 'Hb_whole_br_snr' 'Hb_whole_br_snr_err' 'Ha_whole_br_fwhm' 'Ha_whole_br_fwhm_err' 'Ha_whole_br_sigma' 'Ha_whole_br_sigma_err' 'Ha_whole_br_ew' 'Ha_whole_br_ew_err' 'Ha_whole_br_peak' 'Ha_whole_br_peak_err' 'Ha_whole_br_area' 'Ha_whole_br_area_err' 'Ha_whole_br_snr' 'Ha_whole_br_snr_err'] ['Hb' '1' '580.1602343338035' '231.9551672972695' '1.465051096802534' '481' '396' 'Ha' '1' '272.51355862655873' '70.63922501532213' '1.081403010422852' '461' '252' '9.245707532912206' '0.4574585599029657' '8.489776632417794' '0.0001502894882303707' '0.008387337243509492' '0.0004618821864144085' '1.000509675108674' '0.22080002327528808' '8.47934459708632' '0.0002905587368084994' '0.035643124331113905' '0.006125643189467489' '8.702217280109892' '0.9472241507338639' '8.48916819002536' '1.612976546283562e-05' '0.0007689453715191508' '1.791674822440608e-05' '25.888351462555192' '1.2591596893397945' '8.5090446701279' '0.0' '0.0007689453715191508' '0.0' '75.72512672471277' '2.0741676615408196' '8.518663429283173' '0.0' '0.0007689453715191508' '0.0' '2.428572343227131' '0.4944286364805236' '8.508318938499235' '0.00029658359843143245' '0.0025271367815153526' '0.00030258067341768034' '4.939921804947289' '0.8128840701525444' '8.517937697654508' '0.0' '0.0025271367815153526' '0.0' '17.534537155583507' '2.2005984909529825' '8.790904178291377' '0.00025094765638389305' '0.005176613567212474' '0.0003567777862591684' '17.022651066511685' '2.272952866908188' '8.788551660425895' '0.000237341465169294' '0.010005318680656992' '0.0005011902445280717' '21.99601889518732' '1.3568564398533596' '8.789228559559707' '2.197265501102663e-05' '0.0006633110689751761' '2.2156066899894542e-05' '9.621303753704069' '0.4397544126974706' '8.786976934498671' '0.0' '0.0006633110689751761' '0.0' '28.863911261112207' '0.0' '8.792371870106582' '0.0' '0.0006633110689751761' '0.0' '11.427710999711849' '0.5398458685605938' '8.812368513172988' '0.0' '0.0006633110689751761' '0.0' '11.427710999711849' '0.0' '8.81450665099037' '0.0' '0.0006633110689751761' '0.0' '12076.642928570596' '2993.8530182105924' '6574.7937686037485' '910.1084351825434' '115.20336557626436' '5.6627590863522315' '4863.972353331239' '0.8052260005688368' '1607.9258117180798' '79.45705467012488' '5.26605024034905' '0.2244944841264478' '4957.14570886338' '127.36803464952754' '2576.7962297231816' '60.27191221214139' '360.253467871212' '4.8517514487768665' '6571.660502835031' '1.0506469564152212' '4288.446334712329' '57.71165004842942' '26.648048104161944' '0.9988587932366055']
Step 4: Get all models for the whole spectrum
Continue to look at this section and below if you want to do some further calculations based on the fitting results. Here, we show how to extract different models from our fitting results, such as continuum model, emission line models and host galaxy component. Note that the emission regions of host galaxy template should be blocked, e.g., H$\alpha$ [6540, 6590].
There are two ways to calculate the Fe II flux within given ranges:
Fe_flux_range
in Fit()
to assign ranges.Get_Fe_flux
to calculate the Fe II flux directly after you finished the fitting part.fig, ax = plt.subplots(figsize=(15, 5))
# Plot the quasar rest frame spectrum after removed the host galaxy component
ax.plot(q_mcmc.wave, q_mcmc.flux, 'grey', label='Data')
ax.plot(q_mcmc.wave, q_mcmc.err, 'r', label='Error')
# Skip the error results before plotting
if q_mcmc.MCMC == True:
gauss_result = q_mcmc.gauss_result[::2]
else:
gauss_result = q_mcmc.gauss_result
# To plot the whole model, we use Manygauss to show the line fitting results saved in gauss_result
ax.plot(q_mcmc.wave, q_mcmc.Manygauss(np.log(q_mcmc.wave), gauss_result) + q_mcmc.f_conti_model, 'b', label='Line',
lw=2)
ax.plot(q_mcmc.wave, q_mcmc.f_conti_model, 'c', lw=2, label='Continuum+FeII')
ax.plot(q_mcmc.wave, q_mcmc.PL_poly_BC, 'orange', lw=2, label='Continuum')
ax.plot(q_mcmc.wave, q_mcmc.host, 'm', lw=2, label='Host')
plt.legend()
ax.set_xlim(3500, 8000)
ax.set_xlabel(r'$\rm Rest \, Wavelength$ ($\rm \AA$)', fontsize=20)
ax.set_ylabel(r'$\rm f_{\lambda}$ ($\rm 10^{-17} erg\;s^{-1}\;cm^{-2}\;\AA^{-1}$)', fontsize=20)
#print('optical Fe flux (10^(-17) erg/s/cm^2): ' + q_mcmc.conti_result[q_mcmc.conti_result_name=='Fe_flux_4435_4685'][0])
Fe_flux_result, Fe_flux_type, Fe_flux_name = q_mcmc.Get_Fe_flux(np.array([4400, 4900]))
print('Fe flux within a specific range: \n' + Fe_flux_name[0] + ': ' + str(Fe_flux_result[0]))
Fe flux within a specific range: Fe_flux_4400_4900: 1977.1658429339636
Step 5: Get models for each line complex
All the line parameters are in the gauss_result, it was saved following the order of qsopar.fits
. For each Gaussian, three parameter (scale, ln lambda, ln sigma) are saved.
If you want to filter the line below or above 1200 km/s (ln sigma = 0.00169), the function CalFWHM
can be used as following.
The line_prop is designed to calculate the broad line properties, not for the narrow line. It is saved in the line_result property, but can be computed again.
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
if q_mcmc.MCMC:
gauss_result = q_mcmc.gauss_result[::2]
else:
gauss_result = q_mcmc.gauss_result
# Plot individual line components
for p in range(int(len(gauss_result) / 3)):
if q_mcmc.CalFWHM(gauss_result[3 * p + 2]) < 1200: # < 1200 km/s narrow
color = 'g' # narrow
else:
color = 'r' # broad
ax.plot(q_mcmc.wave, q_mcmc.Onegauss(np.log(q_mcmc.wave), gauss_result[p * 3:(p + 1) * 3]), color=color)
# Plot total line model
ax.plot(q_mcmc.wave, q_mcmc.Manygauss(np.log(q_mcmc.wave), gauss_result), 'b', lw=2)
ax.plot(q_mcmc.wave, q_mcmc.line_flux, 'k')
ax.set_xlim(4640, 5100)
ax.set_xlabel(r'$\rm Rest \, Wavelength$ ($\rm \AA$)', fontsize=20)
ax.set_ylabel(r'$\rm f_{\lambda}$ ($\rm 10^{-17} erg\;s^{-1}\;cm^{-2}\;\AA^{-1}$)', fontsize=20)
"""
Line properties
"""
# The line_prop function is used to calculate the broad line properties
# (defined, by default, as ln sigma > 0.00169 (1200 km/s) )
# OLD WAY: If you want to calculate the paramters of broad Hb
# then find all the broad Hb component, i.e., Hb_br_[1,2,3]_[scale,centerwave,sigma]
# for here q.line_result_name[12:15], q.line_result[12:15] is the broad Hb
# If MCMC=False, this would be:
# fwhm, sigma, ew, peak, area = q.line_prop(q.linelist[6][0], q.line_result[12:15], 'broad')
# NEW WAY: using line_prop_from_name convenience function
fwhm, sigma, ew, peak, area, snr = q_mcmc.line_prop_from_name('Hb_br', 'broad')
print("Broad Hb:")
print("FWHM (km/s)", np.round(fwhm, 1))
print("Sigma (km/s)", np.round(sigma, 1))
print("EW (A)", np.round(ew, 1))
print("Peak (A)", np.round(peak, 1))
print("Area (10^(-17) erg/s/cm^2)", np.round(area, 1))
print("")
# OLD WAY: If you want to calculate the the narrow [OIII]5007
# If MCMC=False, this would be:
# the coresponding parameters are q.line_result_name[21:24], q.line_result[21:24]
# fwhm, sigma, ew, peak, area = q.line_prop(q.linelist[6][0], q.line_result[21:24], 'narrow')
fwhm, sigma, ew, peak, area, snr = q_mcmc.line_prop_from_name('OIII5007c', 'narrow')
print("Narrow [OIII]5007:")
print("FWHM (km/s)", np.round(fwhm, 1))
print("Sigma (km/s)", np.round(sigma, 1))
print("EW (A)", np.round(ew, 1))
print("Peak (A)", np.round(peak, 1))
print("Area (10^(-17) erg/s/cm^2)", np.round(area, 1))
Broad Hb: FWHM (km/s) 6338.7 Sigma (km/s) 6320.6 EW (A) 98.1 Peak (A) 4864.0 Area (10^(-17) erg/s/cm^2) 1375.1 Narrow [OIII]5007: FWHM (km/s) 547.6 Sigma (km/s) 227.2 EW (A) 53.9 Peak (A) 5006.7 Area (10^(-17) erg/s/cm^2) 728.7
Determine number of broad components for a line:
If you want to determine how many Gaussian components are statistically justified for a broad emission line given the S/N, run PyQSOFit
multiple times and use the Bayesian information criterion to determine it.
# Required
data = fits.open(os.path.join(path_ex, 'data/spec-0332-52367-0639.fits'))
lam = 10 ** data[1].data['loglam'] # OBS wavelength [A]
flux = data[1].data['flux'] # OBS flux [erg/s/cm^2/A]
err = 1 / np.sqrt(data[1].data['ivar']) # 1 sigma error
z = data[2].data['z'][0] # Redshift
# Optional
ra = data[0].header['plug_ra'] # RA
dec = data[0].header['plug_dec'] # DEC
plateid = data[0].header['plateid'] # SDSS plate ID
mjd = data[0].header['mjd'] # SDSS MJD
fiberid = data[0].header['fiberid'] # SDSS fiber ID
ngauss_max = 5 # stop at 5 components
bic_last = np.inf
# Number of Gaussians to try loop
for ngauss in range(1, ngauss_max):
print(fr'Fitting broad H$\alpha$ with {ngauss} components.')
# Change the number of Gaussians for the Ha line in the line parameter file
hdu1_new = hdu1.copy()
if 'Ha_br' in hdu1.data['linename']:
hdu1.data['ngauss'][hdu1.data['linename'] == 'Ha_br'] = ngauss
hdu_list = fits.HDUList([primary_hdu, hdu1, hdu2, hdu3, hdu4])
hdu_list.writeto(os.path.join(path_ex, 'qsopar.fits'), overwrite=True)
# Do the fitting
q = QSOFit(lam, flux, err, z, ra=ra, dec=dec, plateid=plateid, mjd=mjd, fiberid=fiberid, path=path_ex)
q.Fit(name=None, nsmooth=1, deredden=True, reject_badpix=False, wave_range=None, \
wave_mask=None, decompose_host=True, host_prior=True, decomp_na_mask=True, npca_gal=5, npca_qso=10, qso_type='CZBIN1',\
Fe_uv_op=True, poly=True, BC=False, rej_abs_conti=False, rej_abs_line=False, \
initial_guess=None, MCMC=False, nburn=20, nsamp=200, nthin=10, linefit=True, \
save_result=False, plot_fig=False, verbose=False)
mask_Ha_bic = q.line_result_name == '2_line_min_chi2'
bic = float(q.line_result[mask_Ha_bic][0])
print(f'Delta BIC = {np.round(bic_last - bic, 1)}')
# Stop condition of Delta BIC = 10 is a good rule of thumb
if bic_last - bic < 10:
print(f'{ngauss-1} components is prefered')
break
bic_last = bic
# Plot the result
if q.MCMC:
gauss_result = q.gauss_result[::2]
else:
gauss_result = q.gauss_result
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
# Plot individual line components
for p in range(len(gauss_result) // 3):
if q.CalFWHM(gauss_result[3 * p + 2]) < 1200: # < 1200 km/s narrow
color = 'g' # narrow
else:
color = 'r' # broad
ax.plot(q.wave, q.Onegauss(np.log(q.wave), gauss_result[p * 3:(p + 1) * 3]), color=color)
# Plot total line model
ax.plot(q.wave, q.Manygauss(np.log(q.wave), gauss_result), 'b', lw=2)
ax.plot(q.wave, q.line_flux, 'k')
ax.set_xlim(6400, 6800)
ax.set_xlabel(r'$\rm Rest \, Wavelength$ ($\rm \AA$)', fontsize=20)
ax.set_ylabel(r'$\rm f_{\lambda}$ ($\rm 10^{-17} erg\;s^{-1}\;cm^{-2}\;\AA^{-1}$)', fontsize=20)
c = 1 # Ha
ax.text(0.02, 0.80, r'$\chi ^2_\nu=$' + str(np.round(float(q.comp_result[c * 7 + 4]), 2)),
fontsize=16, transform=ax.transAxes)
plt.show()
Fitting broad H$\alpha$ with 1 components. Delta BIC = inf
Fitting broad H$\alpha$ with 2 components. Delta BIC = 83.0
Fitting broad H$\alpha$ with 3 components. Delta BIC = 31.8
Fitting broad H$\alpha$ with 4 components. Delta BIC = 8.1 3 components is prefered
Multiproccessing example:
Uncertainty computation will cost some time, so we usually use the multiprocessing to speed up large sample fitting.
We note that running multiprocessing on jupyter notebook may lead to various errors depending on the environment. The example here is only tested in MacOS Big Sur 11.2 with python 3.8.5. Running multiprocessing in non-interactive interpreters is recommended. For more information see https://docs.python.org/3/library/multiprocessing.html
Memory leak warning: For the convenience of display, the code wouldn't close the figure automatically if the plot_fig=True
and save_fig=False
are set, instead, we leave the figure as a variable self.fig
for users. Such case could lead to memory leak by child process if one uses multiprocessing. To avoid memory leak, please set plot_fig
and save_fig
to the same (Both True or False is OK, depending on your need).
def job(file):
print(file + '\n')
data = fits.open(file)
lam = 10 ** data[1].data['loglam']
flux = data[1].data['flux']
err = 1 / np.sqrt(data[1].data['ivar'])
ra = data[0].header['plug_ra']
dec = data[0].header['plug_dec']
z = data[2].data['z'][0]
plateid = data[0].header['plateid']
mjd = data[0].header['mjd']
fiberid = data[0].header['fiberid']
q = QSOFit(lam, flux, err, z, ra=ra, dec=dec, plateid=plateid, mjd=mjd, fiberid=fiberid, path=path_ex)
q.Fit(name=None, nsmooth=1, deredden=True, reject_badpix=False, wave_range=None, wave_mask=None,
decompose_host=True, npca_gal=5, npca_qso=20, Fe_uv_op=True, poly=True, BC=False, rej_abs_conti=False,
rej_abs_line=False,
MCMC=False, plot_fig=False, save_fig=False, save_result=True,
kwargs_plot={'save_fig_path': '.', 'broad_fwhm': 1200},
save_fits_path=path_out, save_fits_name=None, verbose=False)
# Edit the directory before use
from multiprocess import Pool
import os, timeit, glob
if __name__ == '__main__':
start = timeit.default_timer()
files = glob.glob(os.path.join(path_ex, 'data/spec-*.fits'))
pool = Pool(3) # Create a multiprocessing Pool
pool.imap(func=job, iterable=files)
end = timeit.default_timer()
print(f'Fitting finished in : {np.round(end - start)}s')
Fitting finished in : 0.0s ./data/spec-0266-51602-0107.fits ./data/spec-0266-51602-0013.fits ./data/spec-0332-52367-0639.fits
This code can also fit emission line galaxies.
# Remove the broad components from the line parameter file
mask_na = [True if not '_br' in str(row['linename']) else False for row in hdu1.data]
# Or you can set broad component number to zero
# Save line info
hdu1 = fits.BinTableHDU(data=hdu1.data[mask_na], header=hdr1, name='data')
hdu_list = fits.HDUList([primary_hdu, hdu1, hdu2, hdu3, hdu4])
hdu_list.writeto(os.path.join(path_ex, 'qsopar.fits'), overwrite=True)
data = fits.open(os.path.join(path_ex, 'data/spec-0266-51602-0107.fits'))
lam = 10 ** data[1].data['loglam'] # OBS wavelength [A]
flux = data[1].data['flux'] # OBS flux [erg/s/cm^2/A]
err = 1 / np.sqrt(data[1].data['ivar']) # 1 sigma error
z = data[2].data['z'][0] # Redshift
# Optional
ra = data[0].header['plug_ra'] # RA
dec = data[0].header['plug_dec'] # DEC
plateid = data[0].header['plateid'] # SDSS plate ID
mjd = data[0].header['mjd'] # SDSS MJD
fiberid = data[0].header['fiberid'] # SDSS fiber ID
q = QSOFit(lam, flux, err, z, ra=ra, dec=dec, plateid=plateid, mjd=mjd, fiberid=fiberid, path=path_ex)
q.Fit(name=None, nsmooth=1, deredden=True, reject_badpix=False, wave_range=None, \
wave_mask=None, decompose_host=True, host_prior=True, decomp_na_mask=True, npca_gal=5, npca_qso=10, qso_type='CZBIN1',\
Fe_uv_op=False, poly=False, BC=False, rej_abs_conti=False, rej_abs_line=False, MCMC=False, \
save_result=False, plot_fig=True, save_fig=False, save_fits_path=path_out, save_fits_name=None)
data = fits.open(os.path.join(path_ex, 'data/spec-0266-51602-0013.fits'))
lam = 10 ** data[1].data['loglam'] # OBS wavelength [A]
flux = data[1].data['flux'] # OBS flux [erg/s/cm^2/A]
err = 1 / np.sqrt(data[1].data['ivar']) # 1 sigma error
z = data[2].data['z'][0] # Redshift
# Optional
ra = data[0].header['plug_ra'] # RA
dec = data[0].header['plug_dec'] # DEC
plateid = data[0].header['plateid'] # SDSS plate ID
mjd = data[0].header['mjd'] # SDSS MJD
fiberid = data[0].header['fiberid'] # SDSS fiber ID
q = QSOFit(lam, flux, err, z, ra=ra, dec=dec, plateid=plateid, mjd=mjd, fiberid=fiberid, path=path_ex)
q.Fit(name=None, nsmooth=1, deredden=True, reject_badpix=False, wave_range=None, \
wave_mask=None, decompose_host=True, host_prior=True, decomp_na_mask=True, npca_gal=5, npca_qso=10, qso_type='CZBIN1', linefit=False, \
Fe_uv_op=False, poly=False, BC=False, MCMC=False, \
save_result=False, plot_fig=True, save_fig=False, save_fits_path=path_out, save_fits_name=None)