A new Starlink release contains notable updates to the SCUBA-2 configuration files.

The latest Starlink release - 2014A - has been made public. For details please read the release notes provided at: http://starlink.jach.hawaii.edu/starlink/2014A

As part of this new release we want to highlight one significant update and a couple of new additions to the SCUBA-2 reduction arsenal of config files.

Updates to bright_extended config file

This config file 'dimmconfig_bright_extended.lis' has always been intended for reducing data containing bright extended sources. It has remained untouched for a couple of years now despite advances in our understanding of SCUBA-2 reduction of bright regions. The config file now contains the following parameters/links:

   ^$STARLINK_DIR/share/smurf/dimmconfig.lis

   numiter=-40
   flt.filt_edge_largescale=480
   ast.zero_snr = 3
   ast.zero_snrlo = 2

   ast.skip = 5
   flt.zero_snr = 5
   flt.zero_snrlo = 3

In previous Starlink releases i.e. Hikianalia, the bright_extended configuration file only contained the following:

   numiter = -40
   ast.zero_snr = 5

   flt.filt_edge_largescale = 600

New 'FIX' config file

Two new config parameter files have been added These are intended to be used with one of the existing dimmconfig files. They provide new values for selected parameters aimed at solving a particular problem ("blobs" in the final map, or very slow convergence).

• dimmconfig_fix_blobs.lis 
• These parameters attempt to prevent smooth bright blobs of emission appearing in the final map. It does this by 1) identifying and flagging samples that appear to suffer from ringing a soft-edged Butterworth filter in place of the normal hard-edged filter, and 3) rejecting samples for which the separate sub-arrays see a markedly different common-mode signal.

• dimmconfig_fix_convergence.lis
• The parameters defined by this file attempt to aid the convergence process, and should be used for maps that will not converge within a reasonable number of iterations.

CHECK_RMS: comparing SCUBA-2 map noise with the ITC

There's a new SCUBA-2 analysis tool in town, and its aim is to help observers get a feel for whether their maps are reaching the intended sensitivity levels. It comes in two forms:

• PICARD recipe: SCUBA2_CHECK_RMS
• Offline pipeline recipe REDUCE_SCAN_CHECKRMS

At its core, the "CHECK_RMS" analysis estimates the RMS and NEFD for a given observation to compare with values from the SCUBA-2 integration time calculator (ITC). To obtain access to these new recipes, you must be using /stardev in Hilo or update your own version using rsync (see http://starlink.jach.hawaii.edu/starlink/rsyncStarlink).

PICARD recipe: SCUBA2_CHECK_RMS

As usual, the PICARD recipe is designed to be run on processed data. To run it, type:

% picard -log sf -nodisplay SCUBA2_CHECK_RMS myfiles*.sdf

in the directory containing the files to be analyzed. You may notice a warning (cyan text) about missing NEP data: that doesn't affect the main analysis and can safely be ignored. However, if run in Hilo, the recipe can query the archive of log files generated by the QL pipelines for NEP data. These NEP data are used to determine an effective noise and NEFD from the timeseries noise calculations, which can then be compared with the other values.

The recipe can be given calibrated or uncalibrated maps: the default FCF will be applied to uncalibrated data. Note that these maps must be made from a single complete observation, either by hand with makemap or the offline pipeline. Single observations must be used because the recipe relies on FITS header information that will be incorrect for coadded files. Note that maps made by the SUMMIT pipeline should not be given to this recipe; the SUMMIT pipeline will soon start recording the same CHECK_RMS data itself.

The input images are cropped to a circle of radius 90 arcsec from which the median noise and NEFD (the values RMS_MAP and NEFD_MAP in the log file) are calculated. The mean exposure time is estimated over the same map area (TEXP). It's important to note that the RMS_MAP is derived from the variance component in the map, and not from the map itself. This means that the results should be largely independent of the chosen makemap config parameters, and are more likely to reflect changes in the number of samples that go into the map rather than any differences in spatial structure introduced by applying different high-pass filters in the map-maker.

The map size and scan type are read from the FITS header and used to define the observation type for the ITC. The RMS corresponding to this observation type is calculated from the SCUBA-2 ITC (RMS_ITC); the integration time is used to convert that to a corresponding NEFD (NEFD_ITC). However, it is important to note that if you have chosen a non-standard map size for a PONG then this recipe will be unable to calculate values from the ITC, rendering it somewhat less useful.

If the NEP data are available (in Hilo only, and for dates after 20120515), the average NEP for a given observation is calculated from the archived QL pipeline log file (log.nep) corresponding to the date and wavelength of the input map. The FCF (either from the file or the standard value) and mean transmission over the observation is used to convert this to a zenith NEFD (NEFD_NEP). The ITC scaled elapsed time is used to convert that to an RMS noise (RMS_NEP).

The results are written to a log file called log.checkrms, the format of which is described below.

The following recipe parameters are supported:

• KEEPFILES - if set to 1, the _crop files created during processing will be retained on disk, otherwise all intermediate files will be deleted.
• MAP_RADIUS: radius (in arcsec) of the region used for calculating the map noise and NEFD. Default is 90 arcsec. Note that poor results can be derived if the map radius is too small and it may be more useful to use a larger radius for larger maps.
• STATS_ESTIMATOR: Statistical estimator for the NEFD and RMS calculated from the input map. May be mean or median - default is median.

Offline pipeline: REDUCE_SCAN_CHECKRMS

An alternative to the PICARD recipe, this recipe performs all the necessary steps on the raw data to do a complete self-consistent analysis for each observation, from calculating the timeseries noise to processing the data and comparing results with the ITC. The steps performed are:

• Calculate timeseries noise from the entire observation;
• Determine various values from the FITS headers (elapsed time, average airmass, elevation, optical depth, transmission, date etc);
• Make a map using the given config file (default if not specified);
• Calculate the noise (from the variance component), NEFD and exposure time within a 3-arcmin diameter circle at the map centre;
• Use the mapping parameters and elapsed time to derive the expected RMS and NEFD from the ITC;
• Write results to log file.

Since a custom makemap config file can be given (see below), it is possible to use this recipe as your main data reduction recipe (if the standard processing is sufficient), to avoid re-reducing data. Note that the recipe does not yet do any coadding, so currently output maps will have to be combined with MOSAIC_JCMT_IMAGES.

To use it, run it like any other pipeline recipe. Initialize the pipeline in the usual way:

% oracdr_scuba2_XXX -cwd

Add the names of the relevant files to a text file (say, myfiles.lis), and run the pipeline with:

% oracdr -log sf -nodisplay -loop file -files myfiles.lis REDUCE_SCAN_CHECKRMS

All of the recipe parameters supported by the default REDUCE_SCAN recipe can be given to this recipe, the most important of which is likely to be MAKEMAP_CONFIG for specifying an alternative config type for makemap. If not given, the default config file (dimmconfig.lis) is used. Note that the same config file will be used for all data so be sure to use an appropriate config.

Log file format

The log file, called log.checkrms, produced by each method contain the following entries for each map (observation):

• UT - UT date including day fraction
• Source - object name, upper case with spaces removed
• Obs - observation number
• FILTER - filter (wavelength)
• telapsed - elapsed time of observation (sec)
• texp - mean exposure time, derived from EXP_TIME NDF component (sec)
• trans - mean line-of-sight transmission
• nep_av - mean NEP for current observation (W/sqrt(Hz))
• nep_av_err - uncertainty in nep_av (W/sqrt(Hz))
• rms_nep - RMS derived from nefd_nep (mJy/beam)
• nefd_nep - NEFD derived from nep_av and scaled elapsed time (mJy sqrt(sec))
• rms_map - RMS noise in map, obtained from median of error component (mJy/beam)
• nefd_map - NEFD derived from combination of variance and exposure time images (mJy sqrt(sec))
• rms_itc - RMS noise estimated by SCUBA-2 integration time calculator (ITC) (mJy/beam)
• nefd_itc - NEFD derived from rms_itc and scaled elapsed time (mJy sqrt(sec))
• itc_obstype - observation type used by the ITC to derived RMS
• rms_ratio - ratio of rms_map to rms_itc
• El - mean elevation in degrees
• CSO - mean zenith optical depth at 225 GHz, derived from WVM
• Tau - mean zenith optical depth at the current wavelength (given by FILTER above)
• Radius - radius in arcsec of image used in analysis
• pixscale - pixel scale in arcsec
• f - ITC f parameter
• project - project ID

Note that ITC-derived values will be undefined if the PONG map size is non-standard.

So what should I look for?

The log file contains all the results of the analysis. It's worth comparing RMS_ITC with the number you determined at the time your proposal was written. It's important to be aware that the PICARD and ORAC-DR recipes estimate values from the ITC using the actual observed airmass/opacity, rather than an average obtained from the source declination. This will undoubtedly lead to differences in the ITC-derived parameters.

Perhaps the most important number is the rms_ratio which is the ratio of the RMS_MAP to the RMS_ITC. Ideally this value should be 1, though in practice that's rarely the case. Sometimes it's down to the region used for calculating the RMS being too small (in this case, increaseIf the results are very far from the expected values then contact your Friend of the Project to see if it's possible to work out why.

SUMMIT pipeline

While not yet released, the SUMMIT pipeline will also perform its own CHECK_RMS analysis and will write a corresponding log.checkrms with the same format as for the other methods.

A few things to keep in mind when this is released:

• The NEP and derived values will be undefined (NaN) in the log file, as the SUMMIT pipeline does not have access to the QL log file data.
• An entry is written each time a new image is created, which may yield multiple entries for a given observation. However, each entry should be self-consistent.

This blog entry will be updated with the relevant documentation when the SUMMIT pipeline starts to produce these data.

Updates to SCUBA2_CHECK_CAL: 2-component fitting

Happy New Year to everyone reading this (or Hau‛oli Makahiki Hou as we say in Hawai‛i), and welcome to 2014.

The PICARD recipe SCUBA2_CHECK_CAL has been updated to use a two-component gaussian profile to fit the beam when determining an FCF from a calibration observation, as long as the signal-to-noise ratio exceeds 100. Previously, a single-component fit was used which was not constrained to be a gaussian. (If the signal-to-noise ratio is not high enough it will fall back on the old behavior.) This two-component fit is based on the FWHM and relative amplitudes of the two components derived in Dempsey et al. (2013).

This change has no effect on the ARCSEC FCFs, but results in a small (~1%) but consistent reduction in the BEAM FCFs. The effect should be small enough to be negligible, and with this change we now have slightly higher confidence in the resulting FCFs than with the old one-component fits. BEAMMATCH FCFs are likewise unaffected by the change, as testing showed they were best fit with profiles not forced to be gaussian.

You can set the behavior of SCUB2_CHECK_CAL manually using the FIT_GAUSSIAN parameter in your config file. The default value of 2 uses a two-component gaussian fit (when S/N > 100), a value of 1 uses a one-component gaussian fit, and a value of 0 recovers the old behavior of a one-component fit not constrained to be a gaussian.

At the moment, to use this feature you will need to update your Starlink to the current development version. Linux users can simply rsync Starlink from JAC following the instructions at http://starlink.jach.hawaii.edu/starlink/rsyncStarlink

Inclusion of CSO Tau Fits for Determining FCFs

During the first few months of 2013 the WVM at JCMT had several periods of instability where it was unreliable for determining the value of tau. Because of this, members of the science team collaborated to produce a collection of smooth polynomial fits to the tau values from the 225-GHz tau meter at the CSO for the affected nights, which can be used to perform extinction correction in place of the WVM tau values.

In the latest version of starlink (which you can get by rsyncing from /stardev at the JAC) makemap will now first check to see if the date of the observation is one of the (somewhat sporadic) dates when the WVM was unstable. This occurs as long as the ext.tausrc parameter is set to "auto," which it is by default. If the date is one of the affected dates, makemap will look for an available fit to the CSO data in the file specified by the parameter ext.csofit, which is set up by default to refer to the collection of CSO fits produced at the JAC. If makemap cannot find a fit for an observation in the specified path it will print a warning that no fit or WVM data is available and refuse to reduce the observation, though this shouldn't ever happen in normal operation.

If you have observations taken between January 19 and May 14 of this year, using the latest version of starlink rsynced from /stardev at the JAC will help ensure that you get the best extinction data available.

The graph below shows a comparison between the FCFs derived from WVM data and those derived using the new CSO tau fits over the period of January-May 2013. The blue diamonds represent FCFs from the WVM, and the tan circles are FCFs from the CSO tau fits.

New OMP features for projects.

The OMP has recently been updated to display several additional pieces of information on the project pages.

Tau graph

The first thing you will notice is a new graph of tau over the course of the night. Previously, the pages had such a plot that was generated each time the page was loaded. This took a while, and the plot had automatically generated limits on the y-axis which made comparing graphs between two night effectively impossible.

The new graph is made with consistent x- and y-limits that allows for easy comparison between nights. The various weather grades are also delineated, and the value of tau from the CSO WVM is plotted as well.

The new plot of tau vs. time. Click an any picture for a larger view.

The thicker gray horizontal lines running across the graph mark the weather grade boundaries, and the shaded area shows the night hours in Hawaii. The time in UTC is plotted along the bottom axis, with the time in HST along the top.

Pie chart

The pie chart of time spent in each grade.

There is also a pie chart that sums up the amount of time spent in each weather grade throughout the night. The grades are color-coded from green to red to give a quick visual assessment of how good a given night was (green good; red bad). Usually a given night will fall predominately within just one or two grades.


ACSIS standards table

For nights on which data was taken using ACSIS, there will be an HTML table in text form like the example below with some information relating to the ACSIS standards observed.

Obs #	Time	Integ. Int.	Peak Int.
16	19:36:32	251.33	8.82
19	19:50:31	4775.56	6.88
27	20:59:30	4157.76	5.96
36	22:22:36	3520.79	5.08

SCUBA-2 calibrations table

Similarly, for nights where SCUBA-2 was used there will be an HTML table like the example below listing the FCFs for each observation of a calibration object.

20120822			FCFasec				FCFpeak
Obs #	Time (UT)	Source	850µm	err	450µm	err	850µm	err	450µm	err
8	05:41:27	CRL2688	2.41	0.01	4.64	0.03	572.2	1.5	542.9	5.7
37	10:24:23	CRL2688	2.30	0.01	4.66	0.02	513.5	1.2	466.0	3.7
68	17:26:48	CRL618	2.23	0.01	4.28	0.04	487.8	1.4	436.8	4.8

Along with the table on SCUBA-2 nights there will three additional graphs, which will be detailed below. Each of these graphs has two sub-plots, the top one being the 450 micron data and the bottom one being the 850 micron.

NEPs vs observation number graph

This graph shows the min, max, and mean for each of the eight subarrays that make up SCUBA-2, plotted by observation number. The mean is the colored line, and the shaded areas are the filled-in areas between the min and max. By following the lines you can tell which arrays changed significantly, and by watching the shaded areas you can get a feel for the spread in the NEPs. The vertical scale is fixed, and is the same as the scale for the following plot.

NEPs vs time graph

This graph, like the previous one, shows the NEPs, but plots them against time rather than observation number. Each of the eight subarrays is present once again (with the same colors as in the preceding plot). The time (in UTC) is marked along the bottom of the plot, and the observations are marked by number along the top, with vertical lines that descend to help mark when each particular observation began.

NEFDs vs time graph

The final plot is a plot of the NEFDs vs. time. The number of points depends on the night, this particular night only had a few. Like the two previous plots, the vertical axis is fixed to make comparisons between nights easier.

Hopefully these new features will prove useful to users of the OMP, and additional updates or improvements may be forthcoming in the future. Feedback on the new features is welcome.

Checking SCUBA-2 calibrator data

A new PICARD recipe has been added called SCUBA2_CHECK_CAL which can be used to assess the quality of the calibration by calculating flux conversion factors (FCFs) and estimating the beam size to compare with standard values.

It can be run on any processed observations of standard sources (except, obviously, focus data). The input files must be uncalibrated, either the output from running makemap by hand or from the pipeline and then processed with the PICARD recipe UNCALIBRATE_SCUBA2_DATA.

Run it as any other PICARD recipe:

% picard -log sf -recpars myparams.ini SCUBA2_CHECK_CAL calfiles*.sdf

The recipe trims the image, optionally removing a background before it fits the source to estimate the beam.

The next step is calculating the FCFs which are reported to the screen. The map is then calibrated using either a standard FCF or one of those just calculated (controlled by a recipe parameter). The noise in the map is calculated from this calibrated image.

The recipe writes a log file called log.checkcal which contains various parameters. The values of interest include the total flux and uncertainty derived from the uncalibrated (input) map, the noise in the calibrated version of that map, the FCFs and associated uncertainties and the beam FWHM. In the future there will also be an estimate of the error beam. The log file can be read into Topcat so values can be plotted as a function of time, elevation or any other parameter (if there are enough of them).

The total flux is measured using aperture photometry with an aperture radius of 30 arcsec. The aperture size may be adjusted using a recipe parameter, but be aware that comparison with the standard ARCSEC FCF will require correcting for the different aperture. See the SCUBA-2 calibration paper for further details.

By default the recipe can only estimate FCFs from known calibrators. However, it also be used to estimate the FCF from non-standard sources if the flux is known to the user. Recipe parameters exist to allow this to be specified (on a per-source basis if desired). Note that the sources should still be unresolved, point sources.

The FCF calculations can be compared with the standard values (given in the calibration paper) to determine whether or not there may be particular issues with the calibration on a given night. The BEAM FCF is more sensitive to focus than the ARCSEC value, so if the former is out of spec, while the latter is not, it could indicate that the telescope was not as well focussed as it could have been. If both are far out of spec (indicated by red text in the recipe output) then it may suggest that a non-standard FCF is in order. It is also important to look at the trend in the FCFs over the night, and not just a single value.

The main recipe parameters of interest are:

• APERTURE_RADIUS - radius of photometry aperture in arcsec. Default is 30.
• USEFCF - a flag to denote whether the derived FCF should be used to calibrate the data. Default is 0 (use the standard FCF).
• FLUX_850 - total flux of a point source at 850 um. The corresponding 450 um parameter is FLUX_450. Fluxes may be specified for multiple sources by adding the source name (upper case with spaces removed) to the parameter - e.g., FLUX_850.MYSOURCE1, FLUX_850.MYSOURCE2 etc.
• REMOVE_BACKGROUND - a flag to denote whether or not a background should be removed. Default is 0 (do not remove a background). If true, then there are other parameters available to control the background fitting and removal process.

The full list of recipe parameters is described in SUN/265 and will be available in the upcoming Starlink release.

Checking the noise properties of SCUBA-2 science data

Ever wondered if there was an easy way to assess the noise properties of the data going in to your SCUBA-2 maps? Well now there is. A new recipe called ASSESS_DATA_NOISE has been added to the SCUBA-2 pipeline and it performs standard quality-assurance analysis (as done by the quick-look pipeline at the telescope) to give you an idea of the noise properties of the sciences data.

Running ASSESS_DATA_NOISE

To run it, first put all the files you wish to process into a text file (called myfiles.lis in the example below). Create a recipe parameter file if necessary (myparams.ini below). Then initialize the pipeline as usual (not necessary if done previously):

% oracdr_scuba2_XXX -cwd

where XXX is the wavelength of choice. Finally, run the pipeline, adding the name of the recipe at the end of the command line:

% oracdr -loop file -files myfiles.lis -recpars myparams.ini \
  -log sf ASSESS_DATA_NOISE

The recipe displays its results in a KAPPA window, showing the focal-plane layout with the noise for each bolometer in the top left corner, a histogram of the noise in the lower left, and the focal-plane image of the noise equivalent power (NEP) in the lower right. (If you want to turn off this display, add -nodisplay to the ORAC-DR command line above.)

By default the results are calculated for each subscan (i.e., each 30-second file) separately. As described below, options exist to just use the first subscan or to calculate the noise properties for the entire time stream.

The focal-plane mosaics for each subscan are stacked to create a 3-d cube of noise as a function of time (the third axis is MJD). The corresponding NEP images are written into this file under the .more.smurf.nep extension. This file has the suffix _scinoistack.

The results are written out to a log file called log.scinoiseXXX (where XXX is the wavelength as above), which can be loaded into TOPCAT to display the number of bolometers, noise or weighted NEP as a function of time or subscan number for each subarray.

The recipe also writes out a file called log.bolonoise, which contains more results (with the results for each subarray on a separate line), plus index.noise and index.nep. From these files it will be possible to determine which files failed the QA tests, and thus can be excluded from the input list to the map maker.

Of course, the down side to removing data is that it will not be possible to achieve the intended noise level, and may reduce sensitivity to extended structure (as the time series may be broken into smaller chunks). However, it should be easy to determine whether or not using some fraction of the best data gets as close as could be expected.

Customizing

Like most recipes, a number of recipe parameters are available:

• NOISE_CONFIG - the config file you wish to use when calculating the noise. It is recommended that you use the same as the one used when making maps. If not given, the standard noise config file is used.
• NOISE_CALC - may be "quick", "each" or "full" which will use the first thirty-seconds of data (for a quick check), or for every thirty second data file, or use the entire time stream to calculate the overall noise properties. The default is "each".
• FREQRANGE - frequency range over which to calculate the noise. Default is 2-10 Hz.
• FREQLO - lower frequency at which to calculate the background. Default is 0.5 Hz.

See the SMURF documentation for the FREQRANGE and FREQLO before adjusting them (SUN/258).

This is just the first release and its usefulness and output data relies on user testing and input. Please try it out and let me know.

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.

SCUBA-2 reference publications webpage

Please see the following page for appropriate references for the SCUBA2 instrument, calibration and data reduction at:


http://www.jach.hawaii.edu/JCMT/continuum/scuba2/scuba2_references.html

This includes an arxiv link to the Dempsey et al 2012 SPIE paper on SCUBA-2 commissioning which can be used as the reference for SCUBA-2 calibration.
http://arxiv.org/abs/1208.4622

A trilogy of up-to-date SCUBA-2 papers  on the instrument, data-reduction and calibration are currently ready for submission. Please check back to the above link periodically. The new links will be added when they are in press. 

Updates to FCFs and extinction correction

Over the past couple months the flux conversion factors and extinction corrections have been updated and we are now very happy with the answers. The main issue from the user perspective is that the FCF you apply critically depends on which version of SMURF was used to generate the map. At the time of writing this data products downloaded from CADC use an older extinction correction than the value you will find in place for the current SMURF. We plan to update CADC processing shortly but reprocessing the historical archive will take some time.

We have set up a web page at JAC listing the parameters that should be used and instructions on how to determine which version of the software was used to generate your map:

http://www.jach.hawaii.edu/JCMT/continuum/scuba2/scuba2_relations.html

Automatically adding fake sources to real data, part 2

In a previous entry, I described how to go use the new ORAC-DR recipe REDUCE_SCAN_FAKEMAP to automatically add fake sources to the raw time series in order to get a handle on the transfer function of the map-maker. That's all very well, but how about starting with something simpler, such as a Gaussian? Now you can add a Gaussian of given amplitude and FWHM without first creating an image to feed in to the pipeline. As before, the Gaussian properties are defined using recipe parameters:

[REDUCE_SCAN_FAKEMAP]
FAKEMAP_FWHM = 15
FAKEMAP_SCALE = 0.5
FAKEMAP_OFFSET = dx,dy

The FWHM is in arcsec and is the only mandatory parameter (the SCALE and OFFSET are optional). In this case, the FAKEMAP_SCALE parameter is the amplitude of the Gaussian in Jy/beam. The Gaussian can be positioned anywhere in the map by giving an offset (RA, Dec in arcsec).

The pipeline can help you out further by using a couple of useful defaults. If the world "beam" is given as the FWHM, then the pipeline will use the default beam size appropriate for the current wavelength (14" at 850 um, 8" at 450 um). If the amplitude is not specified then (somewhat arbitrary) defaults of 1 and 4 Jy/beam are assumed (at 850 and 450 um respectively).

Note that the FWHM must be greater than half the pixel scale (the pipeline will set it to that value if the given value is less), but in principle there is no upper limit. In practice the upper limit is defined by the array footprint (about 2.5 arcmin for S2SRO data): anything much larger than that will be removed by makemap as part of fitting the common mode.

As usual, update your Starlink installation (with rsync) to try it out.

Automatically adding fake sources to real data

Further to Ed's post on adding fake sources to SCUBA-2 maps, I've added some new functionality to the pipeline which allows this process to be automated.

There is a new recipe called REDUCE_SCAN_FAKEMAP which should be specified on the command-line when running the pipeline. (Normally the recipe is chosen automatically.) This recipe behaves in exactly the same way as the standard recipe, and just passes the given `fake' map (adjusted and scaled as necessary) to the map-maker.

The recipe has a small number of parameters which can be given in addition to those already accepted by the pipeline (e.g. an alternative makemap configuration file). These parameters are all prefixed with "FAKEMAP_":

[REDUCE_SCAN_FAKEMAP]
FAKEMAP_MAP = mymap.sdf
FAKEMAP_SCALE = 0.5
FAKEMAP_REGRID = 1/0
FAKEMAP_OFFSET = dx,dy

The name of the base input map is mandatory (otherwise the pipeline exits immediately). The remaining parameters are optional. The SCALE parameter is the same as the "fakescale" parameter used by the map-maker and defaults to 1 if not given. Two further parameters control how the input map is adjusted before adding to the raw timeseries data.

The OFFSET parameter is a pair of RA,Dec offsets in arcsec which allow the input map to be shifted so that it doesn't coincide with a source in the SCUBA-2 data. The default is not to shift the input map. The offsets are applied to the source: i.e. a shift of 60,60 puts a source originally at (02:30:00.0, 00:00:00.0) at (02:30:04.0, +00:01:00.0) in the output map.

The REGRID parameter is a flag which, if true, regrids the input map to match the pixel coordinates of the output map (i.e. same pixel bounds and pixel size). The default is false, unless an OFFSET is given: in this case the input map must be aligned (regridded) to match the output map so that the pixel bounds are correct. Note that the original input file is not modified: the pipeline makes a copy and operates on that.

Here's an example of calling this recipe on the command line (processing a list of files contained in mydata.lis):

% oracdr -loop file -files mydata.lis -recpars myrecpars.ini REDUCE_SCAN_FAKEMAP

This new capability is available now with an update of your Starlink installation. Please try it out and let me know how well it worked (or didn't!), or send a message to the scuba2dr mailing list.

SCUBA-2 at the ADASS conference

At the start of November I attended the ADASS conference in Boston and presented a poster paper on SCUBA-2 data processing. Readers of this blog probably won't find anything new in here but just in case, the conference paper itself is now on astro/ph at http://arxiv.org/abs/1011.5876. Thanks to Dave Nutter for supplying the Orion B image.

Map-based Despiking, Flag maps, Sub-scan maps, and Filters based on requested angular scales

This post summarizes a few odds and ends that were implemented over the last couple of weeks. The following examples use the 850um scan of Orion (20100216 observation 22) that was provided with the cookbook at http://www.starlink.ac.uk/extra /sc19/sc19_extended_galactic.tgz

Map-based Despiking

First, a new method for despiking the data has been added. Instead of searching throughout the time series for individual samples that are outliers compared to their neighbours (which can accidentally flag large astronomical sources as spikes in some cases), this new method looks at the weighted scatter of samples that land in each map pixel, and removed outliers. This feature is now turned on by default in dimmconfig.lis:

ast.mapspike = 10

which flags all samples that are > 10 standard deviations away from the other samples that contribute to the flux in the same map pixel, after each iteration. In fact, it works so well that we are using it instead of the time-domain de-spiking in most cases.

We first reduce the Orion data using the default dimmconfig_bright_extended.lis,and the ERROR component of the map looks like this:



we then re-run it, but turning off the default map-based despiking, and turning on time-based despiking:

ast.mapspike = 0

spikethresh = 10

spikebox = 50

and the ERROR component of the map now looks like this:



This latter image clearly shows the locations of large spikes that made it through the time-based despiking. Also note that their locations are not correlated particularly with the brightest regions of emission that are obvious at the lower-right (the large scatter in the ERROR component is a reflection of other systematic errors, such as in the flatfield, that become more prominent in regions of strong astronomical signals).

Flagmaps and Shortmaps

We have added two new types of map extensions to help identify problems while map-making. "flagmaps" are images that count the number of samples with QUALITY bits that match a given quality mask, as a function of their position on the sky. This is useful to see, for example, whether DC step correction and/or de-spiking is affected by bright sources. For example, setting

flagmap = (BADBOL,SPIKE,DCJUMP)

will produce a set of .more.smurf.flagmaps extensions corresponding to each continuous chunk of data indicating where all spikes and DC jumps were located. This is the resulting flagmap for Orion:



The "U" shaped patterns primarily indicate the locations of DC steps that were correlated across the array, and a number of the other spots are simply locations of spikes identified by the map-based despiker. To gain further insight, maps can be re-run, setting flagmap to a single quality bit each time.

There is also a modified form of the "shortmaps" extensions. Instead of specifying the amount of time for each map explicitly, you can now set

shortmap = -1

which will produce a new .more.smurf.shortmaps extension each time the telescope completes a full, cross-linked pass through the scan pattern. In our Orion example, we show the full map, followed by the four sub-maps (.more.smurf.shortmaps.ch00sh00000 through .more.smurf.shortmaps.ch00sh00003) clearly showing how the position angle was rotated for each fully-crosslinked pass across the field:



Choosing frequency-domain filters based on angular scale

Finally, we have implemented an alternative method for specifying frequency edges. In the past, the frequency edge would be given explicitly, e.g.

filt_edgehigh = 0.3

flt.filt_edgehigh = 0.3

to remove all power in the bolometer time-series below 0.3 Hz using the pre-processing and iterative filters respectively (this is the approximate 1/f knee). However, in order to evaluate the impact of such filters on real astronomical structure, one would have to divide the scan rate by this number to convert it to an angular scale. The requested scan rate is recorded in the "SCAN_VEL" FITS header, which for this Orion map was 250 arcsec/sec. Therefore a high-pass filter explicitly removes any information in the map on scales > (250/0.3) = 833 arcsec.

This calculation is cumbersome, and also relies on the requested "SCAN_VEL" being correct. We have therefore added an alternative method for specifying filters based on angular scale, e.g.

filt_edge_largescale = 600.

450.filt_edge_smallscale = 2.

850.filt_edge_smallscale = 4.

In this case a high-pass filter will be selected such as to suppress angular scales larger that 600 arcsec, and low-pass filtered will be used to remove scales smaller than roughly the FWHM/4 of the PSF in each band. The calculation is based on the actual mean slew velocity is measured from the pointing data rather than relying on "SCAN_VEL". For example, in the above Orion maps, the default iterative filter is set to

850.flt.filt_edge_largescale=300

and if the map-maker is run with msg_filter=verbose, when the FLT model is evaluated you will see

flt

smf_scale2freq: Based on a slew speed of 243.1 arcsec/sec, setting:

smf_scale2freq: FILT_EDGEHIGH = 0.810 Hz (less than 300.0 arcsec scales)


Even though the 1/f noise varies from observation to observation (and is at a higher frequency at 850um for SRO data), it may be useful to specify the same angular scale irrespective of scan rate (such as ensuring that the attenuation of the filtering is the same for calibrators as for science data, even if they were scanned at different rates).

Related to all of this, we have replaced the "flagstat" parameter, for flagging data when the telescope was slewing too slowly, with

flagslow = 30

450.flagfast = 600

850.flagfast = 980

These default values include only data where the telescope was slewing at least 30 arcsec/sec, and no more than 600 arcsec/sec at 450um, and 980 arcsec/sec at 850um. The lower limit is set by an approximate 1/f knee of 1 Hz to ensure that at scales of order the 850um PSF will be useful in the data. Similarly, the upper speed limits are chosen to avoid smearing given the PSF sizes in each band and the 200 Hz sample rate.

Suppressing large-scale noise with zero-masks

Over the last couple of weeks we have been implementing and testing two new masking features that significantly improve the large-scale noise properties. As mentioned in an earlier post the common-mode signal (and also low-frequencies removed using a high-pass filter in the FLT model) are degenerate with angular scales larger than the array footprint on the sky.

In that earlier post, I described a method of specifying a zero boundary condition in which the regions of the maps with low hits (the edges) are set to zero after each iteration. While this helped significantly with smaller maps (particularly pointing observations), it wasn't much help for significantly larger maps (i.e., covering much more area than the array footprint). This would lead to "saddle" shapes all over the map as the number of iterations was increased. We have therefore added two new masking options.

The first is to define a circular masked region:

ast.zero_circle = (LON,LAT,RADIUS)
or
ast.zero_circle = (RADIUS)

In this case, the portion of the map beyond a circular region centered over LON,LAT with the given radius (all quantities in degrees) is constrained to zero. Alternatively, only a single value may be supplied, in which case it is interpreted as the radius, and the LON,LAT are taken to be the pointing centre of the map (usually the compact source of interest for which this option is generally useful).

The second method automatically determines the masked region based on a SNR threshold.

ast.zero_snr = SNR

In this case, after each iteration of the map, all map pixels above the threshold SNR are included in the "source mask", and map pixels outside this mask are constrained to zero. As bright structures "grow" with respect to the neighbouring zero-constraiend map pixels, each iterations of the source mask will also grow to include more pixels.

For both of these options, the ast.zero_notlast flag applies: for the final iteration the zero constraint will not be applied if this is set.

So from a numerical point of view, what does this constraint do in the regions that were constrained to zero until the last iteration? In its simplest form the map-maker is trying to solve for several signal components in the bolometer time series:

Signal = ASTronomical + COMmon + FiLTered

where AST is the signal produced by astronomical sources, COM is the average signal seen by all of the bolometers, and FLT is generally residual low-frequency noise. By setting AST to zero at any location in the map, the only other places that signal can go are COM and FLT. COM only contains the average signal seen by all of the detectors, and FLT only contains low-frequency signal components. Therefore this operation, perhaps surprisingly, will not remove any compact structures in the zero-padded regions; the effect seems to be roughly equivalent to setting the mean value on array-footprint-sized patches within those regions to zero.

The updated dimmconfig_bright_compact.lis is a good default configuration for small/bright sources; in this case a circular mask is used with a radius of 60 arcsec. Similarly, there is an updated dimmconfig_bright_extended.lis that masks based on a SNR threshold of 5-sigma. We illustrate the effect of the SNR masking using the short observation of Orion mentioned in the previous blog post, s4a20100219_00039.

First, a map produced without any zero constraints and 20 iterations, clearly showing the large-scale degeneracy that results in a strong gradient in the map:



Here is the same reduction, but setting ast.zero_snr=5 and ast.zero_notlast=1:



and the mask indicating which pixels are constrained to zero until the final iteration (white) -- this is the QUALITY component of the final image:



Note that 20 iterations was chosen in this example to illustrate the improved convergence properties when using the zero mask. In fact, only about 8 iterations are needed in this particular case. A better automatic convergence test is still an outstanding issue for SMURF, and at present we recommend checking the maps after each iteration visually.

SCUBA-2 Pointings in the science archive

This is just a quick note to explain what has been happening to pointing observations in the science archive. Pointings are public data so the reduced images of pointing sources account for a significant fraction of publicly available SCUBA-2 reduced data. We have realised that these data are problematic because we were running the standard pointing recipe at CADC which filters out large scale structure and clips the image in order to make it easier to detect the pointing source. This is not obvious to the casual user of the archive and has recently caused some confusion.

The image above (a 20 second 450 micron observation from 2010-02-19 #39) is the standard product for a pointing of OMC-1. As is obvious, a lot of extended emission has disappeared. This week we are going to start reprocessing pointing observations using the standard recipes so that the archive contains science products rather than products that aid telescope observing. As an example the image below is the same observation using the standard recipe and it is much more representative of what SCUBA-2 is capable of:

Spatial noise vs noise maps

I've been comparing the spatial noise in the two SRO SASSY maps (0.5 degrees by 0.5 degrees each) to the noise maps produced by the pipeline and matched-filter picard script. I've been told that this is a known issue, but that I should also share my findings. To do this I made histograms of signal to noise pixel values within the two maps that were processed using the matched-filter picard script. Since the maps are mostly empty and the matched-filter script should remove most of the artificial large scale structure introduced by the mapmaker algorithm, the spatial noise should agree well with the noise maps produced by the pipeline.

850μm maps: For the two SASSY fields (one with 29mJy rms and another with 47mJy rms), the signal to noise map underestimates the noise in the map by about 17% (± about 2% since I'm doing this by eye).

Larger image here

450μm maps (we don't expect to see any astronomical signal from these): For the NGC 2559 field with 1700mJy rms, the noise map underestimates the noise in the map by about 29%. For the cirrus field with 2400mJy rms, the noise in the map is overestimated by about 13%.

Larger image here

Larger image here

Hopefully this is helpful for others who are trying to figure out the amount of noise in their maps and are wondering about the accuracy of the noise maps produced by the pipeline.

S2SRO data processing notes

A bunch of S2SRO people are asking whether they can just download the JSA products and go, without doing any of the data processing themselves. The answer to that is obviously yes, as long as you have convinced yourself that the JSA processing (which is also the default ORAC-DR processing) in combination with your data characteristics will support your scientific conclusions.

In aid of this, Ed Chapin has provided the following note summarising what we are doing to your data (thanks Ed).

The S2SRO data maps

At the time the data were acquired, SCUBA-2 was operating with single subarrays at each of 450um and 850um (the full complement is four subarrays in each band). In addition, there were a number of challenging data artifacts to overcome, including glitches in the bolometer time-series, and severe low-frequency noise caused, primarily, by oscillations in the SCUBA-2 fridge and magnetic field pickup (for further details, see SC/19 Appendix A). Despite these problems, it has been possible to produce science-grade maps for most of the S2SRO projects using one of two standard sets of configuration parameters for the map-making software (see the SMURF manual).

The first set is designed to produce maps of extremely faint fields containing point sources, using strong high-pass filtering to remove most of the low-frequency noise in the data (see dimmconfig_blank_field.lis in SC/19). It should be emphasized that any structures larger than the point spread function (PSF) are removed by the processing. In addition, these maps are run through a matched-filter to assist with the identification of point-sources (see further description below in "Noise Estimates").

The second set of parameters is used to reduced all of the remaining data, including maps of extended and/or bright objects (see dimmconfig.lis in SC19). The filtering is less ambitious in this case, and it is performed iteratively, which reduces ringing around bright sources, but may produce large-scale patchy or saddle-shaped structures in the maps.

For both types of processing, a common-mode signal has also been removed iteratively from the data --- the average signal recorded by all of the detectors at each instant in time. This processing, again, helps to suppress low-frequency noise, but removes all structure in the maps that are larger than the array footprint (approximately 2' x 2').


Calibration

Extinction correction has been performed on all of the data using the line-of-sight water vapour radiometer (WVM) on the JCMT. It is not currently believed that the calibration uncertainties are dominated by noise in these measurements.

Flux conversion factors (FCFs) have also been applied to all of the data. These factors have been established using regular observations of point-like flux calibrators. While there have been large variations observed in these factors, it is not presently understood whether these are real systematic variations (i.e. as a function of time of day, surface quality etc.) or simply reflect measurement uncertainties (see SC/19 Appendix B.1). We therefore apply a single FCF at each wavelength within the JSA pipeline, and assuming the errors are not systematic, estimate the total calibration uncertainties to be about 20% in a random observation.

Unless there are reasons to question the calibration factors for a particular observation, we believe that the current values are the best that can be applied.


Noise Estimates

The map-making process provides an estimate of the noise in each pixel (the VARIANCE component). Its value is calculated by:

• (i) estimating the weighted sample variance for the N bolometer samples that land in the pixel (which are averaged together to estimate the mean flux in that pixel; weighted by the inverse of the noise in the respective bolometers from which they came).
  
  
• (ii) dividing the result of (i) by N to calculate the variance on the mean (i.e. the variance on the estimated flux).
  

On small scales, this noise estimate appears to be accurate. For example, in maps produced using the faint-field processing, which are run through a matched-filter to correlate samples on the scale of the PSF, the noise is well-behaved (see Section 6.1 and Figs.13-14 in SC/19).
However, this noise map does not accurately characterize the uncertainty in aperture-photometry measurements on scales larger than the PSF, especially in the case of bright-field processing, since there is low-frequency noise that is correlated between different pixels. In these cases, the noise does not drop as sqrt(M) as expected, where M is the number of pixels in the aperture. For example, if a cluster of pixels all happen to reside on a local postive noise patch, the uncertainty in the photometry of that cluster is dominated by the error in the baseline of the patch, not the small-scale noise estimated for each pixel. It is therefore necessary to estimate the noise for larger-scale measurements by placing measurement apertures at random locations to estimate the real (and larger) uncertainties.

Finally, we re-iterate the fact that a common-mode signal has been removed from the data, so that in practice the SCUBA-2 maps only contain information on scales ranging from the PSF (6"--14") up to the size of the array (2').

SCUBA-2 products update - fresh batch

Just to let folks know, all the S2SRO data in the JSA has been reprocessed to pick up any observations marked as bad (and also to pick up some infrastructure changes I won't bore you about). Whether it is worth your while or not to re-download these products depends on how diligent you were about marking any bad observations since the last re-run a month or so ago.

Feel free to continue flagging bad observations as you find them - the data will be re-processed again in the future, or you can just drop us a line and ask for a reprocess of your particular project any time.

A few of the groups failed reduction - due to the aforementioned infrastructure improvements we are in a much better position to keep track of these failures, so we are working through the list to resolve the issues. If you are desperate for something you don't see in the archive, don't be afraid to ask to be bumped up the list.

By the way, a few of you have noticed that the JSA "house" products look better than what you are managing at home. Rest assured that there isn't some secret sauce that we are feeding the hamsters; maps will look better when processed on high-memory (16G+) nodes because they will have a higher effective integration time due to the way the mapmaker works. And of course the JSA processing always uses the most recent version of the software, which if you are not one of our bleeding edge Starlink rsync users, will net you a considerable improvement.

SMURF Update (August 5th 2010)

It's been a couple of months since the last SMURF news entry so I thought I'd bring people up to date.

• All the configuration files for the iterative map-maker have been tweaked and some have been renamed. For example the '_faint' config is now called "dimmconfig_blank_field.lis". We also have a new config file for bright calibrators called "dimmconfig_bright_compact.lis".
• The SMO (time series smoothing) model has been improved. A few bugs have been fixed and it's been parallelized and so is much faster. SMO is not enabled in any of the default configuration files.
• The size of the apodization and padding for the FFTs can now be calculated dynamically based on the requested filter. This is now the default behaviour. A new Fourier filter has also been written that does not require apodization (which may be important for very short maps) but is still being tested.
• Quality handling has been revamped inside the map-maker to allow us to report more than 8 different types of flagging. The report at the end of each iteration is now more compact and if you use the SHOWQUAL command to look at exported models you may see that the bit numbers assigned to a particular quality are no longer fixed. Additionally if more than 8 flags are used the exported model will combine related flags (for example PAD and APOD will be merged into ENDS).
• Very noisy bolometers will now be discarded before the iterative map-maker starts. This can help with convergence. See the "noiseclip" config parameter to adjust this.
• The map-maker now compares flatfield ramps taken at the start and end of each observation and disables bolometers that whose calibration has varied too much. This will not help data taken prior to 20100223 where flatfield ramps are not available.
• The step correction algorithm continues to be improved.
• SC2CLEAN will now report quality statistics in verbose mode.
• SC2FFT can now be given a bad bolometer mask.

The cookbook has also been updated and can be read online.