Continuum cookbook

This cookbook contains utilities to detect absorption features in the spectrum and remove them, to determine the continuum level before absorption.

Table of contents

1. Clip flux in spectrum
2. Find lines
3. Continuum from nodes
4. Update line list
5. Find peaks
6. Extract nodes
7. Clean nodes
8. Interpolate nodes
9. Correct spectrum for Lyman-alpha opacity

Clip flux in spectrum

{
  "cookbook": "cb",
  "recipe": "flux_clip",
  "params": {
    "hwindow": "100",
    "kappa": "5",
    "iter": "100",
    "std": "500"
  }
}    

Discriminate absorbed spectrum bins by applying a kappa-sigma clipping within a running window.

The recipes computes the mean of y in a runniung window across the spectrum and saves it in column y_rm. It then computes the deviation \(\Delta y =\) y \(-\) y_rm for all spectral bins and finds the outliers in this distribution.

Since the recipe is meant to operate on quasar spectra, it assumes that the absorption features are narrow and widespread, while the emission feature are wide enough to be considered a part of the continuum. For this reason, the recipe is biased towards detecting outliers with a negative deviation (lower outliers, i.e. absorbed bins) instead of a positive deviation (upper outliers, most likely residuals of sky line or cosmic ray hit subtraction).

In practice, a bin is regarded as an outlier if it has either \(\Delta y < -\)k \(\times\) dy or \(\Delta y > 2\)k \(\times\) dy. Notice the two in the second formula, meaning that a bin must deviate twice as much in the upper direction than in the lower direction to be considered as an outlier.

The clipping is iterated several times, always removing the outliers and re-computing y_rm, until no more outliers are found or a maximum number of iterations iter is reached.

A column y_abs is created and set equal to 1 for all the lower outliers (absorbed bins), 0 elsewhere. Similarly, a column y_em is created for the upper outliers (“emitted” bins) and a column y_cont is created for the remaining bins (continuum bins). The latter is smoothed in the velocity space with a gaussian filter with standard deviation std, and the result is saved in column cont of the spectrum as the current estimate of the emission continuum.

Find lines

{
  "cookbook": "cb",
  "recipe": "lines_find",
  "params": {
    "std_start": "100.0",
    "std_end": "0.0",
    "col": "'y'",
    "kind": "'min'",
    "kappa_peaks": "5.0",
    "resol": "null",
    "append": "true"
  }
}    

Create a line list by convolving a spectrum with different gaussian profiles and finding the peaks in the convolution.

The recipe is an extension of peaks_find, using convolution to detect peaks of different width. The gaussian profiles are applied in velocity space and used to filter out all fluctuations below a given scale.

The recipe proves most effective when a wide range of line widths is spanned, from the widest to the narrowest (typically, with gaussian standard deviation ranging from 100 to 0 km/s). A growing list of detected peaks is created and updated while spanning the range of line widths, eliminating duplicated findings.

The final list of x, xmin, xmax, y, and dy values for the peaks are saved in the current line list. If append is True, the new peaks are appended to the existing line list (if present), otherwise they replace it.

N.B. If the same line is found more than once, using different line widths for the gaussian filtering, it is saved just once in the line lists, and the value of xmin and xmax are the positions of the closest extrema (the closes maxima if the peak is a minimum, and vice-versa) when the line was first detected. This is the reason for spanning line widths in decreasing order: if a wide line is detected at first with a small value of gaussian standard deviation, it may be assigned a xmin-xmax range much smaller than the region it actually covers (because the filter is too fine to smooth out the noise around the peak), and the range is preserved even if the line is detected again with wider gaussians.

Continuum from nodes

{
  "cookbook": "cb",
  "recipe": "nodes_cont",
  "params": {
    "delta_x": "500",
    "kappa_nodes": "5.0",
    "smooth": "0"
  }
}    

Estimate a continuum by extracting, cleaning, and interpolating nodes from regions not affected by lines.

This recipe combines nodes_extract, nodes_clean, and nodes_interp to estimate the continuum level once the lines have been detected.

Update line list

{
  "cookbook": "cb",
  "recipe": "lines_update",
  "params": {
  }
}    

Update the line list with the parameters obtained by fitting the absorption systems.

This recipe is meant to be run after the absorption systems have been fitted with Voigt profiles. It does not work if the list of systems comes from a different line lists.

For each systems, the best-fitting Voigt parameters are propagated to the line list and the FWHM of the lines is computed.

Find peaks

{
  "cookbook": "cb",
  "recipe": "peaks_find",
  "params": {
    "col": "'conv'",
    "kind": "'min'",
    "kappa": "5.0",
    "append": "true"
  }
}    

Find the peaks in a spectrum column.

Peaks are defined as the extrema (minima or maxima) of y column that are more prominent than kappa times the value of dy with respect to neighbouring bins.

The recipe finds all the peaks in the spectrum and saves their x, xmin, xmax, y, and dy values in the current line list. If the peak is a minimum, xmin and xmax are determined as the position of the closest maxima, and vice-versa if it is a maximum. If append is True, the new peaks are appended to the existing line list (if present), otherwise they replace it.

A column line_mask is created in the spectrum (if not present) and set to 1 within the xmin-xmax range of each detected peak.

Extract nodes

{
  "cookbook": "cb",
  "recipe": "nodes_extract",
  "params": {
    "delta_x": "500",
    "xunit": "km / s",
    "mode": "'std'"
  }
}    

Extract nodes from a spectrum, optionally masking lines.

Depending on mode, the recipe can be used either to create a new set of nodes to estimate the continuum (std) or to place nodes across an existing continuum profile (cont). In the latter case, the spectrum must have a cont column beforehand.

In std mode, the recipe splits the spectrum in a number of slices and extracts a node for each slice by averaging the values of x, y, and dy within it. The node list is formatted like a normal spectrum.

Slices are defined in either wavelength or velocity space, as specified by the chosen xunit. They are equally spaced, with their size determined by delta_x.

If the spectrum has a lines_mask column, only the regions where this mask is 0 are considered when etracting the position of the nodes. This means that a slice may not be assigned a node if it is completely masked.

If the spectrum has a deabs column, containing the de-absorbed flux computed by fitting and removing the absorption features, this column is used instead of y to determine the y position of the nodes.

In cont mode, the cont column is used instead of y to determine the y position of the nodes, and the lines_mask is ignored.

Clean nodes

{
  "cookbook": "cb",
  "recipe": "nodes_clean",
  "params": {
    "kappa": "5.0"
  }
}    

Remove outliers from a list of nodes.

The recipe applies the peak finding algorithm of peaks_find to the list of nodes, to detect outliers and remove them iterately. Outliers are defined as the peaks that are more prominent than kappa times the value of dy with respect to the neighbouring nodes.

Interpolate nodes

{
  "cookbook": "cb",
  "recipe": "nodes_interp",
  "params": {
    "smooth": "0"
  }
}    

Interpolate nodes to estimate the emission level.

The nodes in the current nodes list are interpolated with scipy.interpolate.UnivariateSpline. The interpolating function is sampled at all values of x of the spectrum and the line list (if present), and is saved as the current continuum level in column cont.

Correct spectrum for Lyman-alpha opacity

{
  "cookbook": "cb",
  "recipe": "lya_corr",
  "params": {
    "zem": "zem",
    "input_col": "'y'",
    "mode": "'basic'",
    "logN_thres": "100",
    "percentile": "100"
  }
}    

Correct the spectrum flux for the effective Lyman-alpha opacity.

Depending on mode, the correction factor is computed either from a simple integration of the H i growth function (basic), or using the prescriptions by Inoue et al. (2014) (inoue).

The basic mode require an upper integration limit that can be either provided through logN_thres or extracted from the distribution of H i column densities, from the table of fitted absorption systems (if present). In the latter case, the column density corresponding to a given percentile of the distribution is used as the upper integration limit.

The correction factor is multiplied to column input_col of the spectrum and the line list (if present), and the result is saved in column [input_col]_taucorr.