SED fitting of a transient neighbourhood

In this jupyter notebook we go over how to integrate the SED of the transient neighbourhood and remove the contribution from the transient spectrum before fitting with BPASS SEDs using hoki and ppxf.

For a kilonova the region around the transient is unlikely to be the region where the system was born so it is not all that relevant but you could do the same thing for a core collapse supernova in which case you transient is a lot more likely to be near its birth environment (although there are many caveats to that such as the type of CCSN and how distant the galaxy etc…).

# standard python libraries
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pickle

# ppxf stuff for fitting
from ppxf.ppxf import ppxf
import ppxf.ppxf_util as util

# hoki things for fitting
from hoki import load
from hoki import sedfitting as hsed

plt.style.use('../tuto.mplstyle')
c = 299792.458

1 - Separating Transient and Neighbourhood stars SEDs

1.1 - Neighbourhood SED

First let’s plot our galaxy and the situate our transient. I’ve already marked it and given you the pixel coordinate of the “centroid” (or as close as you can get it when the center is a 2 by 2 pixel section).

with open('data/cropped_NGC4993.pkl', 'rb') as pickle_file:
    cube = pickle.load(pickle_file)

f,ax = plt.subplots(nrows=1, ncols=1)
ima = cube.sum(axis=0)
ima.plot(scale='arcsinh', colorbar='v', ax=ax)
ax.scatter([49], [119], color='w', marker='*', label='AT 2017gfo')
ax.legend()

<matplotlib.legend.Legend at 0x7f3e1c2019a0>

A very important thing to note is that the # NOTE THAT THE X AND Y ARE REVERESED IN THE DATA COMPARED TO MATPLOTLIB. So when plotting and handling data we need to keep in mind the axes are transposed and all is well. If you forget it gets confusing quick.

Here we’re going to integrate the SEDs for the pixels that are between 1.5 and 2.5 arcseconds away from the centroid of the transient. That follows Levant et al. 2017 for GW170817 and we didn’t see a need to change it. Play with it for your transient and pick the best option for your particular cases.

# NOTE THAT THE X AND Y ARE REVERESED IN THE DATA COMPARED TO MATPLOTLIB!!!

cube2 = cube.subcube_circle_aperture(center=[119,49], # centroid y,x
                                     radius=2.5,  # this is for the OUTER BOUNDARY
                                     unit_center=None, # the center is given in Pixels so the Units are None
                                     unit_radius='arcsec' # the radii are in arcsec which we need to state here
                                    )

cube15 = cube.subcube_circle_aperture(center=[119,49],
                                      radius=1.5,  # this will be for the INNER BOUNDARY
                                      unit_center=None,
                                      unit_radius='arcsec')

f,ax = plt.subplots(nrows=1, ncols=1, figsize=(4,4))
ima = cube2.sum(axis=0)
maplot = ima.plot(scale='arcsinh', colorbar='v', ax=ax)
ax.set_title('2.5 Arcsec annulus')

Text(0.5, 1.0, '2.5 Arcsec annulus')

f,ax = plt.subplots(nrows=1, ncols=1, figsize=(4,4))
ima = cube15.sum(axis=0)
maplot = ima.plot(scale='arcsinh', colorbar='v', ax=ax)
ax.set_title('1.5 Arcsec annulus')

Text(0.5, 1.0, '1.5 Arcsec annulus')

Now all we need to do is integrate the SED in the 1.5 Arsec radius and subtract it from the integrated 2.5 arcsec flux and we’ll ahve recovered the flux from the pixels between the inner and the outer boundary.

Book keeping: We need a wavelength scale

Like in the voronoi binning notebook we create our 1D wavelength vector from the information we have in the headers.

# The first thing we need to do is create the wavelength scale. We've only got a starting
# wavelength, a Delta lambda and a number of bins in the spectrum (3861 - see above)
# The code below creates the wavelength array

# minimum and maximum wavelength from the header
wl_min = cube.data_header['CRVAL3']
wl_max = cube.data_header['CRVAL3']+cube.data_header['CD3_3']*cube. data.shape[0]

# creating the wavelength array
WL = np.arange(cube.data_header['CRVAL3'], wl_max, cube.data_header['CD3_3'])

Now we go fetch our data

In order to go get our data and extract the flux and the noise to integrate we need to do a bit of reshapeing.

# we need the shape information
shape_2  = cube2.data.data.shape
shape_15 = cube15.data.data.shape
shape_2, shape_15

((3681, 25, 25), (3681, 15, 15))

# then we reshape the data to be 2D
data_2arsec = cube2.data.data.reshape(shape_2[0], shape_2[1]*shape_2[2])
data_15arsec = cube15.data.data.reshape(shape_15[0], shape_15[1]*shape_15[2])

# axis = 0 is the spectrum; axis = 1 is the pixels all in a row
data_2arsec.shape, data_15arsec.shape

((3681, 625), (3681, 225))

# now we sum our spectra over all pixels and we subtract the 1.5 arcsec region from the 2.5 arcsec region
flux_annulus=data_2arsec.sum(axis=1) - data_15arsec.sum(axis=1)

# for the noise I did the same thing but in one line that's harder to look at
# take your pick on how you want to do that stuff in your own work
# also note I sqrt the variance to get the standard deviation
noise_annulus=np.sqrt(cube2.var.data.sum(axis=1).sum(axis=1) - cube15.var.data.sum(axis=1).sum(axis=1))

# Now we need to know what area on the sky (in pixels) we covered so we can
# turn our sum SED into an average SED
arrea_annulus=shape_2[1]*shape_2[2]-shape_15[1]*shape_15[2]

Let’s plot it!

plt.figure(figsize=(5,3))

plt.step(WL,
         flux_annulus/np.nanmedian(flux_annulus),
         zorder=10,
         c='k',
         where='mid'
        )

plt.errorbar(WL,
             flux_annulus/np.nanmedian(flux_annulus),
             yerr=noise_annulus/np.nanmedian(flux_annulus)
            )

plt.title("Annulus 1.5' to 2.5'")

Text(0.5, 1.0, "Annulus 1.5' to 2.5'")

Let’s store this in a neat dataframe.

data_annulus = pd.DataFrame(np.array([WL,
                            flux_annulus/np.nanmedian(flux_annulus),
                            noise_annulus/np.nanmedian(flux_annulus)]).T,
                            columns = ['WL', 'F', 'ERR'])

data_annulus.drop(data_annulus.tail(1).index,inplace=True) # because the last one was messed up

data_annulus

	WL	F	ERR
0	4749.617676	0.869223	0.018855
1	4750.867676	0.835849	0.012889
2	4752.117676	0.850857	0.012228
3	4753.367676	0.877018	0.011718
4	4754.617676	0.886446	0.011199
...	...	...	...
3675	9343.367676	0.624510	0.004571
3676	9344.617676	0.731427	0.005021
3677	9345.867676	0.832288	0.005965
3678	9347.117676	0.900837	0.007194
3679	9348.367676	0.907726	0.008474

3680 rows × 3 columns

1.2 Extracting the Transient data

We’re not done! We can now use this to extract our transient data and remove the galaxy background. This is actually not important for the galaxy SED fitting we’re doing in these tutorials but it is likely to be something that you will want to do when you are studying your shiny new transient.

Here we’re going to focus on the inner 1 arcsecond so a radius of 0.5 arcseconds. This value depends on your seeing when you’re observing as this is what will affect the full width half maximum.

# same as we did before but smaller
cube1 = cube.subcube_circle_aperture(center=[119,49], radius=0.5,  unit_center=None, unit_radius='arcsec')

plt.figure(figsize=(3,3))

ima = cube1.sum(axis=0)
ima.plot(scale='arcsinh', colorbar='v')

<matplotlib.image.AxesImage at 0x7f3e12b5e6a0>

ima.unmask() # we need to do this to access the data - because.

# same as earlier we do some reshaping to get all our ducks in a row, i mean pixels
shape_1  = cube1.data.data.shape
data_1arsec = cube1.data.data.reshape(shape_1[0], shape_1[1]*shape_1[2])
data_1arsec.shape

(3681, 25)

# let's see what they all look like!
plt.figure(figsize=(5,3))

for spec in data_1arsec.T:
    plt.plot(WL, spec, c='k', alpha=0.1)

plt.xlim([4800,7500])

(4800.0, 7500.0)

# let's sum our caterpillars
transient_flux = data_1arsec.sum(axis=1) # spectrum
noise_transient=np.sqrt(cube1.var.data.sum(axis=1).sum(axis=1)) # standard deviation

Let’s plot the transient (uncorrected)

plt.figure(figsize=(10,5))
plt.plot(WL, transient_flux/np.nanmedian(transient_flux), zorder=12)
plt.errorbar(WL, transient_flux/np.nanmedian(transient_flux), yerr = noise_transient/np.nanmedian(transient_flux))
plt.xlim([4800,7500])
plt.ylim([0.8,1.2])
plt.title('Within 0.5 arsecond of the transient')

Text(0.5, 1.0, 'Within 0.5 arsecond of the transient')

Now we need to remove the galaxy background

To do this we need to first figure out how much “weight” the galactic SED has in our transient SED given the number of pixels we have summed together.

The weight is the fraction of the light in the transient spectra that comes from the galaxy based on the surface brightness of the neighbours in the annulus. We need to compare the areas of we’ve been looking at so far:

cube2.shape, cube15.shape, cube1.shape

((3681, 25, 25), (3681, 15, 15), (3681, 5, 5))

area_cube1 = (cube1.shape[1]*cube1.shape[2])
area_cube2 = (cube2.shape[1]*cube2.shape[2])
area_cube15= (cube15.shape[1]*cube15.shape[2])

# fraction we need to multiply to our galactic SED
weight_gal_spec = area_cube1/(area_cube2-area_cube15)

# correcting the transient flux by removing (subtracting) the weighted galactic SED
flux_gal_removed = (transient_flux/np.nanmedian(transient_flux))[:-1] - data_annulus.F*weight_gal_spec

Let’s plot the whole thing!

as it so happens here it didn’t make much of a difference to our transient but it might have! Also note the scales between the transient and the galactic neighrbourhood because everything is normalised to be around 1 (which we need for fitting later).

f, ax = plt.subplots(nrows=2, figsize=(10,7))


ax[0].plot(WL[:-1],flux_gal_removed, label='Galaxy removed')
ax[0].plot(WL, transient_flux/np.nanmedian(transient_flux)-0.12, label='raw - 0.12')
ax[0].set_xlim([4800,9000])
ax[0].set_ylim([0.7,1.2])
ax[0].legend()
ax[0].set_title('AT 2017gfo')


ax[1].plot(WL[:-1],data_annulus.F, label='Galaxy aroudn transient', c='k')
ax[1].set_xlim([4800,9000])
ax[1].set_ylim([0.6,1.2])
ax[1].set_title('NGC 4993 around AT 2017gfo')
ax[1].set_xlabel(r'Wavelength ($\AA$)')

Text(0.5, 0, 'Wavelength ($\\AA$)')

np.nanmedian(flux_annulus)

27978.792442321777

Let’s save our extracted spectra to a data frame

And we can then pickle the data frame (there are other ways to save dataframe data I just used to pickle them and I’ve left it as such here - make your prefered choice)

data = pd.DataFrame(columns=['WL', 'transient', 't_err', 'transient_raw', 'neighbours', 'n_err'])
data.WL = WL[:-1]
data.transient = flux_gal_removed
data.transient_raw = transient_flux[:-1]/np.nanmedian(transient_flux)
data.t_err=noise_transient[:-1]/np.nanmedian(transient_flux)
data.neighbours=data_annulus.F
data.n_err=data_annulus.ERR

data.metadata="WL: Wavelength in angstrom // Transient: Normalised transiant flux GALAXY REMOVED // Trasnient_raw: Normalised transient flux - galaxy NOT removed // neighbours: galaxy flux in anulus 1.5-2 arcsec around transient"

data.to_pickle("data/transient_and_neighbour_flux.pkl")

2 - Fitting the neighbourhood SED with BPASS!

Now that we have extracted our host environment spectral information let’s fit it to see what it is made of.

data.head()

	WL	transient	t_err	transient_raw	neighbours	n_err
0	4749.617676	0.868068	0.025423	0.922394	0.869223	0.018855
1	4750.867676	0.821615	0.015944	0.873856	0.835849	0.012889
2	4752.117676	0.835267	0.015323	0.888445	0.850857	0.012228
3	4753.367676	0.851539	0.014860	0.906352	0.877018	0.011718
4	4754.617676	0.834517	0.014208	0.889920	0.886446	0.011199

2.1 - Create ppxf templates

I recommend checking the ppxf manual for extra detail but I am going to walk you through the essentials

# just to make some of our code less long we'll store some columns in short variable names
WL = data.WL.values
F = data.neighbours.values
noise = data.n_err.values

# Now we need to de-redshift
z = 0.009783
recessional_vel = z*c
WL = (WL - WL*z)

2.1.1 - BPASS models

To make templates we’ll need the BPASS models. To download the models you can check the BPASSv2.2.1 starter kit on Zenodo. Then we’ll need to be able to locate them so you’ll need to change the code below

BPASS_MODEL_PATH = '../../BPASS_hoki_dev/bpass_v2.2.1_imf135_300/'

2.1.3 - Log (natural not log10) rebin the wavelength

This is required by the ppxf magic. Feel free to check out Capellari 2003 and 2017 if you want to read more on how ppxf works and why the log rebin is a really neat mathematical trick ;)

# ppxf.util gives us the log_rebin function
flux, wl_log_bin, velscale = util.log_rebin([WL[0], WL[-1]], F)

# later we'll need the linear scale of our log rebinned spectra to plot our fits
# this step ensures we have a linear wavelength scale for those fits with the
# right number (and location) of wavelenght bins!
wl_fits = np.exp(wl_log_bin)

2.1.4 - Creating templates for ppxf with KVN

KVN is a helper class in hoki.sedfitting that does the dirty work of creating templates that can be fed to ppxf - it also has utilities to decipher what ppxf has found.

hsed.KVN?

kvn=hsed.KVN()
kvn.make_templates(BPASS_MODEL_PATH,
                   fwhm_obs=1.25,
                   wl_obs=WL, # Wavelength bins of the observations (not log rebinned)
                   wl_range_obs=[WL[0], WL[-1]], # wavelength range of the observations
                   velscale_obs=velscale,
                   wl_range_padding=[-50,50], # number of extra Angstroms befrore and after the range
                   z_list=['z010','z020','z030'],
                   verbose=True,
                   )

[---INFO---] TemplateMaker Starting
[--RUNNING-] Initial Checks
[--RUNNING-] Loading  model spectrum
[--RUNNING-] Calulating obs. velocity scale -- No dispersion
[--RUNNING-] Calculating template wavelength (log rebin) and velocity scale
[--RUNNING-] Calculating sigma
[--RUNNING-] Instanciating container arrays
[---INFO---] Using all ages from 6.0 to 10.2 (included) - log_age_cols now set
[--RUNNING-] Compiling your templates
Progress: |██████████████████████████████████████████████████| 100.00% Complete
[-COMPLETE-] Templates compiled successfully

Padding: When creating the templates some rebinning has to occur and if we only load the exact wavelength range this can cause problems with the templates and the fitting later. Then you might ask why we crop our BPASS SEDs at all then, and the answer is simple: size! The BPASS SEDs cover everything from 1 A to 100,000 A. We’d be creating unnecessary work for our computer so it’s better to crop.

Where and how are the templates stored

The SED templates are BPASS SEDs that have been appropriately processed to be used by ppxf for use on our specific data set. Unfortunately processing needs to happen on the templates for different observations with different properties like resolution element and wavelength range. That is why the helper class is so necessary and I can’t just give you a library of template to slap into your fit.

Let’s look at where to find our newly created templates:

kvn.templates.shape

(3766, 42, 3)

The dimensions are [Wavelength, age, metallicity]. But obviously age = 42 and metallicity = 1 is not useful physically, we need to know what that corresponds to. The template properties are stored in the attribute template_properties which has 2D [N templates, properties] where the properties are numpy arrays containing [metallicity, age].

kvn.template_properties[0]

array([0.01, 6.  ])

Tempalte 0 has a metallicity Z=0.010 and an age of 10^6 or 1 million years.

Note that here metallicity comes before the age dimension (in practice you should never need to think about it). The order the templates come in and the order the template properties are arranged is according to what is needed to talk to ppxf: feed it the templates in the order it requries and reading its solutions in the order it provides.

This ordering thing is why we have a helper class so we can abstract it away.

2.2 - Seting up fit with PPXF

2.2.1 - Start guesses for the kinematics

Below we define two start guesses for the recession velocity and the dispersion velocity of the galaxy. ppxf will fit the kinematics at the same time as looking for a best SED template combination so we need to give it a fighting chance.

Here we have already corrected for the average recession velocity so our initial guess is 0. It’s likely the final fitted guess will differ since galaxies rotated and our transient is on the outskirts so likely spining towards or away from us.

In terms of dispersion velocity you can make guesses using breadth of lines you know to be narrow, or you can eyeball it given the galaxy type and then tweak the number if the fits come out garbage. It doens’t need to be spot on to be decent usually.

start = [0, 180] #180km/s from previous fits

2.2.2- Good pixels

In ppxf you can provide information on which pixels to take into account and which to ignore. Let’s see an example here (some more artificial than others but here for the pedagogy):

CUTTING BASED ON FLUX

1. Imagine if we have an emission line we don’t want to fit (NOTE: ppxf has some stuff to take into account emission features which you can read into in their manual). A stupid but often effective way to handle this is to crop anything that rises to high! Here we don’t have have any to worry about but when fitting all voronoi bins we do because the center of NGC 4993 is a LIER region (see Levan et al. 2017). Say we want to crop everything that rises above 2.5 (determined through visual anlysis in other spectra)

2. Our flux should be positive - if it’s not something is probably messed up so let’s cut everything that drops below 0.

CUTTING BASED ON WAVELENGTH

3. Let’s say we want to hide specific line because we thing our templates won’t be able to fit it well; here we pick the first absorption region; we can also do this by telling ppxf the pixels between that range are not good.

# 1
FLUXMAX =2.5
# 2
FLUXMIN = 0.0
# 3
WLMIN = 4875
WLMAX = 4975

# This looks convoluted but this is the format ppxf wants
# A list of which wavelenght pixels are good (THEIR INDEX not their value)

goodpixels = np.argwhere((flux < FLUXMAX) & (flux > FLUXMIN) # Flux condition
                         & ((WL < WLMIN) | (WL > WLMAX)) # Wavelength condition
                         ).T[0] # Make it a 1D horizontal vector

2.2.3 - Fitting with ppxf

The method was originally described in 2004 but the most extensive paper on ppxf is Capellari 2017. When you download ppxf there are also exmaple scripts provided with the code and what I use is adapted from those scripts. You can also consult the online PyPi manual.

A few quick start things to know:

• moments: if set to 2, will only fit the kinematic parameters [recession velocity, dispersion]. We set it to 4 because we are also fitting for SFH and metallcity

• degree: adding a polynomial can help if your spectrum is not super well calibrated or if there is a little bit of reddening are not fitting for with the reddening utilities (which I’m not here because there is not a lot of dust and I don’t care about the nature of the dust).

• clean: this does sigma clipping and can help or hinder. See the Supplementary Information in Stevance et al. 2023 (UPDATE) to see how that played a role in our fitting of NGC 4993.

# From the ppxf manual:
# Reference velocity in km/s (default 0). The input initial guess and the output velocities
# are measured with respect to this velocity. This keyword can be used to account for
# the difference in the starting wavelength of the templates and the galaxy spectrum.
dv = c*np.log(kvn.wl_tem[0]/kvn.wl_obs[0])

### FEEDING PPXF AND FITTING
pp = ppxf(kvn.templates,           # these are the BPASS SED templates
          flux,                    # observed flux vector
          noise,                   # the standard deviation vector
          velscale,                # the velocity scale calculated by the log rebin ppxf utility
          start,                   # the start guesses for recession vel and dispersion
          goodpixels=goodpixels,   # our goodpixels
          plot=False,              # I set this to false because I made my own in hoki.sedfitting
          moments=4,               # order of the Gauss-Hermit moment to fit (see above)
          degree=2,                # order of the polynomial to fit on top of the SED mixture
          vsyst=dv,                # see above
          clean=False,             # whether to perform sigma clipping and itterative fitting
          lam=WL                   # the observed wavelength vector
         )

plt.ylim([0.4,1.2])


# FITTING THE RESULTS!
res = hsed.ppxf_plot_fixed(pp,
                           wl_fits=wl_fits,
                           WL=WL, F=F)

 Best Fit:       Vel     sigma        h3        h4
 comp.  0:       -36       164     0.005     0.064
chi2/DOF: 23.05; degree = 2; mdegree = 0
method = capfit; Jac calls: 2; Func calls: 12; Status: 2
linear_method = lsq_box; Nonzero Templates (>0.1%):  5 / 126

2.3 - Saving and exploring the results with KVN

2.3.1 ``make_results`` and SED weights

KVN is not only a template creator, it also deciphers the results from ppxf to give them to us in a neat table and (as we’ll see later) we can save KVN as a pickle file so that our templates and fitting results are stored in one place and can be called back into other codes or shared with colleagues without needed to rerun the fitting.

First we need to ‘make the results’ by putting the ppxf.ppxf.ppxf object we created (here called pp) during fitting into our KVN

kvn.make_results(pp)

kvn.results

	met	age	weights
0	0.02	8.9	0.039217
1	0.03	9.0	0.077932
2	0.01	9.5	0.282404
3	0.01	9.7	0.316477
4	0.01	9.9	0.094677

The results are stored in a pandas.DataFrame. * met is the metallicity (absolute, not with respect to solar Z=0.010 is 1% metals) * age is the log10 age. That’s the way ages are binned in BPASS there is not other option * weights are the weights ppxf gives those respective SEDs. They don’t necessarily add up to 1. If you want percentages or fractions you’ll need to rescale.

The only SED properties shown in the table are the SEDs ppxf actually uses for the fit.

2.3.2 - The matching_ spectra and polynomials

When we make the results KVN also goes to fetch the SEDs used to match the observed SED so we can have a convenient access to them: * kvn.matching_raw_spectra: the raw templates pre-fit (NOT THE RAW BPASS SEDs - they are processed for ppxf) * kvn.matching_spectra: the spectra processed with the kinematics found by ppxf * kvn.matching_apolynomial: the additive polynomial (if you included one) * kvn.matching_mpolynomial: the multiplicative polynomial (if you included one)

MAJOR NOTE: The kvn.matching_spectra need to be plotted with wl_fits NOT with kvn.wl_tem or kvn.wl_obs and it can’t be either any of the logarithming wavelength scales we have flowting around otherwise it will have weird offsets compared to NGC4993 (or your galaxy). There are a lot of different wavelength scales we deal with in this process, remembering which goes with what is one of the main beginner pit falls so if your spectra don’t line up this should be your first port of call - what wavelength scale am I using and is it the right one?

CORRESPONDING WAVELENGTH SCALE FOR PLOTTING

WL scale	Plotted spectra
kvn.wl_tem	kvn.matching_raw_spectra
wl_fits	kvn.matching_spectra
wl_fits	kvn.matching_apolynomial
wl_fits	kvn.matching_mpolynomial

Let’s have a look at what these various templates look like. We’ll also plot the fit.

# TEMPLATE (RAW POST-FIT)
plt.plot(kvn.wl_tem,
         kvn.matching_raw_spectra[0],
         color='orange', alpha=0.8,
         label='A RAW template'
        )

# TEMPLATE (WITH KINEMATICS)
plt.plot(wl_fits,
         kvn.matching_spectra[0],
         color='red', lw=1.5,
         label='A template kinematically processed'
        )

# TEMPLATE (WITH KINEMATICS)
plt.plot(wl_fits,
         kvn.matching_apolynomial,
         color='maroon', lw=1.5,
         label='The added polynomial'
        )

# OBSERVED NGC 4993
plt.plot(wl_fits,
         flux,
         color='k', lw=1, alpha=0.9,
         label='NGC 4993'
        )

# FIT TO NGC 4993
plt.plot(wl_fits,
         kvn.ppxf.bestfit,
         color='blue', lw=0.9,
         label='The fit'
        )

plt.legend(fontsize=8, ncol=2)

<matplotlib.legend.Legend at 0x7f3e0aada8e0>

Note that to get the full fit that has put together the SEDs with the right weights and included the kinematics and the polynomial we call kvn.ppxf.bestfit. Here we are directly looking into the ppxf object which contains all the results. The kinematically process tempaltes and polynomial can also be found in there alongside a lot more info (consult the manual for full details) such as the chi2 and a record of all the parameters we gave like the velscale.

Now a final point we need to address is What’s the age and metallicity of that template we just plotted?. I don’t know that it will necessarily match the order in which things are reported in the results table so you need a way to link your matching templates to their template properties - for that we have kvn.match_indices. As the name suggestes they are the indices identifying the template spectra in our full long list of teemplate spectra.

kvn.match_indices

array([ 88,  92, 105, 111, 117])

So the properties are:

kvn.template_properties[kvn.match_indices]

array([[0.02, 8.9 ],
       [0.03, 9.  ],
       [0.01, 9.5 ],
       [0.01, 9.7 ],
       [0.01, 9.9 ]])

Which matches what’s in the result table at the minute, but if you had processed that table and re-ordered things (which you might do in some circumstances) you’d be in trouble so this is the most sure way to get your matching template properties.

2.3.3 Pickling KVN

The final thing we’re going to do is pickle our KVN object so we can save our work and then load it again next time we want to plot our fit or look at the star formation weights found by ppxf. This way we have everything in one place: the templates, the full fit solutions (with the ppxf object contained within KVN), and some convenient summaries of our results like the results table.

kvn.save('data/pwiz_at2017gfo.pkl')

To load the pickle file we need create a KVN object and then load

pwiz2 = hsed.KVN()
pwiz2.load('data/pwiz_at2017gfo.pkl')

# just to show it's the same
sum(pwiz2.ppxf.bestfit-kvn.ppxf.bestfit)

0.0

BONUS - Example of a bad fit

A common source of bad fit with ppxf is something wrong with the wavelenghts or velocities, either the templates weren’t created right because a log wavelength was given where a linear wavelength vector was expected (or vice versa),or the dv are wrong.

The result will be a spectrum that CLEARLY doesn’t have the right kinematics: it’s too broad. Let’s see an example below so you can identify that problem quickly and double check your steps.

pp = ppxf(kvn.templates,           # these are the BPASS SED templates
          flux,                    # observed flux vector
          noise,                   # the standard deviation vector
          velscale,                # the velocity scale calculated by the log rebin ppxf utility
          start,                   # the start guesses for recession vel and dispersion
          goodpixels=goodpixels,   # our goodpixels
          plot=False,              # I set this to false because I made my own in hoki.sedfitting
          moments=4,               # order of the Gauss-Hermit moment to fit (see above)
          degree=2,                # order of the polynomial to fit on top of the SED mixture
          vsyst=0,                 # CHANGED THIS TO 0
          clean=False,             # whether to perform sigma clipping and itterative fitting
          lam=WL                   # the observed wavelength vector
         )

plt.ylim([0.4,1.2])


# PLOTTING THE RESULTS!
hsed.ppxf_plot_fixed(pp, wl_fits=wl_fits, WL=WL, F=F)

 Best Fit:       Vel     sigma        h3        h4
 comp.  0:     -2000      1000    -0.300    -0.081
chi2/DOF: 45.83; degree = 2; mdegree = 0
method = capfit; Jac calls: 10; Func calls: 52; Status: 2
linear_method = lsq_box; Nonzero Templates (>0.1%):  3 / 126

array([-0.11513668, -0.14540687, -0.12683344, ..., -0.11569206,
       -0.04607372, -0.03827896])