HeloiseS / hoki

Bridging the gap between observation and theory
https://heloises.github.io/hoki/intro.html
BSD 3-Clause "New" or "Revised" License
47 stars 8 forks source link

[BUG] Template convolution to instrumental resolution incorrect when `dispersion_obs` is used #105

Open Knusper opened 4 months ago

Knusper commented 4 months ago

Reading the pPXF papers and documentation, I understand that the templates have to be convolved with a Kernel that matches the instrumental resolution (line spread function) of the observations to that of the template. It is established practice, that both are described by Gaussians, so that the matching kernel is also a Gaussian with sigma = sqrt(sigma_ins^2 - sigma_tem^2), where both sigmas are wavelength dependent (see Sect. 2.2 in Cappelari 2017). (Note also, that for higher redshifts we have to de-redshift the galaxy spectrum, and thus also compress the LSF width accordingly.)

In another bug report, I stated that the sigma_tem for the BPASS models should be 0, since these are synthetic models. Perhaps there is a confusion between spectral sampling and LSF width, since also the tutorial of the MUSE data uses fwhm_obs=1.25, which is the spectral sampling of the MUSE data - not the LSF width. The latter is roughly 2 sampling elements, but varies as a function of wavelength (see, e.g., Sect. 4.2.2 in Bacon et al. 2023.)

Now, in l. 408 of the kvn.py I found an strange expression

self.sigma = self.fwhm_dif/2.355/self.fwhm_tem

Of course, this does not matter at all, since you fixed fwhm_tem to 1. But there is no need for this extra division. This is just cosmetics, but it made me look deeper...

fwhm_dif is calculated earlier (l.361ff of kvn.py):

       if dispersion_obs is not None:
            if verbose: print(f"{dialogue.running()} Calulating obs. velocity scale")
            # Calc. velocity scale
            _frac = self.wl_obs[1]/self.wl_obs[0]
            _dwl_obs = (_frac - 1)*self.wl_obs
            self.velscale_obs = np.log(_frac)*c #speed of light # Velocity scale in km/s per pixel

            if verbose: print(f"{dialogue.running()} Calulating FWHM")
            # Calc. FWHM galaxy (res.)
            _fwhm_obs = 2.355*dispersion_obs*_dwl_obs
            self.fwhm_obs = np.interp(self.wl_tem, self.wl_obs, _fwhm_obs)
            self.fwhm_dif = np.sqrt((self.fwhm_obs**2 - self.fwhm_tem**2).clip(0))

        elif dispersion_obs is None and fwhm_obs is not None:
            self.wl_range_obs = wl_range_obs
            self.fwhm_obs = fwhm_obs
            if verbose: print(f"{dialogue.running()} Calulating obs. velocity scale -- No dispersion")
            __, __, self.velscale_obs = ppxf_util.log_rebin(self.wl_range_obs, self.wl_obs)
            #self.velscale_obs=c*np.log(self.wl_obs[1]/self.wl_obs[0])
            self.fwhm_dif = np.sqrt((self.fwhm_obs**2 - self.fwhm_tem**2))

Either dispersion_obs or fwhm_obs are given in Angstrom, but incorrectly dispersion_obs is converted to pixel before subtracting the template resolution in quadrature. This conversion is not made, when supplying a single value fwhm_obs.