0
$\begingroup$

I am trying to derive the kinematics of a galaxy spectra. To test my results I want to implement the Penalized Pixel Fitting method. For the galaxy spectrum I use a central line with sufficient S/N. I cut of part of the edges to use the same part of the spectra I used with a different method. The following code shows my initial try, but I am unsure about the definition of the velscale keyword or maybe I have to work in the same rest frame piece for galaxy and template? My obtained results, seem to have too low sigma and too large negative h3 and h4 (by comparing with results from other papers). Is there anything wrong about the way I tried to implement this?

import ppxf as ppxf_package
from ppxf.ppxf import ppxf
import ppxf.ppxf_util as util
import matplotlib.pyplot as plt


def my_ppxf():
    # data access
    # template
    temp_data = np.loadtxt('C:/Users/reich/OneDrive/Dokumente/Uni/Bachelorarbeit/Python/bins/template.txt')
    loglam_temp = temp_data[0, :]
    flux_temp = temp_data[1, :]
    # galaxy central line
    gal_data_center = np.loadtxt('C:/Users/reich/OneDrive/Dokumente/Uni/Bachelorarbeit/Python/bins/0_central.txt')
    loglam_gal = gal_data_center[0, :]
    flux_gal = gal_data_center[1, :]

    # 1.) cut off pixels at both ends off arrays before removal of continuum
    # subtract 1 to get uneven number of pixels
    cut_off_pixel = 90
    flux_temp = flux_temp[cut_off_pixel:-(cut_off_pixel - 1)]
    flux_gal = flux_gal[cut_off_pixel:-(cut_off_pixel - 1)]
    loglam_gal = loglam_gal[cut_off_pixel:-(cut_off_pixel - 1)]
    loglam_temp = loglam_temp[cut_off_pixel:-(cut_off_pixel - 1)]

    # 2.) normalize spectra
    flux_gal /= np.median(flux_gal, 0)
    flux_temp /= np.median(flux_temp, 0)

    c = 299792.458                  # speed of light in km/s
    frac = loglam_gal[1]-loglam_gal[0]    # Constant lambda fraction per pixel
    velscale = c*frac

    noise = np.full_like(flux_gal, 0.003)  # Assume constant noise per pixel here

    z = 0.00414811      #galaxy redshift

    # Here the actual fit starts. The best fit is plotted on the screen.
    # set initial guesses for the parameters, h3 and h4 by default 0
    vel = c*np.log(1 + z)   # eq.(8) of Cappellari (2017)
    start = [vel, 200.]  # (km/s), starting guess for [V, sigma]

    pp = ppxf(flux_temp, flux_gal, noise, velscale, start, moments=4, plot=True)
    # Print the results
    print('V:', pp.sol[0])
    print('Sigma:', pp.sol[1])
    print('h3:', pp.sol[2])
    print('h4:', pp.sol[3])
    plt.show()

    # Plot the best-fit model and the residuals
    plt.figure()
    plt.plot(loglam_gal, flux_gal, label='Galaxy spectrum')
    plt.plot(loglam_gal, pp.bestfit, label='Best-fit model')
    plt.plot(loglam_gal, flux_gal - pp.bestfit, label='Residuals')
    plt.xlabel('Wavelength (Angstroms)')
    plt.ylabel('Flux')
    plt.legend()
    plt.show()

#------------------------------------------------------------------------------

if __name__ == '__main__':

    my_ppxf()```

The output:
```Templates weights:
     0.272
V: 1233.8063289282345
Sigma: 127.00695685589174
h3: -0.06568111895579586
h4: -0.3

And the plot:enter image description here

My spectra before the plot look like this, only applying the normalization

flux_gal /= np.median(flux_gal, 0)
flux_temp /= np.median(flux_temp, 0)

Plot: enter image description here

So here the spectral lines of the template are actually stronger than the ones of the galaxy. After the ppxf fit they appear much weaker, I thought this has to do something with the galaxy models used by ppxf and the corresponding polynomials, that can change the line strengths.

Manually changing the degrees of the multiplicative and additive polynomials gives a much better fit, but still not perfect:

pp = ppxf(flux_temp, flux_gal, noise, velscale, start, moments=4, plot=True, mask=mask, mdegree=5,degree=5,vsyst=dv)

With the better results: enter image description here

$\endgroup$
8
  • $\begingroup$ You have just the one template spectrum, is that right? (Could be an example of template mismatch; the fit looks pretty bad.) $\endgroup$ May 9, 2023 at 15:57
  • $\begingroup$ Yes i am only using one template to compare my results to a different method. I know that the influence of template mismatch can not be that bad, because with this template-galaxy combination much better fit results have been obtained. So it must be something about the Implementation. I was not sure if I have to use the same rest frame piece for galaxy and template for the ppxf method, could it maybe be because of that? $\endgroup$
    – trynerror
    May 9, 2023 at 18:17
  • $\begingroup$ Convolving the template spectrum with the LOSVD (mostly just convolving with a Gaussian) smooths it, making the absorption lines broader and weaker. $\endgroup$ May 10, 2023 at 15:38
  • 1
    $\begingroup$ What happens if you turn off the additive polynomial term? (ppxf(…, degree=-1, …)) $\endgroup$ May 10, 2023 at 15:40
  • $\begingroup$ I played a bit with both the degrees of multiplicative and additive polynomials, siginificantly increasing the quality of the fit. At the moment I chose 5 for both. The fit is good in the central region, but becomes worse in the outer regions. Are both of the polynomial types sensitive for that region ? $\endgroup$
    – trynerror
    May 10, 2023 at 18:31

1 Answer 1

1
$\begingroup$

The velscale keyword defines the spacing of the spectrum's wavelength vector, in velocity (km/s) units. Look in the ppxf code (in ppxf.ppxf_util) for the log_rebin function, which computes and returns this (among other things).

Your data and template spectra need to have the same spacing in $\ln \lambda$ (your loglam_gal and loglam_temp): in other words, they have to have the same velscale. (It's not clear from your code if this is the case.)

If the data and template spectra start at different (log) wavelengths, you need to compute this difference in the form of a velocity offset, and provide that value to ppxf in the vsyst keyword. I usually compute it this way (using linear wavelength vectors):

vsyst = c_kms * math.log(lambda_template[0]/lambda_galaxy[0]) .

Using your variables, that would probably be something like

vsyst = c_kms * (loglam_temp[0] - loglam_gal[0])

You don't need to use "the same rest-frame piece". The template spectrum should be rest-frame, otherwise you won't get good velocity measurements. The data spectrum should be as observed. The vsyst keyword is meant to take care of any differences between the initial wavelengths of the (observed) data and the template.

I would also suggest starting off simply: don't try to fit for all four moments, just fit for velocity and dispersion when starting off.

I asked about template mismatch because in the plot the best-fit template spectrum looks very different from the data, as though it's been vertically compressed -- e.g., the absorption lines are much shallower than they are in the data.

Edited to add:

Here's an example of what a fit should look like. This is an S0/a galaxy spectrum from a VLT-MUSE datacube, using most of the cube added together rather than a single spaxel -- so the S/N is high. I did not normalize the data spectrum; it's really not necessary (and makes it easier to go back and redo the fit with proper noise estimates if you don't normalize). (The template spectra were normalized to begin with.)

Best Fit:       Vel     sigma        h3        h4
comp.  0:      1327       159     0.017     0.037

This is the figure automatically generated by pPXF (ppxf(..., plot=True, ...)), trimmed down to approximately the same wavelength range as your spectrum (pixel 0 = 8.4765, pixel 800 = 8.6467 in natural-log of wavelength in Angstroms).

enter image description here

Black = data, red = best-fit model spectrum (weighted combination of 25 template spectra); green = residuals. Note that there is a weak emission line at $\sim 220$ pixels, which is redshifted [O III] 5007; H$\beta$ is probably also present in emission, though weaker.

This is the same thing, but allowing pPXF to use only the single highest-weighted template spectrum from the previous fit -- this is more analogous to your "fit using a single template spectrum" approach. It's clearly a worse fit (and sigma is too high), but it doesn't show the mismatch in scaling (i.e., template absorption lines too weak) that appears in your plot.

Best Fit:       Vel     sigma        h3        h4
comp.  0:      1328       193     0.020     0.112

enter image description here

[Data = VLT-MUSE observation of NGC 4643 from TIMER Project (Gadotti et al. (2019); reanalyzed in Erwin et al. (2021).]

$\endgroup$
9
  • $\begingroup$ My fluxes and wavelengths are logarithmically rebinned so my keyword frac is the same for all differences loglam_gal[x_(n+1)]-loglam_gal[x_n]``. I included the vsyst` keyword, but I do not understand why you define it as log of already log functions? I think maybe one problem is that my template weight is 0.227, so the influence of the additive and multiplicative polynomials is much higher, which could in turn look like template mismatch? How can I force ppxf to have template weight one, since it is only one template and I know that it should be a good stellar model for the galaxy? $\endgroup$
    – trynerror
    May 10, 2023 at 8:29
  • $\begingroup$ When I fit only v and sigma, my results actually get worse ( they are systematically even lower than the required results). $\endgroup$
    – trynerror
    May 10, 2023 at 8:30
  • $\begingroup$ Furthermore, I am not sure about the construction of the constant noise value. How can I calculate this (what parameters do I need to estimate it?). I only calculated the S/N of this spectrum. $\endgroup$
    – trynerror
    May 10, 2023 at 8:36
  • $\begingroup$ I would maybe concentrate on the template spectrum, which does not look like a good match to the galaxy spectrum (the absorption lines should have the same depth!). Admittedly, I have never tried using pPXF with just a single template spectrum, so I don't now if that's part of the problem (probably not...?). $\endgroup$ May 10, 2023 at 9:10
  • $\begingroup$ OK, you're right about the computation of vsyst; I was copying a snipped from my own Python code, which uses the linear wavelength vector. I've updated my answer to correct this. $\endgroup$ May 10, 2023 at 9:37

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .