Fitting SEDs to sources¶
Typical workflow¶
Once the data has been prepared following the description in Data format, and that you have a set of models in the format described in Creating model packages the typical workflow to follow when using the SED fitting tool is:
Fit a given set of models to the data
Optionally split the output between well- and badly-fit sources
Make plots of the SED fits, or of the parameter distribution
Each of these steps can be performed by one ore more functions in the SED Fitting tool, which are described below in detail:
Note on quantities and units¶
The Python SED fitting tool makes use of the astropy.units
package to handle
units and unit conversions. Several of the options that need to be specified in
the functions described below require Quantity
instances. Defining quantities is straightforward:
from astropy import units as u
# Define scalar quantity
q1 = 3. * u.kpc
# Define array quantity using a list
q2 = [1., 2., 3.] * u.arcsec
# Define array quantity using a Numpy array
q3 = np.array([1., 2., 3.]) * u.cm ** 2 / u.g
Fitting the models to the data¶
In order to fit the data, you will need to make use of the
fit()
function, which can be imported using:
from sedfitter import fit
At a minimum, you will need to call the function with the following input:
The file containing the data (see Data format).
The list of filters and aperture radii for which the data is given (see Common filter names).
The directory containing the models (see Available Model Packages and Creating model packages).
The name of the output file
The range of Av and distances to assume for the fits
An extinction law, given as a opacities to extinction (in cm^2/g or equivalent) versus wavelength.
The fit()
page describes in detail how each of these inputs
should be specified, and also lists additional options to further fine-tune the
fitting and the output.
The following is a complete example showing how to use the
fit()
function:
from astropy import units as u
from sedfitter import fit
from sedfitter.extinction import Extinction
# Define path to models
model_dir = '/Volumes/Data/models/models_r06'
# Read in extinction law)
extinction = Extinction.from_file('kmh94.par', columns=[0, 3],
wav_unit=u.micron, chi_unit=u.cm**2 / u.g)
# Define filters and apertures
filters = ['2J', '2H', '2K', 'I1', 'I2', 'I3', 'I4']
apertures = [3., 3., 3., 3., 3., 3., 3.] * u.arcsec
# Run the fitting
fit('data_glimpse', filters, apertures, model_dir,
'output.fitinfo',
extinction_law=extinction,
distance_range=[1., 2.] * u.kpc,
av_range=[0., 40.])
Note
in the filter list, you can also specify wavelengths as Astropy
Quantity
instances. If you do this, the SED
wavelength closest to that specified will be used in the fitting.
Note
if you do not specify the columns and units when reading in the
extinction, the first two columns are read and are assumed to be in
c.g.s.. If you have previously used the Fortran version of the SED
fitter, you will need to specify columns=[0, 3]
to choose the
first and fourth column.
Plotting SEDs¶
Once you have fit the data, you will likely want to plot the resulting SED
fits. To do this, you will need to make use of the plot()
function, which can be imported with:
from sedfitter import plot
The plot()
requires the output file from the
fit()
function as well as the name of an output directory. For
example, continuing the example above, you can do:
from sedfitter import plot
plot('output.fitinfo', 'plots_seds')
By default, only the best-fit parameter is shown, but this can be changed by
using the select_format
option, which is described in more detail in
Model selection syntax. For example, to write out all the models with a
\(\Delta\chi^2\) value per data point (relative to the best fit) of less
than 3, you can do:
plot('output.fitinfo', 'plots_seds', select_format=('F', 3))
In addition, there are many options available to
customize the format and appearance of the plots. For more information about
these options, see the plot()
page.
Plotting parameters¶
Functions are available to make 1- and 2-d parameter plots:
from sedfitter import plot_params_1d, plot_params_2d
As when Plotting SEDs, one needs to specify the output file from the
fit()
function, the output directory, and the name of the
parameters to plot:
from sedfitter import plot_params_1d, plot_params_2d
# Make histograms of the disk mass
plot_params_1d('output.fitinfo', 'MDISK', 'plots_mdisk',
log_x=True)
# Make 2-d plots of the envelope infall rate vs disk mass
plot_params_2d('output.fitinfo', 'MDISK', 'MDOT', 'plots_mdot_mdisk',
log_x=True, log_y=True)
By default, only the best-fit parameter is shown, but this can be changed by
using the select_format
option, which is described in more detail in
Model selection syntax. In addition, there are many options available to
customize the format and appearance of the plots. For more information about
these options, see the plot_params_1d()
and
plot_params_2d()
pages.
Splitting well- and badly-fit sources¶
After computing the fits with fit()
, it is possible to split
the output file on the basis of the \(\chi^2\) value for the best-fit. This
is done using the filter_output()
function which is imported
with:
from sedfitter import filter_output
For example, to split the above file into well- and badly-fit sources based on the absolute \(\chi^2\) of the best-fit, you can do:
filter_output('output.fitinfo', chi=3.)
This will produce files named output_good.fitinfo
and
output_bad.fitinfo
by default (although you can also specify custom
names for the output files). It is also possible to split the fits based on the
\(\chi^2\) value per datapoint using the cpd
option. More information
about the available options is available in filter_output()
.
Extracting the fit and model parameters¶
The output files produced above are in binary format and are not
human-readable. To produce ASCII files of the output, you can use the
write_parameters()
and
write_parameter_ranges()
functions. The former is used to
write out all the parameters of all the models requested, while the latter will
only write out the minimum and maximum for each parameter. The functions are imported with:
from sedfitter import write_parameters, write_parameter_ranges
To use these functions, you will need to specify the input binary file, and the output ASCII file name:
from sedfitter import write_parameters, write_parameter_ranges
# Write out all models with a delta chi^2-chi_best^2 per datapoint < 3
write_parameters('output.fitinfo', 'parameters.txt',
select_format=('F', 3.))
# Write out the min/max ranges corresponding to the above file
write_parameter_ranges('output.fitinfo', 'parameter_ranges.txt',
select_format=('F', 3.))
More information about the available options is given in
write_parameters()
and
write_parameter_ranges()
.