Fitting an Absorption Spectrum¶

Loading an Absorption Spectrum¶

To load an absorption spectrum created by (:ref: AbsorptionSpectrum), we specify the output file name. It is advisable to use either an .h5 or .fits file, rather than an ascii file to save the spectrum as rounding errors produced in saving to a ascii file will negatively impact fit quality.

f = h5py.File('spectrum.h5')
wavelength = f["wavelength"][:]
flux = f['flux'][:]
f.close()

Specifying Species Properties¶

Before fitting a spectrum, you must specify the properties of all the species included when generating the spectrum.

The physical properties needed for each species are the rest wavelength, f-value, gamma value, and atomic mass. These will be the same values as used to generate the initial absorption spectrum. These values are given in list form as some species generate multiple lines (as in the OVI doublet). The number of lines is also specified on its own.

To fine tune the fitting procedure and give results in a minimal number of optimizing steps, we specify expected maximum and minimum values for the column density, doppler parameter, and redshift. These values can be well outside the range of expected values for a typical line and are mostly to prevent the algorithm from fitting to negative values or becoming numerically unstable.

Common initial guesses for doppler parameter and column density should also be given. These values will not affect the specific values generated by the fitting algorithm, provided they are in a reasonably appropriate range (ie: within the range given by the max and min values for the parameter).

For a spectrum containing both the H Lyman-alpha line and the OVI doublet, we set up a fit as shown below.

HI_parameters = {'name':'HI',
        'f': [.4164],
        'Gamma':[6.265E8],
        'wavelength':[1215.67],
        'numLines':1,
        'maxN': 1E22, 'minN':1E11,
        'maxb': 300, 'minb':1,
        'maxz': 6, 'minz':0,
        'init_b':30,
        'init_N':1E14}

OVI_parameters = {'name':'OVI',
        'f':[.1325,.06580],
        'Gamma':[4.148E8,4.076E8],
        'wavelength':[1031.9261,1037.6167],
        'numLines':2,
        'maxN':1E17,'minN':1E11,
        'maxb':300, 'minb':1,
        'maxz':6, 'minz':0,
        'init_b':20,
        'init_N':1E12}

speciesDicts = {'HI':HI_parameters,'OVI':OVI_parameters}

Generating Fit of Spectrum¶

After loading a spectrum and specifying the properties of the species used to generate the spectrum, an apporpriate fit can be generated.

orderFits = ['OVI','HI']

fitted_lines, fitted_flux = generate_total_fit(wavelength,
    flux, orderFits, speciesDicts)

The orderFits variable is used to determine in what order the species should be fitted. This may affect the results of the resulting fit, as lines may be fit as an incorrect species. For best results, it is recommended to fit species the generate multiple lines first, as a fit will only be accepted if all of the lines are fit appropriately using a single set of parameters. At the moment no cross correlation between lines of different species is performed.

The parameters of the lines that are needed to fit the spectrum are contained in the fitted_lines variable. Each species given in orderFits will be a key in the fitted_lines dictionary. The entry for each species key will be another dictionary containing entries for ‘N’,’b’,’z’, and ‘group#’ which are the column density, doppler parameter, redshift, and associate line complex respectively. The i ^th line of a given species is then given by the parameters N[i], b[i], and z[i] and is part of the same complex (and was fitted at the same time) as all lines with the same group number as group#[i].

The fitted_flux is an ndarray of the same size as flux and wavelength that contains the cummulative absorption spectrum generated by the lines contained in fitted_lines.

Saving a Spectrum Fit¶

Saving the results of a fitted spectrum for further analysis can be accomplished as follows using the h5 file format.

f = h5py.File('fitted_lines.h5')

for species in orderFits:

    f.create_dataset("{0}/N".format(species),
        data=fitted_lines[species]['N'])
    f.create_dataset("{0}/b".format(species),
        data=fitted_lines[species]['b'])
    f.create_dataset("{0}/z".format(species),
        data=fitted_lines[species]['z'])
    f.create_dataset("{0}/complex".format(species),
        data=fitted_lines[species]['group#'])

f.close()