<p>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 <code>log_rebin</code> function, which computes and returns this (among other things).</p>
<p>Your data and template spectra need to have the <em>same</em> spacing in <span class="math-container">$\ln \lambda$</span> (your <code>loglam_gal</code> and <code>loglam_temp</code>): in other words, they have to have the <em>same</em> <code>velscale</code>. (It's not clear from your code if this is the case.)</p>
<p>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 <code>vsyst</code> keyword. I usually compute it this way (using linear wavelength vectors):</p>
<pre><code>vsyst = c_kms * math.log(lambda_template[0]/lambda_galaxy[0]) .
</code></pre>
<p>Using your variables, that would probably be something like</p>
<pre><code>vsyst = c_kms * (loglam_temp[0] - loglam_gal[0])
</code></pre>
<p>You don't need to use &quot;the same rest-frame piece&quot;. The template spectrum should be rest-frame, otherwise you won't get good velocity measurements. The data spectrum should be as observed. The <code>vsyst</code> keyword is meant to take care of any differences between the initial wavelengths of the (observed) data and the template.</p>
<p>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.</p>
<p>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.</p>
<p><strong>Edited to add:</strong></p>
<p>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 <em>not</em> 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.)</p>
<pre><code>Best Fit:       Vel     sigma        h3        h4
comp.  0:      1327       159     0.017     0.037
</code></pre>
<p>This is the figure automatically generated by pPXF (<code>ppxf(..., plot=True, ...)</code>), 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).</p>
<p><a href="https://i.sstatic.net/p7ypN.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/p7ypN.png" alt="enter image description here" /></a></p>
<p>Black = data, red = best-fit model spectrum (weighted combination of 25 template spectra); green = residuals. Note that there is a weak emission line at <span class="math-container">$\sim 220$</span> pixels, which is redshifted [O III] 5007; H<span class="math-container">$\beta$</span> is probably also present in emission, though weaker.</p>
<p>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 &quot;fit using a single template spectrum&quot; 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.</p>
<pre><code>Best Fit:       Vel     sigma        h3        h4
comp.  0:      1328       193     0.020     0.112
</code></pre>
<p><a href="https://i.sstatic.net/sKR2i.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/sKR2i.png" alt="enter image description here" /></a></p>
<p>[Data = VLT-MUSE observation of NGC 4643 from TIMER Project (<a href="http://adsabs.harvard.edu/abs/2019MNRAS.482..506G" rel="nofollow noreferrer">Gadotti et al. (2019)</a>; reanalyzed in <a href="https://ui.adsabs.harvard.edu/abs/2021MNRAS.502.2446E" rel="nofollow noreferrer">Erwin et al. (2021).</a>]</p>

CC BY-SA 4.0
2023-05-10T08:29:19.410Z
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?
trynerror
CC BY-SA 4.0
2023-05-10T08:30:16.427Z
When I fit only v and sigma, my results actually get worse ( they are systematically even lower than the required results).
trynerror
CC BY-SA 4.0
2023-05-10T08:36:10.580Z
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.
trynerror
CC BY-SA 4.0
2023-05-10T09:10:37.457Z
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...?).
Peter Erwin
CC BY-SA 4.0
2023-05-10T09:37:08.890Z
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.
Peter Erwin
CC BY-SA 4.0
2023-05-10T09:37:56.723Z
I've also updated my answer to include examples of good (if not necessarily best) fits, showing how the convolved, shifted template spectrum should match the data.
Peter Erwin
CC BY-SA 4.0
2023-05-10T09:40:39.510Z
I wouldn't worry about the constant noise value. Ideally, you want reasonably accurate noise estimates -- and you'll have to rescale them properly if you normalize the galaxy spectrum as you are doing -- but that's a second- or third-order issue (the fits I show in my updated answer used constant noise values).
Peter Erwin
CC BY-SA 4.0
2023-05-10T11:47:30.253Z
I included a plot of the (normalized) spectra before I apply the ppxf method. Are the strength of the template lines to high, since they do not really have the same depth? After the ppxf method the galaxy model obtained appears to have much less deep lines then the observed galactic spectrum. I thought this change of line strengths could have been due to the galaxy model of the ppxf algorithm (is this maybe connected to the low template weigth of 0.272?) The central line I used has S/N > 50, so it should be enough for the kinematic analysis?
trynerror
CC BY-SA 4.0
2023-05-10T11:54:46.310Z
Maybe there is a way to force ppxf to use template weigth 1, since it is the only template I want to use for comparison?
trynerror
CC BY-SA 4.0
2023-05-09T22:18:21.277Z
2023-05-10T09:42:43.470Z
2023-05-10T09:42:43.470Z
Peter Erwin
<p>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?</p>
<pre><code>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
</code></pre>
<p>And the plot:<a href="https://i.sstatic.net/jVd8d.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/jVd8d.png" alt="enter image description here" /></a></p>
<p>My spectra before the plot look like this, only applying the normalization</p>
<pre><code>flux_gal /= np.median(flux_gal, 0)
flux_temp /= np.median(flux_temp, 0)
</code></pre>
<p>Plot:
<a href="https://i.sstatic.net/T3fnL.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/T3fnL.png" alt="enter image description here" /></a></p>
<p>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.</p>
<p>Manually changing the degrees of the multiplicative and additive polynomials gives a much better fit, but still not perfect:</p>
<pre><code>pp = ppxf(flux_temp, flux_gal, noise, velscale, start, moments=4, plot=True, mask=mask, mdegree=5,degree=5,vsyst=dv)
</code></pre>
<p>With the better results:
<a href="https://i.sstatic.net/HPBEn.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/HPBEn.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2023-05-09T15:57:41.377Z
You have just the one template spectrum, is that right? (Could be an example of template mismatch; the fit looks pretty bad.)
Peter Erwin
CC BY-SA 4.0
2023-05-09T18:17:20.733Z
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?
trynerror
CC BY-SA 4.0
2023-05-10T15:38:55.787Z
Convolving the template spectrum with the LOSVD (mostly just convolving with a Gaussian) smooths it, making the absorption lines broader and weaker.
Peter Erwin
CC BY-SA 4.0
2023-05-10T15:40:36.910Z
What happens if you turn off the additive polynomial term? (`ppxf(…, degree=-1, …)`)
Peter Erwin
CC BY-SA 4.0
2023-05-10T18:31:07.720Z
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 ?
trynerror
CC BY-SA 4.0
2023-05-10T18:45:44.057Z
To be honest, I usually don't use additive polynomials. I've found multiplicative degree=8 useful for MUSE spectra, but that's for a fairly large wavelength range (4800-7200 Angstrom, usually).
Peter Erwin
CC BY-SA 4.0
2023-05-10T18:48:37.247Z
If I dont use additive polynomials the fit becomes much worse. Isn't that indicating that it actually is a case of template mismatch, since the additive polynomials correct for that ? Seem to get the best results for both polynomial type of degree 11
trynerror
CC BY-SA 4.0
2023-05-10T18:49:55.217Z
which ppxf output can I use to actually quantify the quality of the fit and compare?
trynerror
<observational-astronomy><galaxy><data-analysis><spectra><kinematics>
What am I doing wrong in this ppxf application?
trynerror
# What am I doing wrong in this ppxf application?

(Asked by: username_0 on 2023-05-09)

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:
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:

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:

**Tags:** observational-astronomy, galaxy, data-analysis, spectra, kinematics
> **Comments:**
> **username_1**: You have just the one template spectrum, is that right? (Could be an example of template mismatch; the fit looks pretty bad.)  
> **username_0**: 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?  
> **username_1**: Convolving the template spectrum with the LOSVD (mostly just convolving with a Gaussian) smooths it, making the absorption lines broader and weaker.  
> **username_1**: What happens if you turn off the additive polynomial term? (`ppxf(…, degree=-1, …)`)  
> **username_0**: 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 ?  
> **username_1**: To be honest, I usually don't use additive polynomials. I've found multiplicative degree=8 useful for MUSE spectra, but that's for a fairly large wavelength range (4800-7200 Angstrom, usually).  
> **username_0**: If I dont use additive polynomials the fit becomes much worse. Isn't that indicating that it actually is a case of template mismatch, since the additive polynomials correct for that ? Seem to get the best results for both polynomial type of degree 11  
> **username_0**: which ppxf output can I use to actually quantify the quality of the fit and compare?  

---

## 1 answers

### Answer by username_1 on 2023-05-09. Score: 1 | Accepted answer

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).

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


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

> **Comments:**
> **username_0**: 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?  
> **username_0**: When I fit only v and sigma, my results actually get worse ( they are systematically even lower than the required results).  
> **username_0**: 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.  
> **username_1**: 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...?).  
> **username_1**: 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.  
> **username_1**: I've also updated my answer to include examples of good (if not necessarily best) fits, showing how the convolved, shifted template spectrum should match the data.  
> **username_1**: I wouldn't worry about the constant noise value. Ideally, you want reasonably accurate noise estimates -- and you'll have to rescale them properly if you normalize the galaxy spectrum as you are doing -- but that's a second- or third-order issue (the fits I show in my updated answer used constant noise values).  
> **username_0**: I included a plot of the (normalized) spectra before I apply the ppxf method. Are the strength of the template lines to high, since they do not really have the same depth? After the ppxf method the galaxy model obtained appears to have much less deep lines then the observed galactic spectrum. I thought this change of line strengths could have been due to the galaxy model of the ppxf algorithm (is this maybe connected to the low template weigth of 0.272?) The central line I used has S/N > 50, so it should be enough for the kinematic analysis?  
> **username_0**: Maybe there is a way to force ppxf to use template weigth 1, since it is the only template I want to use for comparison?  

---

