Chapter 8
Post-processing Reduction Steps

 8.1 Flux conversion factors
  8.1.1 Aperture flux
  8.1.2 Peak flux
  8.1.3 Determining your own Flux conversion factors
  8.1.4 Applying the FCF
 8.2 Cropping your map
 8.3 Co-adding multiple maps
  8.3.1 Registering maps
 8.4 Sensitivity
  8.4.1 Getting the noise
  8.4.2 Map statistics
  8.4.3 Viewing the noise histogram
  8.4.4 Examining the error map with GAIA
 8.5 Regridding your data
 8.6 Displaying masks
 8.7 Point-source extraction: the matched filter
 8.8 Clump finding
 8.9 Map provenance & configuration parameters

8.1 Flux conversion factors

When your data comes out of the map-maker it is in units of picowatts (pW). A flux conversion factor, or FCF, needs to be applied to scale your data from units of pW to janskys (Jy). For more information on calibrating SCUBA-2 data see Dempsey et al. (2013) [8].

Below are our default FCFs.






APERTURE
PEAK




FCFarcsec (Jy/pW/arcsec2)
FCFbeam (Jy/pW/beam)




  450 μm 850 μm   450 μm 850 μm




4.71 ± 0.5 2.34 ± 0.08 491 ± 67 537 ± 24





IMPORTANT NOTES:

8.1.1 Aperture flux

To get the flux density of extended sources with aperture photometry you should apply the FCFarcsec. You can then sum the emission in an aperture. FCFarcsec was determined using a 60-arcsec diameter aperture. If your aperture differs from this you should scale your flux accordingly—the scaling factor can be read off the curve of growth (see Appendix F). This graph gives the ratio of aperture flux to total flux for a range of aperture diameters.

FCFarcsec is determined using 1-arcsec pixels. For different pixel sizes you will also need to multiply by the pixelsize squared.

8.1.2 Peak flux

If you want to read off the peak flux from your map, you should apply the FCFbeam (also known as the peak FCF). When you open your map in Gaia the value of the brightest pixel will be the peak flux of your source. (If your source is point-like, the peak value in the map is the total flux density). Applying FCFbeam will result in a map with units of Jy/beam. For point-like or compact sources smaller than the beam (with a Gaussian profile), this peak value will be the flux density of your source.

8.1.3 Determining your own Flux conversion factors

Calibration observations are taken at various points throughout the night and we recommend you compare the FCF calculated from your nearest calibrator to the standard FCF. You can use this FCF to scale the the standard value accordingly. You should chose the calibrator nearest in time to your science observations, unless the weather has changed significantly.

(1)
Reduce the calibrator with the map-maker using dimmconfig_bright_compact.lis as the configuration file.
(2)
Determine the FCF value by passing the output map to the Picard recipe SCUBA2_CHECK_CAL.
  % picard SCUBA2_CHECK_CAL 850calibrator.sdf

This will produce a log file (log.checkcal) which records the both FCFbeam and FCFarcsec. Check that these values closely approximate the standard FCF values given above.

(3)
Re-reduce the calibrator using the map-maker with the same configuration file you used for your science observations.
(4)
Again determine the FCF using SCUBA2_CHECK_CAL.
(5)
Apply this FCF to your reduced science data using cmult—see the next section. Remember to apply the pixel size squared factor if using FCFarcsec.

8.1.4 Applying the FCF

You can apply the FCF using the Picard recipe CALIBRATE_SCUBA2_DATA. By default this will multiply your map by 1000×FCFbeam, obtained from the start of Section 8.1. This produces a calibrated map with units of mJy/beam.

You can supply a parameter file if you wish to use a different value for the FCF or use FCFarcsec. The recipe will write out the calibrated file with an _cal suffix and will change the units in the map header.

  % picard CALIBRATE_SCUBA2_DATA  mapinpW.sdf

Note that the recipe will also take account of the pixel size when applying FCFarcsec.

8.2 Cropping your map

The nature of the scan patterns results in SCUBA-2 maps significantly larger than the requested size. The high noise towards the outer edges is a consequence of the scanning pattern. Although this excess data are down-weighted during reduction by the map-maker, you may wish to remove it before either combining maps (see Section 8.3) or publishing your map.

You can crop your map to the map-size set in the data header or to any requested size of box or circle using the Picard recipe CROP_SCUBA2_IMAGES. The centre of the cropping area will always be the centre of your map.

  % picard -recpars mypar.lis CROP_SCUBA2_IMAGES map_cal.sdf

The example above includes a parameter file specifying the radius of the circle to be extracted (in arcsecs). The format for the parameter file is shown below.

  [CROP_SCUBA2_IMAGES]
  MAP_RADIUS = 1800.0
  CROP_METHOD = circle

If this parameter file is omitted it will default to a box of sides equal to the map size in the header (as requested in the MSB). The output from CROP_SCUBA2_IMAGES is a file with the suffix _crop. Full details of this recipe can be found in the Picard website1

Tip:
The default crop shape will be a square. Avoid losing good data by specifying a circle using the parameter file.

8.3 Co-adding multiple maps

You may have multiple maps of the same source which you would like to co-add. Picard has a recipe called MOSAIC_JCMT_IMAGES that co-adds maps while correctly dealing with the exposure time and weights NDF extensions. The images are combined using inverse-variance weighting and the output variance is derived from the input variances.

  % picard -recpars mypar.lis MOSAIC_JCMT_IMAGES 850map*_cal_crop.sdf

This creates a single output file based on the name of the last file in the list, and with a suffix _mos.

There are a number of options associated with MOSAIC_JCMT_IMAGES (see the Picard manual for a full description). However, the main one is choosing between wcsmosaic (default) and the Ccdpack option makemos for the combination method. For more information on makemos and advice on choosing the best method see SUN/139.

The example parameter file below chooses makemos using a 3-σ clipping threshold.

  [MOSAIC_JCMT_IMAGES]
  MOSAIC_TASK = makemos
  MAKEMOS_METHOD = sigmas
  MAKEMOS_SIGMAS = 3

Currently there is no advantage in terms of data quality to reducing all observations simultaneously or separately. However, the latter does allow the option of assessing the individual maps before co-adding.

Tip:
The list of files can be the output from cat. Remember to include the back quotes. For example, % picard MOSAIC_JCMT_IMAGES `cat myfiles.txt`.

8.3.1 Registering maps

You can register a series of SCUBA-2 maps to a common reference position using the Picard recipe SCUBA2_REGISTER_IMAGES. This is only possible if a there is a common, known source that is present in all of the input maps. This should be done before combining your maps.

  % picard -recpars myparams.ini SCUBA2_REGISTER_IMAGES ‘cat listoffiles.txt‘

Here the parameter file contains the equatorial position of the reference source as in the example below. See SUN/265 for more details.

  [SCUBA2_REGISTER_IMAGES]
  REGISTER_IMAGES = 1
  REGISTER_X  = HH:MM:SS.S
  REGISTER_Y  = DD:MM:SS.S

REGISTER_X and REGISTER_Y may also be galactic longitude and latitude respectively, both measured in decimal degrees.

8.4 Sensitivity

8.4.1 Getting the noise

You can use the Picard recipe SCUBA2_MAPSTATS get the noise. This recipe estimates the RMS from both the map, the NEP, and the RMS predicted by the Integration Time Calculator. It then writes out a series of results in a log file called log.mapstats. The parameters written to this file are listed in Appendix C.

  % picard SCUBA2_MAPSTATS map.sdf

This recipe will report the noise in the same input units provided.

Note: SCUBA2_MAPSTATS is only designed to work on reductions of a single observations. On coadded observations it could produce misleading results, or even fail completely to work.

If multiple files are run through SCUBA2_MAPSTATS, either in a single call of PICARD or by repeatedly running PICARD in the same terminal on different files, the results will be appended to the existing log.mapstats file. The final columns — project, recipe and filename — are given to ensure it is clear to users which line of the logfile corresponds to which input file.

Steps for getting your coadded-map noise

After applying any necessary FCF, SCUBA2_MAPSTATS simply executes the following steps to get the map noise.

(1)
Crop the map to remove the noisy edges.
  % picard CROP_JCMT_IMAGES map_cal.sdf
(2)
Run the Kappa command stats to extract the median value from the error array.
  % stats map_cal_crop comp=err order

Tip:
Use the error array to avoid contamination of the noise distribution from bright sources.

8.4.2 Map statistics

The Kappa commands histat and stats are very similar and both return a range of statistics describing any NDF. In addition to the main data array they can be passed the error (comp=err), variance (comp=var) or quality (comp=qua) arrays (if available).

The reported statistics include the pixel maximum and minimum, standard deviation, number of pixels used and omitted, along with pixel mode, and mean. If you supply the order keyword, the median is also shown.

  % stats comp=err map_cal_crop order

Note: the standard deviation of the data array will give a similar result to the mean/median of the error array except with additional contamination from sources.

8.4.3 Viewing the noise histogram

You can view a histogram of the error array with the Kappa command histogram. Again comp=err must be specified.

  % histogram map_cal_crop comp=err numbin=200 style="color=white"

The output is shown in the Figure 8.1. For more information on the options for histogram see SUN/95.


pict

Figure 8.1: The error array viewed with Kappa command histogram.


8.4.4 Examining the error map with GAIA

It is also useful to view the error map itself. Open your reduced map in Gaia, then select the Error button on the Select NDF in container file window—see Figure 8.2. You will need to adjust the scaling to view the error map properly.

To assess the noise using Gaia, go to the toolbar on the main window and click on Image-Analysis Image regions. Next select the region shape you would like to check and draw it on your map by clicking and dragging the mouse. Click the Stats selected button in the Image regions window to get a report of the statistics in the selected region.

Tip:
You can write out the error array of your map into a new NDF using the Kappa command ndfcopy. For example, % ndfcopy map comp=err map_err


pict

Figure 8.2: The error array displayed with Gaia. After clicking the Error button (circled in red) you will have to rescale your map.


8.5 Regridding your data

To change the pixel size use the Kappa command compave. The following example increases the pixel size from 4 arcsec to 8 arcsec by using a compression factor of 2.

  % compave map map_regrid 2

Tip:
Remember you can use ndftrace if you are unsure of the pixel size.

8.6 Displaying masks

SCUBA-2 maps created using “AST-masking” (see Section 3.5) will contain a Quality array indicating the background pixels that were masked (i.e. forced to zero). In itself this is fairly simple - background pixels have a non-zero Quality value and source pixels have a Quality value of zero. However, it is also possible to have independent masks for the FLT and COM models, in addition to the AST mask. For instance, maps created using dimmconfig_bright_extended or dimmconfig_jsa_generic will have Quality arrays that contain both a FLT and an AST mask, and so some care needs to be used when interpreting the Quality array.

Each value in the Quality array is restricted to taking integer values between 0 and 255, and so can be thought of as 8 separate bits. Each “bit plane” within the Quality array holds a single mask - AST, FLT or COM. The showqual command can be used to find out which mask is held by which bit plane:

  % showqual fred.sdf
     AST (bit 1) - "Set iff AST model is zeroed at the output pixel"
     FLT (bit 2) - "Set iff FLT model is blanked at the output pixel"

This means that the AST mask is stored in bit 1 (the least significant bit), the FLT mask is stored in bit 2, and there is no COM mask. Note, if a map was produced using FLT masking but no AST masking, then the FLT mask would be stored in bit 1.

The decimal integer value of any element of the Quality array is equal to the binary value formed from the bits listed by showqual. So in the above case the maximum Quality value is 3 (the decimal equivalent of binary “11” — both bits set). Remembering that a bit is set (i.e. is 1) for background pixels and cleared (i.e. is 0) for source pixels, it follows that the four possible decimal Quality values in the above case (0-3) are:

(1)
- neither bit set, so the pixel is inside both the AST and the FLT mask (a source pixel).
(2)
- bit 1 set but not bit 2, so the pixel is outside the AST mask but inside the FLT mask (a border-line pixel).
(3)
- bit 2 set but not bit 1, so the pixel is inside the AST mask but outside the FLT mask (a border-line pixel).
(4)
- both bits set, so the pixel is inside neither mask (a background pixel).

Figure 8.3 shows a simple map of the Quality array values — the black areas have value zero and are thus inside both masks, the dark brown areas have value 2 and are thus inside the AST mask but outside the FLT mask. The light brown areas have value 3 and are inside neither mask. In this particular case, there are no areas with a quality value of 1, so the FLT mask is contained entirely within the AST mask.


pdfpict
Figure 8.3: The quality component from a typical map.


To produce a plot like this, first open your map in Gaia. Select the top level NDF from list in the pop-up window then click the Quality button. You will need to rescale the map to bring out the masks–see Figure 8.4 for another example showing the full Gaia window (using rather brighter colours this time!).


pict

Figure 8.4: Using Gaia to display the mask applied by the map-maker. Select the QUALITY component of your map.


There are several commands within Kappa that manipulate Quality arrays in various ways. For instance, the setbb command allows pixel data values to be set bad if the associated Quality value has a specified collection of set bits. Thus:

  % setbb fred 1

will set all pixels bad in fred.sdf except for those inside the AST mask. Likewise,

  % setbb fred 2

will set all pixels bad except for those inside the FLT mask. Note, the change made by setbb is temporary - it can be undone by doing:

  % setbb fred 0

To display the fred.sdf map and then overlay the AST mask in blue and the FLT mask in red, do:

  % gdclear
  % display fred mode=perc percentiles=\[2,98\]
  % setbb fred 1
  % contour fred clear=no mode=good labpos=! style=’colour=blue’
  % setbb fred 2
  % contour fred clear=no mode=good labpos=! style=’colour=red’
  % setbb fred 0

The resulting plot is shown in Figure 8.5.


pdfpict
Figure 8.5: The AST and FLT masks contoured together over an image.


Alternatively, Gaia can be used to contour the Quality array over an image, but you cannot then distinguish the AST mask from the FLT mask. First you will need to copy out the QUALITY component of the data to a separate file using ndfcopy:

  % ndfcopy comp=qua map map_mask

Then open your map in Gaia and contour the mask NDF on top. Select Contouring... from the Image-Analysis menu, and supply the name of the mask NDF either by entering its name in the Other image: box, or selecting with Choose file... file browser. See Figure 8.6.

For more information on the use of masks by makemap and the parameters that affect them see sectyion 3.5.


pict

Figure 8.6: Using Gaia to display your map with the AST mask used by the map-maker contoured on top.


8.7 Point-source extraction: the matched filter

This effectively fits a single Gaussian point spread function (PSF), centered over every pixel in the map, and applies a background suppression filter to remove any residual large-scale noise.

Cosmology maps usually contain very faint sources that often need extra help extracting. The Picard recipe SCUBA2_MATCHED_FILTER can be used to improve point-source detectability.

The matched filter works by smoothing the map and PSF with a broad Gaussian, and then subtracting from the originals. The images are then convolved with the modified PSF. The output map should be used primarily for source detection only. Although the output is normalised to preserve peak flux density, the accuracy of this depends on how closely the real PSF matches the telescope beam size. In the case of nearby sources, each ends up contributing flux to both peaks.

  % picard -recpars mypar.lis SCUBA2_MATCHED_FILTER 850_map_cal_crop.sdf

As in the example parameter file below we have requested the background should be estimated by first smoothing the map and PSF with a 15-arcsec Gaussian.

  
  [SCUBA2_MATCHED_FILTER]
  SMOOTH_FWHM = 15
  

This is a fairly common technique used throughout the extra-galactic sub-millimetre community to identify potential sources. A full description of the matched filter principle is given in Appendix D, while the Picard manual gives full details of all the available parameters.

8.8 Clump finding

The Cupid application findclumps can be used to generate a clump catalogue. It identifies clumps of emission in one-, two- or three-dimensional NDFs. You can select from the clump-finding algorithms “FellWalker”[2], “Gaussclumps”, “ClumpFind” or “Reinhold” and must supply a configuration file specific to each method. See the Cupid manual for descriptions of the various algorithms.

The result is returned as a catalogue in a text file and as a NDF pixel mask showing the clump boundaries.

  % findclumps in=S2map.sdf out=clumpmap.sdf outcat=clumps.FIT logfile=clumps.log \
    config=^config.dat method=fellwalker rms=25 shape=polygon

The shape option allows findclumps to create an STC-S description (polygonal or elliptical) for each clump. These are added as extra columns to the output catalogue.

Polygon Each polygon will have, at most, 15 vertices. For two-dimensional data the polygon is a fit to the clump’s outer boundary (the region containing all good data values). For three-dimensional data the spatial footprint of each clump is determined by rejecting the least significant 10% of spatial pixels, where "significance" is measured by the number of spectral channels that contribute to the spatial pixel. The polygon is then a fit to the outer boundary of the remaining spatial pixels.
Ellipse All data values in the clump are projected onto the spatial plane and "size" of the collapsed clump at four different position angles—all separated by 45—is found. The ellipse that generates the same sizes at the four position angles is then found and used as the clump shape.


pict

Figure 8.7: STC shapes displayed over the input map with Gaia. To overlay your clumps in this way display your map with Gaia. Under the Image-Analysis menu select Positions... followed by Import CUPID catalogue... In the pop-up window select the .FIT file generated by findclumps and ensure the STC shape box is checked before importing.


8.9 Map provenance & configuration parameters

You may want to check a reduced file to determine which raw files went into it and what configuration parameters were used with the map-maker. There are a number of useful Kappa commands to help you with this:

Map provenance
provshow will list all the NDFs and operations that were used in the creation of your map, while hislist will show a truncated list of input files and operations but will list every configuration parameter used by the map-maker. This includes all the hidden default values as well as those specified in your supplied configuration file. These commands are executed like this:
  % provshow 850map
  % hislist 850map

The output of provshow can be very long. If you merely want to know what were the original file, the ones with no parents, select the root ancestors.

  % provshow 850map show=roots
Map-maker parameters
The task configmeld will compare two maps and highlight any differences between the configuration parameters that were used to create the maps. A configuration file can be supplied in place of a map, in which case the parameters used to create the map are compared with those in the configuration file. configmeld requires a visual-differences tool to be installed on your machine; those currently recognised are: meld opendiff, diffmerge, kdiff3, tkdiff, and diffuse, and are searched for in that order. The follow example illustrates two ways to run configmeld:
  % configmeld 850map1.sdf 850map2.sdf
  % configmeld 850map1.sdf ^mydimmconfig2.lis

Another useful command is the Kappa command configecho. This is very versatile and will display the name and value of one or all configuration parameters either from a configuration file or from the history of an NDF.

The first example below will return the value of numiter from the map 850map.sdf. The second example will display the values of all the parameters in 850map.sdf and will prefix them with a ‘+’ if they differ from the the default parameter values defined in $SMURF_DIR/smurf_makemap.def.

  % configecho name=numiter config=! ndf=850map.sdf
  % configecho name=! config=^$SMURF_DIR/smurf_makemap.def ndf=850map.sdf

1http://www.oracdr.org/oracdr/PICARD