Radial velocity from Gaussian fit not as expected

[Knusper]

According to https://lime-stable.readthedocs.io/en/latest/outputs/outputs1_measurements.html
the v_r output parameter contains the radial velocity (note that there are brackets missing in the definition).

However, the following puzzles me, since I don't get the same radial velocity when using the manual formula and the lime output (the pickle dictionary with spectrum, variance, and wavelength axis is attached).

from astropy import constants as cons
import numpy as np
import pickle
import lime

c = cons.c.to('km/s').value

with open('em_line_fit_debug.pickle', 'rb') as f:
    in_dic = pickle.load(f)

xax0 = in_dic['xax0']  # wavelength axis (rest-frame / air wavelengths)
ln_spec_in = in_dic['ln_spec_in']  # input spectrum (continuum subtracted)
ln_spec_var = in_dic['ln_spec_var']  # variance

# lime bands (for test just Balmer-series & [OIII] doublet
particle_list = ('H1', 'O3')
bands = lime.line_bands(wave_intvl=(xax0[0], xax0[-1]),
                        particle_list=particle_list,
                        units_wave='Angstrom', decimals=None, vacuum=False)

# lime spectrum
lspec = lime.Spectrum(input_wave=xax0, input_flux=ln_spec_in, input_err=np.sqrt(ln_spec_var),
                      units_flux='1e-20*FLAM', redshift=0)
lspec.fit.continuum(degree_list=[3,5,7], emis_threshold=[3,3,3], plot_steps=False)
lspec.line_detection(bands, sigma_threshold=5, plot_steps=False)
lspec.fit.frame(bands=bands)

# check Balmer alpha line fit
l = lspec.frame.loc[['O3_5007A']]

# radial velocity from LiME
v_r_lime = l.v_r.iloc[0]
v_r_manual = (l.center.iloc[0] / l.wavelength.iloc[0] - 1 ) * c 

print('\n')
print('v_r_lime:' + str(v_r_lime))
print('v_r_manual:' + str(v_r_manual))

The output is:

Line fitting progress:
[==========] 100% of 4 lines (H1_6563A)

v_r_lime:-18.850342740006468
v_r_manual:6.74700738846691

Now I would have expected that v_r_manual and v_r_lime are the same - v_r_manual is in this example also the value that appears correct given a previous result on this spectrum.

What is going on here?

(Btw. is there a way to turn of the "line fitting progress" message...)

em_line_fit_debug.pickle.gz

[Vital-Fernandez]

Thank you @Knusper very much for your comment. I have checked your code and there is indeed a difference.

This has to do with the reference wavelength used to calculate the radial velocity. In the case of a single lines, it is the peak wavelength. In your code this would be:

v_r_manual = (l.peak_wavelength.iloc[0] / l.wavelength.iloc[0] - 1 ) * c

where l.peak_wavelength is the wavelength of the maximum flux pixel on the line band, which depending on the spectral resolution may be different from the central wavelength.

In complex blended profiles, specially in IFU data sets, the intensity of the peak may change considerably due to ionization conditions rather than kinematic... which may make it misleading to keep track of the kinematics using the peak wavelength... so Instead I use the theoretical leading line wavelength times (1 + z_input).

I would rather have the same logic for every line type... but I believe most users using longslist spectra prefer to use the peak wavelength and remove the redshift dependence...

Do you think the documentation should be more clear? It maybe good to add an sketch....