SCUBA-2 Data Reduction Manual 

A new cookbook, sc21,  is available describing the reduction of SCUBA-2 data. This cookbook focusses on the Dynamic Iterative Map-Maker (or map-maker), which performs the processing of raw data into a science grade map, along with common post-processing steps.

 http://www.starlink.ac.uk/docs/sc21.htx/sc21.html

An explanation of how the map-maker works is given, including details of the default parameters which control each stage. Specialised configuration files which alter these parameters to suit different science goals are also introduced.

The cookbook guides users through all the post-processing options, from cropping and co-adding maps using PICARD, to applying an FCF and calculating the noise. The science pipeline is also discussed, with instructions for running it on a local machine.

Chapter List:
1. Introduction
2. SCUBA-2 Overview
3. Raw SCUBA-2 Data
4. The Dynamic Iterative Map-Maker
5. Reducing your Data
6. Tweaking the Configuration File
7. Examples of Different Reductions
8. The SCUBA-2 Pipeline
9. SCUBA-2 Data Calibration

Processing CLS data with ORAC-DR

A dedicated recipe has been written to process data from the Cosmology Legacy Survey using a "jack-knife" method of estimating and removing residual large-scale noise. The recipe is based on a script written by Ed for S2SRO data (and used to produce the figures in the SMURF paper). In this article I will describe how to use the recipe, what steps it performs and what parameters may be used to control the processing. In the future, this technique may be made more widely available for applying to other data sets (not just blank field ones).

How do I use it?

To use this recipe, initialize the pipeline, create a text file with the names of all the files to process and a file with recipe parameters if necessary. Then add the recipe name REDUCE_CLS to the usual command line:

% oracdr -loop file -files filenames.lis -recpars params.ini 
  -log sf -nodisplay REDUCE_CLS

By default the recipe uses the "blank field" makemap config file (though this may be overridden with the MAKEMAP_CONFIG recipe parameter; see below). The recipe creates a slew of output files with the following suffices:

1. _fmos - signal map (1 for each observation)
2. _mappsf - map-filtered PSF (1 for each observation). This is the same as the above signal map but with an artifical gaussian added to the time-series (and located at the map centre).
3. _wmos - coadded signal map
4. _wpsf - coadded map-filtered PSF
5. _jkmap - jack-knife map
6. _whiten - whitened signal map
7. _whitepsf - whitened, map-filtered PSF map
8. _cal - calibrated, whitened signal map
9. _mf - above map processed with matched-filter (using the _whitepsf map as the PSF)
10. _snr - signal-to-noise map created from above map (_mf)

Note that the per-observation files begin with s, and the coadds (and subsequent products) begin with gs. The ones you are probably most interested in are the _mf and _snr files.

What does it do?

The recipe works as follows.

1. Each observation is processed separately to produce a signal map (_fmos). In addition, each observation is re-processed with an artificial gaussian source added at the map centre (this will be used to create the "map-filtered PSF" image, _mappsf).
2. The signal maps are combined using inverse-variance weighting to create a coadded signal map (_wmos). The images with the artificial gaussians added in are also coadded to produce the "map-filtered PSF" (_wpsf).
3. The data are split into two groups made from alternating observations. Each group is coadded, and the two coadds are subtracted to produce the jack-knife map (_jkmap).
4. The SMURF command sc2filtermap is used to estimate the radial angular power spectrum (circular symmetry is assumed) within a region defined by twice the minimum noise in the coadded signal map. (The size of this region is written to the FITS header of the output map under the keyword WHITEBOX.)
5. The inverse of this angular power spectrum is applied to the coadded signal and PSF images (_whiten and _whitepsf) to remove residual low spatial-frequency noise.
6. The amplitude of the fake source in _whitepsf is compared with the input value to derive a corrected FCF. This new FCF is used to calibrate the whitened signal map to create the _cal map.
7. A matched filter is applied to _cal using _whitepsf as the PSF to create _mf.
8. Finally, a signal-to-noise ratio image is created (_snr).

What parameters are available?

The following recipe parameters can be used to control the processing:

• FAKEMAP_SCALE - Amplitude of the fake source (in Jy) added to the timeseries to assess the map-making response to a point source.
• MAKEMAP_CONFIG - Name of a config file for use with the SMURF makemap task. The file must exist in the current working directory, $MAKEMAP_CONFIG_DIR, $ORAC_DATA_OUT, $ORAC_DATA_CAL or $STARLINK_DIR/share/smurf.
• MAKEMAP_PIXSIZE - Pixel size in arcsec for the output map. Default is wavelength dependent (4 arcsec at 850 um, 2 arcsec at 450 um).
• WHITEN_BOX - Size of the region used to calculate the angular power spectrum for removing residual low-frequency noise in the data. Default is a square region bounded by the noise being less than twice the minimum value.

What if I want to run it again?

In addition to the full pipeline recipe, there is a PICARD recipe called SCUBA2_JACKKNIFE which performs all of the post-map-making steps. This can be used to examine the influence of including or omitting individual observations (say ones with visible artefacts that the pipeline is not able to trap), or investigate the effect of varying the size of the whitening region, or trim the images to a specific size before re-running the jack-knife steps.

Note that the pipeline must have been run once to produce all the necessary files which go in to this recipe. And note that the full pipeline recipe must be run again if a different config file or input gaussian amplitude is to be used.

The majority of the flexibility in controlling the processing occurs after the individual signal maps have been created. As a refresher, running this recipe would mean typing:

% picard -log sf -recpars params.ini SCUBA2_JACKKNIFE myfiles*.sdf

All of the control is through the recipe parameters in the file params.ini. The most important item to note is that you must provide a map-filtered PSF (otherwise a default will be used). However, you already have one of these: the file ending _mappsf.sdf listed above.

• PSF_MATCHFILTER - the name of the map-filtered PSF file
• WHITEN_BOX - Size of the region used to calculate the angular power spectrum for removing residual low-frequency noise in the data. Default is a square region bounded by the noise being less than twice the minimum value.

This is only the first version of the recipe, and its limitations have yet to be determined. However, one case I have found where the recipe fails is in a large, multi-field map with highly-variable noise. In this case it might be better to process each field separately before creating the mosaic (though this has the disadvantage that different filtering may be applied to each field). Try it out and see what you get. All feedback welcome.

FIT1D: A new SMURF command for ACSIS data

More time ago than I am willing to admit, I started coding a Starlink routine to fit spectral lines in ACSIS cubes. Got a long way until SCUBA-2 commissioning and calibration put a halt to it, but I have finally managed to finish the program, technically a beta-version, as part of the upcoming Kapuahi release of Starlink.

Usage:

% fit1d  in  out  rms  [config]  [userval]  [pardir]  [parndf]  [parcomp]

What distinguishes FIT1D from other profile fitting routines is that it specifically attempts to deal with two issues: non-gaussian profile shapes and the fact that ACSIS data-cubes have many, potentially very different, profiles to fit. Regarding the latter, there are many fitting routines that produce fitted profiles for data-cubes, but FIT1D also produces cubes with the fitted parameters themselves and can use such files as input to give control over the fit of, in principle, each individual profile in the cube. Thus it is e.g. possible to fit broad lines on the nucleus of a galaxy and narrow lines everywhere else. More about that below.

A section has been added to the SMURF documentation in SUN/258 about FIT1D, which I will try to summarize here. FIT1D is generic in that it can fit profiles along any axis of an up to 7-dim hyper-cube, but will be discussed here in the context of a default RA-Dec-Vel ACSIS cube. Note that the routine assumes that data have been baseline-subtracted, using e.g. MFITTREND, i.e. that the profiles have a zero-level at 0. 

Gauss-Hermite shapes as a function of the 3rd-order
      skewness coefficient 'h3' and the 4th-order the kurtosis (peakiness)
      coefficient 'h4'. The red box indicates the limits on acceptable
      values for h3 and h4 as defined in the defaults configuration file. Note
      that the fitted profile by default is restricted to positive values
      and will omit the shown negative features.

Non-gaussian Profiles.

FIT1D essentially re-implements the fitting code for non-gaussian profiles from the GIPSY package (Kapteyn Institute, Groningen, The Netherlands). Function types that can be fitted are Gaussian, Gauss-Hermite, and Voigt profiles. In particular, Gauss-Hermite functions are a powerful extention when fitting profiles that are skewed, peaky, or only approximately gaussian. The figure above shows Gauss-Hermite profiles as a function of the skewness coefficient h3 and the kurtosis (peakiness) coefficient h4. See  for SUN/258 further details, but note that the default setting in the configuration file is for FIT1D to suppress the negative features in the fitted profiles and to leave only the positive part of Gauss-Hermites.

Because of their ability to fit distorted shapes, Gauss-Hermites are particularly well suited to "capture" the maximum amount of emission from a cube. The fits can be remarkably accurate as is shown in the the figure below showing a 3-component fit (i.e. up to 3 spectral lines) using gausshermite2 functions (i.e. fitting both h3 and h4). Collapsing the resulting cube with fitted profiles can thus result in an accurate and almost noise-free white-light or total-emission map.

Fit1d - Black: original profiles; Red: results of a
    3-component Gauss-Hermite2 fit (fitting both h3 and h4)

FIT1D derives its ability to fit a complex line-shape both from the Gauss-Hermite function but also from that it can fit multiple (sub) components to get the best match possible. However, that can make the interpretation of the fits in terms of the physical characteristics and quantities difficult, hence for those you may also want to make a fit of the line-shape using a single standard Gaussian function. 

Component Parameter files

Besides a data-cube with the fitted profiles FIT1D also outputs so-called Component parameter files as NDF extensions in the header of the output file. These can also be copied out as independent data-cubes. There is a file for each component (i.e. line) that was fitted along the profile up to the number of components requested by the user. Each plane of a Component parameter file has an image of the value of a fitted parameter across the field-of-view. For instance, the one resulting from a gaussian fit has images respectively showing the fitted Amplitude, Position (velocity), and FWHM as well as a plane with an id-number of the function used.

Much of the (anticipated) use of FIT1D derives from the fact that Component parameter files can be used as input as well: either to provide initial estimates or fixed values to the fitting routine.  The difference between values specified in the Component parameter files
and ones declared in a User parameter values file is that the former can vary across the field-of-view whereas the latter will result in the same value being used for all profiles. E.g. for use with spectral-line surveys the User parameter values file can be used to provide initial estimates of the frequencies or velocities at which lines are expected or to fix fits at those frequencies.

By manipulating Component parameter files e.g. resulting from an initial fit, the user can customize or correct subsequent fits. In extrema, a Component parameter file could be made from scratch based on a model and be used to create a spectral-line data-cube with that model (config option: model_only=1) or be used as initial estimates for a fit. Of more practical use, Component parameter files can be used to correct problems associated with a fit since the art of fitting is not in the fitting algorithm, but in providing accurate initial estimates. For instance, the left image below shows a section of an Amplitude plane of a fit where there are problems in a few locations. Setting these location to bad values and using FILLBAD to interpolate over them, the corrected Component parameter file was used as initial estimate for a subsequent fit, resulting in the image on the right

Fit1d - Left: Section of a parameter file showing
      originally fitted amplitudes; Right: Amplitudes after using a
      corrected parameter file from the original fit as initial estimates
      for a subsequent fit.

More creative options are possible: after an initial fit with a gaussian, the function id can be changed to a gausshermite1 in part of the field and the resulting file used as initial estimates for a subsequent fit to account for skewed profiles there. Similarly, the initial guess of the FWHM can be made wide on e.g. nucleus of a galaxy while leaving it more narrow outside. As another example, the fit of multiple components can be limited to only part of the field by setting the parameter file for the second and higher components to bad values outside the relevant region (multiple component parameter files can be used as input: one for each component to be fitted).

In conclusion: please remember that this is a beta-release and that you may run into unanticipated issues. Also chosen limits in the configuration file may need tweaking. If an initial fit looks poor, try adjusting minamp (in units of rms!) or, in particular, minfwhm (in units of pixels!) in the configuration file (see: $SMURF_DIR/smurf_fit1d.def). Also use range to limit the fit to a relevant region of the spectrum.

The released implementation of FIT1D can fit up to 7 components per profile per run, but the output of multiple runs each covering a range in velocities or frequencies can be combined. The fit itself is fully multi-threaded and will be much faster on a modern multi-core computer: a 3-component gausshermite2 fit of 1.1 million spectra (a 2 Gb input file) took 15 minutes on a dual-core, 16 Gb memory machine versus 4 minutes on one with 12 cores and 75 Gb of memory.

Happy fitting!

Remo