PyWiFeS: A Rapid Data Reduction Pipeline
for the Wide Field Spectrograph (WiFeS)

The first step in CCD data reduction is to extract the science pixels from the raw data and convert that data from ADUs to real photon counts (i.e., electrons). This is done by measuring and subtracting the overscan level from the overscan regions of the data, then multiplying the overscan-subtracted ADU values by the gain of the readout electronics. This is done seamlessly for all epochs of WiFeS data with the subtract_overscan routine, using detector characteristics stored in the WiFeS metadata file.

The second step of image pre-processing is to remove any cosmetic defects in the detector. The first generation of detectors (Fairchild CCDs) in WiFeS were free of any cosmetic defects. In 2013 a second generation of detectors (E2V detectors) were installed in March and May of 2013 for the red and blue detectors, respectively. These detectors provided much improved performance in terms of long-term stability and read noise, but each have a few columns (1 for red, 2 for blue, out of over 4000) of the CCD which are defective in most rows. Correction of these bad pixels is implemented by simply interpolating across the bad columns using the PyWiFeS routine repair_blue_bad_pixels (and its red counterpart). These bad pixels are also flagged in the data quality extensions when the data is separated into the multi-extension FITS format. Fortunately these bad columns are typically far from any important galaxy emission lines in the red ($\lambda\sim$5550Å for R3000 and $\lambda\sim$5590Å for R7000) and despite potential impact in the blue ($\lambda\sim$4940Å for B3000 and $\lambda\sim$5040Å for B7000), these columns comprise only a few hundred km s^-1 gap in velocity coverage.

Any residual two-dimensional structure in the CCD bias level is removed using bias frames. Typically, this step is performed by median combining several bias frames into a superbias. This method is limited by the stability of the detector bias level, which was found to vary both spatially and temporally (on a time scale of minutes to hours at the level of a few e^-) in the first generation of detectors in WiFeS. We illustrate these variations in Fig. 2, by comparing 11 blue biases observed over a 13 hours interval corresponding to one observing night (including afternoon calibrations). As illustrated in the top panel, the shape of WiFeS biases is highly non-linear along the x dimension. The mean bias level is subject to variations of  1 e^- within the first hour, and again at the end of the night. The largest shape variations of the mean bias level occur for x positions around 3500 and beyond.

These spatial and temporal variations limit the precision of a standard superbias subtraction method. A solution to this problem, first implemented on WiFeS data (in a 4-amplifiers read-out mode) by Rich et al. (2010), is to fit a multi-dimensional surface to the bias taken closest to any given exposure, and use this locally reconstructed bias instead of a global superbias. This bias fitting method was later on officially added to the original IRAF data reduction pipeline for WiFeS, and is implemented in PyWiFeS with a slightly different algorithm.

A typical blue, raw superbias (y-binning of 2 and 1-amplifier read-out mode) is shown if Fig. 3 (top frame). The underlying structure, although somewhat masked by noise, is clearly visible. Red biases display a similar behavior, and are not illustrated here. Typically, a large number of bias frames are co-added (using the imcombine routine in PyWiFeS) to create a raw superbias, from which the bias structure is fitted in a two step process:

1. 1.

   The raw bias is collapsed and median averaged along the y dimension for all x position (see the top panel of Figure 2). This step was introduced first, to correct for the largest spatial variations which occur along the x dimension.

2. 2.

   The spatial variations along the y dimension are obtained by performing a 2D smoothing (using a symmetric gaussian kernel of 50 pixels) of the raw superbias minus the 1-dimensional correction obtained in step 1, and subtracted on a row-by-row basis.

In Fig. 3, we show in the middle panel the reconstructed superbias based on the raw superbias in the top panel. In the bottom panel, we show the residual after subtracting the reconstructed bias from the raw bias. As expected, the large scale spatial variations have been removed and the residual is flat at a sub noise-level.

In the default reduction script (see Section 4.2), this master reconstructed superbias is used as the default to perform the bias-subtraction step on any given exposure. It is usually constructed from a series of biases acquired during the afternoon. Individual bias frames acquired by the observer throughout the night can be associated with particular science frames (see Section 4.1), and will be used by the reduction script to create local reconstructed biases which are then subtracted from the requested science frame. Although these local reconstructed biases rely on individual bias frames, they are in fact noise free, and represent the most accurate way to perform the bias subtractions step with the WiFeS instrument with first-generation detectors, and best account for the spatial temporal variability of their bias levels.

It is unclear at this stage how the bias levels of the recently installed second-generations detector behave. Once they will have been characterised, PyWiFeS will be updated (if required) to ensure the most appropriate bias subtraction method is being used.

Cosmic ray (CR) rejection for WiFeS data is accomplished in PyWiFeS by means of a custom implementation of the Laplacian kernel technique devised by van Dokkum (2001). The PyWiFeS CR rejection routine performs additional operations on the data to account for the slanted (sometimes curved) wavelength solution of the instrument, as well as the multiplicity of operating on twenty five slitlets of data.

A schematic representation of the CR rejection steps is shown in Figure 5. The first step in identifying CRs is to subtract a smooth model of the sky background. This is complicated for WiFeS data by the slanted wavelength solution, so the sky model is derived in three steps: (1) resample the data to a rectilinear wavelength grid, (2) smooth with a box-median kernel along the cross-dispersion ($y$) direction, then (3) transform the smoothed sky spectrum back to the original pixel wavelength sampling. This sky model is then subtracted from the original data, and that subtracted data is convolved with a Laplacian kernel as devised by van Dokkum (2001). A threshold based on noise statistics of the data and detector read noise is then applied to the convolved data and those pixels above the threshold are identified as being contaminated by CRs. A slightly lower threshold is then applied to neighboring pixels to those identified as CR-contaminated. The data in all CR-contaminated pixels is then replaced with by interpolating from nearby CR-free pixels, and the entire process is repeated for the requested number of iterations (the default is three). The CR rejection procedure is performed for all 25 WiFeS slitlets, and can take advantage of parallel processing capabilities by simply passing the keyword argument multithread=True to the CR rejection function.

The format of data generated by integral field spectroscopy is flux values as a function of wavelength at each spatial location on the sky. This sampling of the flux is commonly referred to as a “data cube” as it samples flux in three dimensions. This data format requires the three coordinates $(x,y,\lambda)$ be derived for each pixel on the detector. While the $x$ coordinate is simply the slitlet number, determining a consistent $y$ coordinate for all slitlets requires a measurement of the spatial zeropoint in each slitlet. This is accomplished using a “wire” frame observation, which we describe in detail in Section 2.4.1.

The $(x,y)$ coordinate pairs calculated within the instrument field of view do not correspond to the same physical coordinates on the sky $(\alpha,\delta)$ at all wavelengths. This is due to the influence of atmospheric differential refraction (ADR – Filippenko, 1982), where the wavelength-dependent index of refraction of the atmosphere deflects light of different wavelengths by subtly different values. Calculation of this effect is a solved astrophysical problem, and we describe the implementation in PyWiFeS in Section 2.4.2.

For most applications, the irregular spatial and wavelength sampling realised by the actual WiFeS detector pixels is too complex to easily extract physical quantities over the full instrument FOV. Instead the observed fluxes are resampled onto a regular (rectilinear) grid of $(x,y,\lambda)$ for the final data cube, and this step for PyWiFeS is outlined below in Section 2.4.3.

Once WiFeS data have been corrected for pixel-wise variations of detector efficiency and spatial variations of the instrument throughput, the wavelength dependence of overall throughput of the instrument must still be corrected. This is achieved through observation of astronomical sources with known spectral energy distributions, referred to as spectrophotometric standard stars (see Bessell, 2005, for a thorough review). This process of measuring the absolute sensitivity of the instrument and correcting observational data to this absolute scale is known as flux calibration, and is achieved in PyWiFeS using routines defined in the wifes_calib sub-module. Flux calibration is typically accomplished in three steps:

1. 1.

   Standard star spectra are extracted from data cubes using the extract_wifes_stdstar routine.

2. 2.

   The instrument sensitivity function is derived by comparing observed standard star spectra to their reference spectra using the routine derive_wifes_calibration.

3. 3.

   The flux calibration solution is applied to uncalibrated WiFeS data cubes using the calibrate_wifes_cube routine.

Standard star spectra are extracted using a simple aperture extraction technique employed by the extract_wifes_stdstar routine. This is illustrated schematically in the left panel of Figure 11, where the sky background spectrum is estimated as the mean spectrum of spaxels outside some fiducial radius. After this sky spectrum is subtracted from each spaxel, the flux within some object aperture is summed to obtain the final object spectrum. The center of the extraction aperture is determined by fitting the flux centroid in the $(x,y)$ plane, which is generally a robust measure of the star’s location. The extraction and sky radii are set by default to 5″ and 10″, respectively, though these values (and the extraction center) can be modified by passing the desired values as keyword arguments to the extraction routine.

After the standard star spectra have been extracted, the instrument sensitivity function can be derived by comparing these to their reference spectra. In PyWiFeS this is accomplished in a (mostly) automated fashion using the derive_wifes_calibration routine. This function takes a list of standard star cubes as input (or optionally a list of extracted spectra) and extracts all spectra according to the algorithm described above.

The derive_wifes_calibration routine then attempts to identify the name of the standard star from the ’OBJECT’ header field of the data cube file (alternatively the user can pass a list of object names to this routine). These names are compared to a list defined in the wifes_calib.py file, comprising a Python dictionary containing the reference spectrum file name for each standard star. Future planned upgrades to PyWiFeS include the ability to cross-reference standard star lists by the object coordinates. Once the reference star name has been identified, its reference spectrum is loaded and the ratio of observed flux to reference flux is calculated for each star in the input list. This ratio is then corrected for the smooth atmospheric extinction using the default Siding Spring Observatory extinction curve measured by Bessell (1999).

The final flux calibration solution is fitted from the ’counts-to-flux’ ratio values (i.e., the sensitivity curves) obtained for all stars in the list passed to derive_wifes_calibration. Points falling in regions of strong telluric absorption are not included in the fit. The final sensitivity curve can be calculated with a single star or with a list of several stars. The sensitivity curves from all stars can be normalised to one another using the ’norm_stars’ keyword, which normalises all stars based on the middle 200 pixels, scaled to the maximum sensitivity curve. This allows a higher signal-to-noise measurement of the instrument sensitivity while not being hindered by grey extinction variations throughout the night.

Once the flux calibration solution has been derived, it can be applied to uncalibrated WiFeS data cubes using the calibrate_wifes_cube routine. This routine interpolates the derives flux calibration solution to the wavelength values of the input data cube, and multiplies the observed flux (in photon counts) to convert it to the true flux (in $F_{\lambda}$). This routine also applies the extinction correction, again based on the Bessell (1999) canonical Siding Spring extinction curve.

We note that the current flux calibration routine in PyWiFeS does not allow for a fit of residual extinction (i.e., deviation from the input extinction curve), though that capability is planned to be included in future upgrades to PyWiFeS. However, the star spectrum extraction routine is capable of saving the extracted spectrum in a format readable by the IRAF onedspec package, and the flux calibration routine can accept IRAF-style flux calibration solutions and extinction curves. Thus, users who require stringent flux calibration and extinction corrections can work with existing IRAF routines by employing the appropriate formatting keywords when calling the PyWiFeS extraction and flux calibration functions.

The broad smooth absorption of light by the Earth’s atmosphere is corrected using the atmospheric extinction curve in the flux calibration step in PyWiFeS. However, narrow structured absorption features caused by molecular species in the atmosphere (primarily molecular oxygen and water) persist in all object spectra after this step. These telluric features are most efficiently corrected by observing sources with smooth spectra at similar airmass and time of the night as a science target (e.g., using the smooth star division technique of Bessell, 1999). Alternatively, measuring telluric absorption in a number of smooth spectrum objects over a range of airmass can provide an acceptable estimate of the average telluric absorption behavior for a given night.

In PyWiFeS, telluric features are measured using the routine derive_wifes_telluric and corrected using the routine apply_wifes_telluric. As with the flux calibration routine, the derive_wifes_telluric routine takes a list of star cubes (or extracted spectra) as input. Regions of the object spectra not affected by telluric absorption are fitted with a low order polynomial, as demonstrated in Figure 12. The ratio of observed flux to the predicted smooth flux defines the telluric absorption for each spectrum.

Because telluric absorption strength depends on airmass, the absorption from each each spectrum in the list is corrected to airmass $secz=1.0$ and the mean absorption across all objects defines the final telluric absorption correction functions. The airmass dependence is a product of the saturation level of the absorption, which often differs between O₂ (which is more saturated) and H₂O (see, e.g., Buton et al., 2013). The logarithmic power law index for these two species can be set by the user as keyword arguments to the telluric fitting routine, and are stored along with the respective correction curves in a Python pickle (.pkl) file. From an observing night in 2013 we measured these indices to be 0.40 and 0.72 for O₂ and H₂O, respectively, which are both intermediate between the optically thick (saturated) value of 0 and the optically thin limit of 1. A more detailed study of the airmass dependence of telluric absorption features is planned for future work, but can be easily modified in practice using the keywords built into the routine.

Finally, once the telluric correction functions for O₂ and H₂O have been determined, they can be applied directly to data (with the appropriate airmass dependence) by calling the apply_wifes_telluric routine.

Discovering and fitting emission lines automatically from a large volume of data without user interaction requires a robust set of algorithms. For each row of each slitlet, we first identify all pixels which potentially have emission line flux in them with the following technique. We measure the mean and RMS of the background flux level (i.e., from pixels without emission line flux) via an iterative outlier rejection technique. Similarly, we measure the mean and RMS derivative of flux along the dispersion axis. Pixels must then pass a series of cuts to be identified as a potential emission line center: the flux value must be below the saturation level and above a certain flux threshold, the flux in the adjacent pixels must also exceed the flux threshold, and the derivative of adjacent pixels must exceed a threshold set from the flux derivative statistics (by default thresholds are the mean plus $3\sigma$ for flux and flux derivative as determined from the background statistics). Any adjacent sets of pixels which pass these cuts are grouped together and the potential line center is identified as the middle pixel value of the group.

Once potential emission lines have been identified, their exact centers must be measured. PyWiFeS provides two algorithms which can be chosen by setting the ’find_method’ keyword of the derive_wifes_wave_solution routine. The first, called ’loggauss’, is fast but less accurate, and is intended for rapid wavelength solution fits. It fits the logarithm of the emission line profile as a quadratic function (equivalent to a Gaussian in linear flux) but naturally must excise the pixels with zero flux from the fit. The second technique, called ’mpfit’, is somewhat slower but extremely accurate, and should be considered the standard for science-grade data reduction. This method performs a Gaussian fit to the emission line flux profile using the mpfit Python package (a Python port of the IDL implementation of the non-linear least squares fitting program MINPACK Markwardt, 2009). Though this procedure takes up to thirty minutes on a standard laptop, it is designed to naturally take advantage of parallel processing capabilities and generally takes about two minutes on the Linux clusters at RSAA. A visual comparison of these two methods is presented in Figure 13, which illustrates the sub-pixel errors in line centres found using the ’loggauss’ method as compared to the more robust ’mpfit’ method.

Once a series of emission lines have been identified in the data, the true wavelength of the atomic transition producing that line must be identified. Classically, this process is done interactively where an observer must visually inspect the arc lamp spectrum and compare it to a reference spectrum where lines have been identified. After identifying a few lines in the center of the longslit, a low order polynomial fit to the dispersion relation (i.e., wavelength versus pixel number) is used to automatically identify the remaining lines, and these lines are cross-identified in all remaining rows of the data. This visual inspection and line identification procedure was deemed to be cumbersome, especially given the excellent wavelength stability of the WiFeS instrument (thermal drift of less than a few angstroms).

In practice, emission lines in WiFeS arc lamp spectra are associated with a reference wavelength in a multi-step process. First a baseline guess for the wavelength at each found emission line position (in $x,y$ on the detector) is calculated from a known reference wavelength solution. A coherent shift of the wavelength guess for each row is applied by performing a cross-correlation of the spectrum for that row with a reference spectrum (this is optional but strongly recommended). Next, a comparison of the found line wavelengths to reference line wavelengths is performed (for each detector row), so that each found line is assigned a reference if and only if that reference wavelength is the closest wavelength in the reference list and the found wavelength is the closest found line to that reference line (i.e., a nearest neighbors requirement from both the found and reference lists). A baseline cut is applied at this stage so that associated line pairs with highly discrepant wavelength values (default is 5Å) are rejected as false associations.

Finally, once lines have been reliably measured and associated with a reference wavelength, the fitting of the wavelength solution proceeds. The technique employed in PyWiFeS is built on a global model of the instrument optics, which we describe in the following Section.

The manner in which light of different wavelengths is dispersed within the WiFeS instrument can be described analytically using the geometric optics of VPH gratings (Barden et al., 2000). Light from all spatial positions in the WiFeS FOV passes through the same VPH grating, but light from each spaxel has a different angle of incidence and physical point of entry due to the slit geometry from the image slicer (Fig 14a). The slicer separates the slitlets by reflecting each spatial slice with a differential in the vertical plane – spreading the light out into a pseudo-slit. Each slitlet is also offset in the horizontal direction by widths of the neighboring slits producing the staircase slit geometry. Slitlets are further offset in the vertical direction by their own length to allow the necessary storage space on the CCD for nod-and-shuffle observations.

The light from each slitlet is dispersed following the grating equation:

with a grating line density $a$ lines per mm. The off-axis angle, $\gamma$, introduces additional optical path difference for off-axis slits, resulting in significant spectral curvature towards the edge of the CCD. The staircase slit geometry also results in a slit-number dependence for angle of incidence $\alpha$ introducing a shift to the diffraction pattern with slit-number in addition to the classic spectral curvature.

The WiFeS optical model is described by a total of 43 parameters. The 18 primary parameters include values such as the grating line density and blaze wavelength, the primary angle of incidence of the optical beam, the centre of the beam on the detector (in $x$ and $y$), radial distortion terms, focal length of the camera, tilt of the detector (in $x$ and $y$), as well as other geometric terms. In addition to these, each of the 25 slitlets has some very small deviation (which we label $\Delta\alpha_{i}$) from uniform angular separation $\alpha$ between adjacent slitlets, caused by limits in the manufacturing precision of the image slicer. Though incredibly small ($\lesssim 10^{-5}$ radians), these offsets produce a coherent shift in wavelength of order 0.05Å for each slitlet.

The 43 WiFeS optical model parameters are fitted by minimising the residual offset between the optical model wavelengths and reference wavelengths for the successfully found emission lines. In practice this is done in a series of steps which each fit for a specific group of optical model parameters (e.g., the 25 $\Delta\alpha_{i}$ values are fitted together as the final group of parameters). Residuals from this optical model method are typically of order 0.05Å for the R=7000 gratings, and 0.1Å for the R=3000 gratings. Figure 15 shows the final diagnostic plots for the optical model wavelength solution fit for a R7000 grating Ne-Ar arc lamp exposure. Structure in the residuals is largely removed, with a final residual dispersion of 0.053Å.

In principle, many of the parameters in the WiFeS optical model should be constant over time (the $\Delta\alpha_{i}$ values, the beam centre and detector tilt, detector focal length, etc.). Analysis of large volumes of WiFeS data is currently underway to determine the stability of these parameters, and future versions of PyWiFeS will likely reflect knowledge of the static optical model parameters.

WiFeS is located on the Nasmyth focus of the 2.3m ANU telescope. This design ensures that the instrument is subject to a constant gravity vector, and makes it a very stable spectrograph overall. WiFeS can however be subject to temperatures changes of the order of $\sim$5-15 degrees during an observing night. When subject to a change in temperature, the WiFeS gratings will expand or contract, increasing or decreasing their resolution. To account for and correct this effect, observers are required to acquire arc frames during the night to monitor the changes of the wavelength solution. In Figure 16, we illustrate the effect of temperature drift on the wavelength solution.

In the top two panels, we show the spectra of a Ne-Ar arc lamp observed with the B3000 and R7000 gratings. We reduce the arc lamp exposure similarly to a science frame using PyWiFeS, and plot the 25 spectra extracted along the middle slice (along the y direction) of the final data cube. The arc frame has been acquired at the end of a winter observing night on 16th August 2012. The wavelength solution is derived from a similar arc exposure acquired early in the preceding afternoon. The reference wavelengths for the different arc lines are marked with vertical dashed lines. As expected, using an arc far-away from a given observation results in a bad wavelength calibration of the data set, with errors of the order of 1-2Å(typically 40-60 km/s). We note that the error results in a blueshift of the spectra for the blue frame, and a redshift of the spectra for the red frame. In the bottom plot, we show the same Ne-Ar spectra, but using the frame itself to derive the wavelength solution. As expected, the arc lines are in this case in perfect agreement with their reference wavelengths. We urge WiFeS observers to keep this point in mind when planning their observing run, and to take arc frames at regular time interval to monitor the temperature changes of the grating. Even if one’s science goals do not require an absolute wavelength calibration, regular arc frames are still required to obtain a wavelength solution consistent between the red and blue frames. Our experience suggests that an interval of $\sim$1 hour between arcs is a reasonable choice, but we leave it to every observer to decide which cadence is most appropriate to their science goals.

High-level information characterizing important properties of data are typically referred to as metadata. For spectroscopic observations such as those collected with WiFeS, examples of metadata include the observation type, the target being observed, or whether multiple exposures should be associated with one another. This information is then used to decide which data reduction procedures should be used, and is perhaps most commonly handled manually. Because spectroscopic data reduction often follows a repeated pattern, we desired an abstract metadata structure which could be assigned to a particular set of data then handled generically by a standard data reduction script.

In practice, the PyWiFeS metadata structure for a given set of data is stored in a Python dictionary which is saved in a Python pickle file. This dictionary is defined in the scripts save_blue_metadata.py and save_red_metadata.py (the metadata structures must be defined separately for the red and blue cameras), which we will call the “metadata saving” scripts. These are included in the PyWiFeS distribution and each script simply defines the metadata dictionary and saves it to a pickle file.

The metadata dictionary has several keywords which store information about the various types of calibration images, as well as science observation groups, and observations of standard stars. Generic instrumental calibration files – such as flat field observations, arc lamp exposures, and bias frames – which can be applied to calibrating all science and standard star data are stored in simple lists. These files are considered to be the “master” calibration files for a given data reduction run, and are applied in the calibration of all files unless specific calibration files are specified for a given observation.

Science and standard star observations are each stored as a list of Python dictionaries. Each dictionary describes an observation group, including all science frames which are ultimately co-added, as well as associated calibration frames (“local” calibrations) taken explicitly for use in calibrating that particular science target (e.g., a bias frame and arc frame taken immediately after observing a science target). Each keyword of the dictionary is a list, with the ’sci’ keyword containing a list of all observations to be co-added as the final science frame, the ’sky’ keyword containing a list (typically of length 1) of sky frames to be subtracted during data reduction, ’arc’ containing the “local” arc frame, and so on. Standard star dictionaries also have the special keyword ’type’ whose value is a list consisting of one or both of the strings ’flux’ and ’telluric’, used to distinguished if that star is used in the flux calibration and/or telluric correction procedures.

For example, if the only science observation is a pair of science frames followed by a sky frame, then a bias and an arc frame, the pertinent part of red metadata saving script would appear as this:

   sci_obs = [
       {’sci’  : [’r0001’, ’r0002’],
        ’sky’  : [’r0003’],
        ’bias’ : [’r0004’],
        ’arc’  : [’r0005’],
        },
   ]

This illustrates the naming convention employed in the metadata structure, where the observation root name is used to identify the data so that ’r0001’ refers to the observation whose raw data is in the file “r0001.fits”. At the end of the metadata script, the list ’sci_obs’ is saved as the ’sci’ keyword of the metadata dictionary. A summary schematic diagram of the PyWiFeS metadata structure is shown in Figure 17.

The metadata structure is saved to a Python pickle file by running the save_red_metadata.py script from the command line, and this pickle file is passed as input to the data reduction script. The metadata saving script is a critical part of the data reduction scheme, and should be carefully prepared by the user before each reduction.

For convenience, we include with PyWiFeS a Python script which inspects all FITS files in a directory and attempts to automatically generate the metadata saving scripts (for blue and red). We caution that this script is very simplistic, and identifies calibrations by the ’IMAGETYP’ header keyword set by the WiFeS controller. All calibration files are by default assigned to the “master” calibration lists, and any calibrations which the user took for specific science frames should be assigned by the user to those frames as in the example provided above. The PyWiFeS metadata structure is designed to be as abstract as possible to enable ease of data reduction, but requires careful construction of the metadata by the user in order to ensure data reduction is performed correctly.

The PyWiFeS data reduction scripts reduce_blue_data.py and reduce_red_data.py take a PyWiFeS metadata file as input and perform the data reduction procedures outlined by the user. The “standard” version of these scripts are included in the PyWiFeS distribution, but can be edited to add or modify data reduction steps as desired. In this Section we describe the default structure of the reduction script.

The data reduction procedure performed by the reduction script is defined in a Python list proc_steps, which is a list of the data reduction steps. Each step is defined by a Python dictionary which has four keywords:

• •

  step – This is the main label for each reduction step and identifies which function (described below) should be called. For example, the first step “overscan_sub” calls the function run_overscan_sub, which in turn calls the appropriate PyWiFeS functions to operate on the data.

• •

  run – A simple Python Boolean (True or False) which sets whether that step in the reduction is run. This keyword is useful when individual steps of the reduction need to be re-run, e.g., if the wavelength solution does not converge to the user’s satisfaction. It is important to note that a step which has been set to ’run:False’ will be interpreted by the reduction script as having been successfully executed, and the following step will look for the output of the “switched-off” step. If users wish to completely remove a step from the reduction chain, it should either be commented out or deleted.

• •

  suffix – This defines how the output of the reduction step will be labeled. For example, an image whose root name is ’r0001’ (as defined in the metadata file) and is processed through the first reduction step with ’suffix’:’01’ would have output “r0001.p01.fits” where the ’p’ indicates the file has been through the processing step labeled by ’01’. The reduction script knows to look for the suffix of the previous step to define the input image to that step. Following the same example, ’r0001’ passed to the next step with ’suffix’:’02’ would result in that step looking for “r0001.p01.fits”, processing it, and saving the output in “r0001.p02.fits”.

• •

  args – This keyword is set as a dictionary of arguments which should be passed to the main function being called in that data reduction step. For example, in the “cube_gen” step which generates final data cubes, one can set the minimum wavelength of the cubes by setting the args keyword as ’args’:\{’wmin_set’:5400.0\}. This will pass the argument on so that every time generate_wifes_cube is called in this step, it will be called with the keyword argument wmin_set=5400.0. This functionality allows the user to adjust the output of the PyWiFeS functions easily from the processing step list.

Each reduction step calls a “master reduction function” (such as run_overscan_sub in the example above) which is defined in the reduction script (which can be modified by the user). These functions must be defined to take as input the metadata structure (the Python dictionary described in Section 4.1), the suffix from the previous reduction step (prev_suffix), the suffix to be used in output from the current step (curr_suffix), and finally the additional arguments to be passed the pertinent PyWiFeS functions. A loop at the end of the reduction script automatically parses the processing step list (proc_steps) and calls the master reduction functions with the correct suffixes.

Each master reduction function parses the metadata in a specific way and calls the pertinent functions from PyWiFeS to perform that reduction step. For example, the run_overscan_sub function calls pywifes.subtract_overscan on all raw frames, run_obs_coadd calls pywifes.imarith_mef on all groups of science and standard star frames, and run_extract_stars calls pywifes.extract_wifes_stdstar on all standard star data cubes.

Some steps in the processing list do not perform operations on all data, but instead prepare master calibration files such as the flatfield response function or the wavelength solution. The master reduction function which generates these master calibration files also checks to see if any science observation groups require the use of “local” calibration files, and derives the appropriate calibrations from those files as well. The names of the master files are defined early in the reduction script (and again can be changed as desired by the user) and are global variables which are referenced when calibrations must be applied in the processing steps, but use of local calibrations is also handled automatically. For example, in the data cube generation step, science frames are transformed into data cubes using the master wire and master wavelength solution by default, but using the local wavelength solution for observation groups where a “local arc” was defined.

The full list of the default processing steps is:

Step 00 

Overscan subtraction (must happen before any other step).

Step 01 

Bad pixel repair.

(CAL) 

Creation of “superbias”, and bias model from superbias and any local biases.

Step 02 

Bias frame subtraction.

(CAL) 

Creation of “superflats”, both lamp and sky.

(CAL) 

Measurement of the slitlet profiles.

(CAL) 

Desag (optional)

(CAL) 

Flat frame cleanup.

(CAL) 

MEF flat creation, both lamp and sky.

Step 03 

MEF creation for all frames.

(CAL) 

Derivation of the wavelength solution.

(CAL) 

Derivation of the wire (spatial) solution.

(CAL) 

Derivation of the flatfield response function (requires wavelength solution).

Step 04 

Cosmic ray rejection.

Step 05 

Subtraction of sky frames.

Step 06 

Co-adding of frames in observation groups.

Step 07 

Application of flatfield response function.

Step 08 

Generation of data cubes.

(CAL) 

Extraction of uncalibrated standard star spectra.

(CAL) 

Derivation of flux calibration solution from standard star spectra.

Step 09 

Flux calibration of all data cubes.

(CAL) 

Extraction of calibrated, but not telluric-corrected, standard star spectra.

(CAL) 

Derivation of telluric correction from standard star spectra.

Step 10 

Correction of telluric absorption.

Step 11 

Reformatting into 3D data cube format.

Steps which perform processing operations on all science frames (and other frames requiring that reduction step) are labeled by number in the above list, while steps labeled as (CAL) denote processing steps which prepare or operate on master calibration files.

We also note that nod-and-shuffle (“N+S”) frames are handled automatically in the PyWiFeS reduction scripts. Special treatment of N+S data needs to occur when the MEF files are first created, when subtraction of sky frames occurs, and any intermediate steps between. In the default scheme described above, MEF files are created at the ’p03’ step, and N+S frames are automatically processed with the N+S version of the slitlet separation function. For example, if the raw WiFeS image was “r0001.fits”, the N+S version of this step would save the object pixels as the MEF file “r0001.p03.fits” and the nodded sky pixels are saved as the MEF file “r0001.s03.fits”. By default, these files are both processed in the cosmic ray step to produce “r0001.p04.fits” and “r0001.s04.fits” files, then in the sky subtraction step the nodded sky pixels in “r0001.s04.fits” are subtracted from the object pixels in “r0001.p04.fits” to produce the sky-subtracted “r0001.p05.fits”. All of these steps are executed automatically without the user needing to identify N+S frame, as these are identified from the ’WIFESOBS’ FITS header keyword.

Finally, we note that all PyWiFeS data processing functions are designed to work seamlessly with WiFeS data taken with detector binning, or in half-frame readout mode. However, the data reduction script and metadata creation script are not currently designed to handle multiple observing modes at the same time (though that functionality is planned for future upgrades). For example, attempting to subtract a full-frame bias from a half-frame image will cause the reduction to fail. We caution that these failures modes are not Python exceptions built specifically into PyWiFeS, but rather are the results of normal Python errors being raised (e.g., due to incompatibility of data size). This means that it is possible to perform reduction steps which are wrong but do not raise Python errors, such as subtracting a full-frame bias with binning of 2 from a half-frame image with binning of 1, or flux calibrating R7000 data with a calibration solution derived for the R3000 grating. Again, users should be careful in compiling their metadata so that only compatible data are being sent to the data reduction script.