3. Spectral line fitting and radial velocity measurements

The main purpose of MUSEreduce is the spectral line fitting in order to create template spectra from observations, which are used to measure the radial velocities of stars and gas, especially in the absence of working spectral template libraries, e.g., for pre-main-sequence stars in the optical.

The radial_velocities.RV_spectrum class and the spectral fitting module line_fitter are the heart of this package using pyspeckit as basic spectral fitting routine and ppxf for the spectral cross-correlation. Together with a Monte Carlo bootstrap technique, it is possible to measure radial velocities with an accuracy of \(1-3\,\rm{km}\,\rm{s}^{-1}\) in the absence of any template library. For a detailed description on how this technique works, we refer to Zeidler et al. 2019.

3.1. Spectral line fitter

This module is the work horse for fitting spectral lines from an input catalog to a given 1D spectrum. The line fitting itself is performed using the module pyspeckit.

pyspeckit uses an iterative method to fit the spectral lines and the continuum simultaneously. line_fitter automatically adjusts the input parameters to pyspeckit in the manner to optimize the fitting (continuum order and wavelegth limits)

The line_fitter determines whether the fit fullfils the set parameters and was succesfull. The line_fitter also allows to fit blends by fixing a maximum ratio between the primary line (the one the user is interested in) and the secondary line (blend). This ensures that the blend does not become the dominant line.

In order to accomodate hyper-velocity stars an option is given that the wavelength limits are automatically adjust for each iteration \(n\) based on the solution \(n-1\) of the primary line taking into account \(\Delta \lambda / \lambda\).

If the spectra with high radial velocities, or if the object is red shifted an initial guess of the radial velocity needs to be provided. This may be done by the Kwargs rv_sys of the radial_velocities.RV_spectrum. This is important if the shift is much larger than the wavelength limits, where the automatic adjustment fails. Additionally, the more accurate this rv_sys is provided the faster the fit converges.

line_fitter.line_fitter(self, linecat, line_idx, niter, input_resid_level, max_contorder, max_ladjust, adjust_preference, input_continuum_deviation, llimits, max_exclusion_level, blends, autoadjust, fwhm_block)

This is module fits the spectral lines using pyspeckit. It automatically determines the goodness of fit and decides the best solution for the line and continuum fit. It handles blended lines in a way a maximum lines ratio can not be exceeded to not make the weaker blend the dominant lines. Additionally, the limits in wavelength range in iteration \(n\) can be adjusted automatically based on iteration \(n\) to accomodate for larger wavelength shifts.

In this way it is possible to fit spectral lines using a line catalog as only input.

Args:

linecatnumpy.array()

Array with the spectral lines and their wavelengths

line_idxstr

Name of the primary line

niterint

Number of iterations

input_resid_levelfloat

The maximum MAD for the fit residuals for a succesfull fit

max_contorderint

The maximum polynominal order of the continuum to have

max_ladjuststr

The maximum number of wavelength range adjustments in steps of 5 Angstrom

adjust_preferencestr

contorder: continuum order is adjusted first

wavelength: wavelength range is adjusted first

input_continuum_deviationfloat

by how much the continuum is allowed to deviate from a running median estimate. This is set to prevent lines mimicking a continuum

llimitslist

the limits for the wavelength fit as set in ppxf

max_exclusion_levelfloat

The exclusion level for lines to be excluded from the next baseline estimate as set in pyspeckit

blendsascii-file or None

A file with primary lines that contain blends to provide a maximum amplitude ratio of the primary and the blend to prevent that the blend becomes the dominant line in the fit.

autoadjustbool

True: the wavelength limits llimit will be adjusted to the fit of the previous iteration. All other wavelength range are adjusted accordingly taking into account the proper velocity corrected shift \(\Delta \lambda/\lambda\). This is especially important to detect hyper-velocity stars.

False: no adjustment to the limits done

fwhm_blockbool:obj:

True: The minimum fwhm of the voigt profiles of the fitted lines is the instrument’s dispersion

False: The minimum fwhm of the voigt profiles of the fitted lines is zero

Returns:

line_idxstr

Name of the primary line

temp_lfloat

fitted wavelength of the primary line

temp_afloat

fitted amplitude of the primary line

temp_slfloat

fitted Lorentzian gamma of the primary line

temp_sgfloat

fitted Gaussian sigma of the primary line

spec_select_idxnumpy.array()

indices of the used part of the spectrum

template_fnumpy.array()

template spectrum shifted to the rest-frame

continuumnumpy.array()

the fitted continuum

lstartfloat

first used wavelength bin (might differ from input if it was adjusted during fitting)

lendfloat

last used wavelength bin (might differ from input if it was adjusted during fitting)

contorderint

The order of the polynominal used for the continuum

fit_farray

the fitted spectrum

significancefloat

the line strength over the continuum

fit_failedbool

True: if the fit failed for some reason. This line will be excluded from further analyses

fit_f_highresfloat

the fitted spectrum and the oversampling resolution

spec_select_idx_highresnumpy.array()

indices of the used part of the over sampled spectrum

template_f_highresnumpy.array()

over sampled template spectrum shifted to the rest-frame

continuum_highresnumpy.array()

the fitted over sampled continuum

Warning

The automatic handling of absorption and emission lines implemented in vers. 0.1.2 has not been tested yet.

3.2. Radial Velocities

This is the main class for measuring the radial velocities. A basic example is descibed in Examples. A detailed description can be found in Zeidler et al. 2019.

Throughout the document the primary line is the spectral line of interest for which the line fitting will be executed. Secondary lines are spectral lines in blends.

A detailed demonstration on how to use the radial velocity fitter is provided in Examples including template files.

3.2.1. The radial velocity fitter

class radial_velocities.RV_spectrum(spec_id, spec_f, spec_err, spec_lambda, loglevel='INFO', templatebins=100000, specbinsize=1.25, dispersion=2.4, linetype='absorption', rv_sys=0.0)

This class contains the 1D spectrum, as well as all fitting parameters and output values including the radial velocity catalog.

Args:

spec_idstr

unique identifier of the spectrum

spec_fnumpy.array()

flux array of the spectrum in erg/s/cm:math:^2/Angstrom

spec_errnumpy.array()

flux uncertainty array of the spectrum

spec_lambdanumpy.array()

wavelength_array of the spectrum in Angstrom

Kwargs:

loglevelstr (optional, default: INFO)

DEBUG: all functions run on a single core to obtain proper output in the correct order for proper debugging.

templatebinsfloat (optional, default: 100000)

Number of wavelength bins for the oversampled spectrum used to fit the spectral lines.

specbinsizefloat (optional, default: 1.25)

The spectral bin size. The default is set to fit the MUSE dataset

dispersionfloat (optional, default: 2.4)

The dispersion of the spectrograph. The default is set to the nominal MUSE instrument dispersion

linetype: str (optional, default: absorption)

absorption: All spectral lines are absorption lines

emission: All spectral lines are emission lines

both : spectral lines can be absorption or emission lines. This mode has not been tested yet !!!

rv_sysfloat (optional, default: 0)

systematic RV shift in km/s Should be provided for large velocity offsets or redshifted objects

catalog(initcat=None, save=False, load=None, printcat=False)

This module handles the catalog that holds the fit results of the primary lines.

Kwargs:

initcatascii (optional, default: None)

ascii file containing the primary lines format: name, lambda, start, end

savebool (optional, default: False)

True: The catalog is written to a file.

loadstr (optional, default: None)

Loads a catalog from catalog file.

printcatbool (optional, default: False)

True: The catalog is printed to the terminal

clean()

This modules resets all of the output. This must be executed before repeating the spectral fit without re-initiating the class.

line_fitting(input_cat, line_idxs, niter=5, n_CPU=-1, resid_level=None, max_contorder=2, max_ladjust=4, adjust_preference='contorder', input_continuum_deviation=0.05, llimits=[-2.0, 2.0], max_exclusion_level=0.3, blends=None, autoadjust=False, fwhm_block=False)

Initializing the line fitting by calling line_fitter

Args:

input_cat: numpy.array()

input spectral line catalog with wavelengths in Angstrom

line_idxs: numpy.array()

str numpy.array(): of the line names for which RV fits should be performed, must be identical to “name” in init_cat

Kwargs:

niterint (optional, default: 5)

The maximum number of iteration if convergence is not reached

n_CPUfloat (optional, default: -1)

Setting the number of CPUs used for the parallelization. If set to -1 all available system resources are used. Maximum number of CPUs is the number of spectral lines the fit is performed on.

resid_level: float or None (optional, default: None)

The maximum MAD for the fit residuals for a succesfull fit

max_contorderint (optional, default: 2)

The maximum polynominal order of the continuum

max_ladjustint (optional, default: 4)

Ahe maximum number of wavelength range adjustments in steps of 5 Angstrom

adjust_preference: str (optional, default: contorder)

contorder: continuum order is adjusted first

wavelength: wavelength range is adjusted first

input_continuum_deviation float (optional, default: 0.05)

Fraction by how much the continuum is allowed to deviate from a running median estimate. This is set to prevent lines mimicking a continuum

llimits: list (optional, default: [-2., 2.])

the limits for the wavelength fit as set in ppxf

max_exclusion_level float (optional, default: 0.3)

The exclusion level for lines to be excluded from the next baseline estimate as set in pyspeckit

blends: ascii or None (optional, default: None)

A file with primary lines that contain blends to provide a maximum amplitude ratio of the primary and the blend to prevent that the blend becomes the dominant line in the fit.

autoadjust bool (optional, default: False)

True: the wavelength limits llimit will be adjusted to the fit of the previous iteration. All other wavelength range are adjusted accordingly taking into account the proper velocity corrected shift \(\Delta \lambda/\lambda\). This is especially important to detect hyper-velocity stars.

fwhm_block bool (optional, default: False)

True: The minimum fwhm of the voigt profiles of the fitted lines ais the instrument’s dispersion. This prevents unphysical lines.

False: The minimum fwhm of the voigt profiles of the fitted lines is zero.

plot(oversampled=False)

Plots the spectrum and the regions of the primary lines to a file including the spectral fit and the template.

Kwargs:

oversampledbool (optional, default: False)

True: The oversampled spectral fit and the template are used in the plot.

rv_fit(guesses, niter=10000, line_sigma=3, n_CPU=-1, line_significants=5, RV_guess_var=0.0)

The module runs the radial velocity fit using ppxf and the Monte Carlo radial velocity determination.

Args:

guessesnumpy.array()

The initial guesses for the the radial velocity fit guesses in the form [RV,sepctral_dispersion]

Kwargs:

niterint (optional, default: 10000)

number of iterations to bootstrap the spectrum

line_sigma: int (optional, default: 3):

sigma for the RV clipping for the individual lines

n_CPUfloat (optional, default: -1)

Setting the number of CPUs used for the parallelization. If set to -1 all available system resources are used. Maximum number of CPUs is the number of spectral lines the fit is performed to.

line_significants: int (optional, default: 5)

The sigma-level for the spectral line to be above the continuum in order to be considered valid

RV_guess_varfloat (optional, default: 0)

The maximum variation the RV guess will be varied using a uniform distribution.

rv_fit_peak(line_sigma=3, line_significants=5)

This module determines the radial velocity solely on the fitted peaks of the spectral lines

Kwargs:

line_sigmaint (optional, default: 3)

The sigma-level for the sigma clipping before determining peak radial velocity measurement

line_significants: int (optional, default: 5)

The sigma-level for the spectral line to be above the continuum in order to be considered valid

3.2.2. Monte Carlo radial velocity determination

This is the module that performs the Monte Carlo boot strapping to measure radial velocities. It is automatically called by the radial_velocities.RV_spectrum.rv_fit attribute of radial_velocities.RV_spectrum class and there is no need to repeat this step manually. Nevertheless, we show this part of the code since it may be useful for other applications to measure RVs.

ppxf_MC.ppxf_MC(log_template_f, log_spec_f, log_spec_err, velscale, guesses, nrand=100, degree=4, goodpixels=None, moments=4, vsyst=0, sigma=5, spec_id=None, RV_guess_var=0.0, n_CPU=-1)

This module runs the Monte Carlo ppxf runs, which is needed for the RV measurements. Most of the input parameters are similar to the standard ppxf parameters (see Cappellari and Emsellem 2004) for a more detailed explanation).

Args:

log_template_fnumpy.array()

The logarithmically binned template spectrum

log_spec_fnumpy.array()

The logarithmically binned source spectrum

log_spec_errnumpy.array()

The logarithmically source spectrum uncertainties

velscale: float

The velocity scale of the source spectrum

guessesnumpy.array()

The initial guesses for the the radial velocity fit guesses in the form [RV,sepctral_dispersion]

Kwargs:

nrandint (optional, default: 5)

The maximum number of iteration if convergence is not reached

degreeint (optional, default: 4)

The degree of the additive polynomial to fit offsets in the continuum. A low order polynominal may be needed to get better results

goodpixelsnumpy.array() (optional, default: 5)

A numpy.array(): of pixels that are used by ppxf to fit the template to the source spectrum

momentsint (optional, default: 4)

The moments to be fit (v, sigma, h3, h4)

vsystfloat (optional, default: 0.)

A systematic velocity. This may be needed if the system move at high velocities compared to the rest frame. If the guess of vsyst is good, the fit runs faster and more stable.

sigmaint (optional, default: 5)

The sigma used to clip outliers in the histogram determination.

spec_idstr (optional, default: None)

An ID number of the source spectrum. This becomes handy when fitting many individual sources because the output files will be named with the ID.

RV_guess_varfloat (optional, default: 0)

The maximum variation the RV guess will be varied using a uniform distribution.

n_CPUint (optional, default: -1)

The number of cores used for the Monte Carlo velocity fitting. If n_CPU=-1 than all available cores will be used.

3.3. History

Added in version 0.1.0: working radial velocity fitting package

Added in version 0.1.0: no RV calculation for lines that are below the significance level

Added in version 0.1.0: instrument dispersion added as keyword

Added in version 0.1.1: moved to pep-8

Added in version 0.1.2: now handles absorption and emission lines. Not tested yet, though

Added in version 0.1.3: During the boot strap procedure, the initial guesses of the velocities are varied between certain limits, to ensure more stability to the fit in the case of an “unlucky choice of initial parameters. Added Kwarg to RV_spectrum.rv_fit is RV_guess_var, which is by default 0. It describes the min/max variation of the initial RV guess for each fit.

Added in version 1.0: The release version as originally published in Zeidler et al. 2019.

Added in version 1.1: Introducing the Kwargs rv_sys, which should be used for large initial RV offsets from rest frame or for redshifted spectra. rv_sys should provide an estimate of that shift.