Pipelines and Archives

SCUBA-2 Calibration: REDUX.

2012-02-03T10:57:00.023-10:00

The short story:

The heater coupling factors have been adjusted to more realistic values. In practice, this does not change the performance of the instrument - however it does change the absolute value of the FCFs. These values were adjusted in the software in mid-December.

The WVM tau algorithm has been fixed and improved. This will not affect you directly: though the nightly plots now look extremely good and are officially used for weather band determination.

This has allowed new, and better calculation of the relation between the 225GHz tau derived from the WVM and the opacities at the two SCUBA-2 filter-bands. They are now as follows:

TAU_[850] = 4.6 * (TAU_[225] - 0.0043)
TAU_[450] = 26.0 * (TAU_[225] - 0.019)

The FCFs (flux conversion factors) have been derived for both wavelengths from an extensive reduction of calibrator sources observed over eight months of SCUBA-2 commissioning and science verification observations. They are as follows:

850um:
FCF_[arcsec] = 2.42 +/- 0.15 Jy/pW/arcsec**2
FCF_[peak] = 556 +/- 45 Jy/pW/beam
Beam area = 229 arcsec**2

450um:
FCF_[arcsec] = 6.06 +/- 0.32 Jy/pW/arcsec**2
FCF_[peak] = 606 +/- 55 Jy/pW/beam
Beam area = 97 arcsec**2

Reminder on how to calibrate your data:

Other posts discuss how best to reduce your data (and what recipes are needed). The latest software releases (since January 2012) all include extinction correction (with the relations above) and the changed coupling factors. If you reduced your data prior to this, you will need to reduce them again to account for these changes. Applying the FCFs reported here to old reductions of your data will be wrong.

The arcsec FCF: (when you want integrated fluxes)

The arcsec FCF is the factor by which you should multiply your map if you wish to use the calibrated map to do aperture photometry.

The peak FCF: (the FCF-formerly-known-as-beam):

This FCF is the number by which to multiply your map when you wish to measure absolute peak fluxes of discrete sources.

The whys and the wherefores (the gory details):

Reductions of the calibration observations:

Uranus and Mars were used as the primary calibrators for these results. In addition, CRL 618, CRL 2688 were also predominant secondary calibrators. All of the calibrators were reduced in January 2012 using the updated dimmconfig_bright_compact.lis. The improvements in the recipe (mostly in how it chooses to stop iterating) have resulted in extremely flat maps with nearly no evidence of the 'bowling' seen around strong sources during S2SRO reductions. To improve the accuracy of the peak-fitting and aperture photometry, the maps were reduced with 1 arcsecond pixels at both wavelengths.

Once the maps were reduced they were then analysed using the PICARD script SCUBA2_FCFNEFD. There have been some changes to this script: some few bugs have been fixed, the average FCF's have been adjusted, as well as the (now emprically derived) beam area, and a few reference fluxes have been adjusted. FCF_beamequiv has been removed entirely, and all calculations of integrated values are now done using AUTOPHOTOM, with a defined aperture and annulus for background subtraction.

Following extensive analysis of the most optimal parameters, all calibrations were reduced using a 60" diameter aperture (at both wavelengths) with an annulus between 90" and 120" from the source position.

Questions? Let's provide answers to a few we've already seen:

"These FCFs are very different to the old numbers quoted!"

Yes they are. The heater coupling factor change and the new tau relations play a significant part in this. But in addition, the large sample of observations has allowed for a much more accurate determination of the beam area. At 450um in particular, optical effects of the telescope show that the error beam is large and the beam is not gaussian. This results in an effective FWHM that is much broader than the 7.5" quoted previously (though that is the approximate FWHM of the fit to the centre of the beam) - it is more like 9.5", taking into account the error beam. Therefore, the measured (and fitted) peak is relatively lower, requiring a higher FCF to calibrate the peak flux in your data.

"There is a lot of scatter when I calculate the 'beam' or peak FCFs for my calibrators (particularly at 450um)"

No kidding. Peak values are obtained either from reading off the peak of the map (in gaia or by another method) or by fitting to the peak using beamfit (as is done in PICARD). The beam shape (particularly at 450um) can be extremely susceptible to changes in focus and atmospheric instability, amongst other things. The integrated value (FCF_arcsec) is more robust against such changes. If you are measuring a peak fit from a calibrator and see a strong deviation from the expected value things to check are:

- how 'focussed' does the image look? If you see distortion in the shape of a source that should be point-like, or distortion or 'shoulders' in the beam then it is likely that the peak value will be unreliable.

- was the observation taken early in the evening? Focus and atmospheric effects are known to be worst in the early evening hours and sometimes in the morning after sunrise. If you are looking at calibrators, try and look at ones taken later in the night and see if there is improvement.

A 'trap' has been set in PICARD to warn you if the attempted fit to the peak misses the actual peak value by more than 10%. Looking at the fit to the shape also helps in this instance. In any case, the quoted peak FCF value at the top of the post is derived from the arcsec FCF and the empirical beam area derived from nearly 500 observations at both wavelengths and this number has been shown to be robust.

"How stable are these FCFs? (read: do I need to reduce my own calibrators?)"

Very. The absolute errors at both wavelengths are within 5% and no significant trends have been seen in the last six months. Instrument performance is being monitored very closely and any deviations are likely to be noted specifically. However, we do not discourage you taking calibrators from the nights your data was taken and reducing them yourself - we appreciate the sanity checks! Another handy rule: if you do it to your data, do it to your calibrator. If you have specific methods you plan to use on your data, apply the same methods to your calibrator in order to ensure your calibration is correct. We are now happy to say though, that these FCFs look stable and correct, so using these numbers should provide you with well-calibrated data.

"What happened to FCF_beamequiv?"

FCF_beamequiv is a seductive, evil little value that tempted us to stray to the dark side, albeit temporarily. In essence it was created to use as a comparison to SCUBA performance, but should never have been used to actively calibrate SCUBA-2 data as it assumed a perfect gaussian beam. The statements above explain that this is patently untrue, especially at 450um. The beamequiv number was quoted previously, and incorrectly, as the true FCF, and it is largely the reason that the new (and correct) numbers seem so much larger. We have banished it from PICARD and it shall now be known as the FCF-that-shall-not-be-named.

Supplying an externally-generated zero mask

2011-11-16T12:15:00.007-10:00

In an earlier post I described zero-masking , in which the map is forced to a value of zero in certain locations to help with convergence (i.e. reduce saddles, negative bowls around sources etc.). The two methods that have been available until now are: defining circular regions beyond which the map is set to zero (good for known compact sources); and using a threshold on the local SNR (an automatic way to identify source and blank sky).

Since SMURF still handles non-continuous pieces of data independently (i.e. multiple maps of the same field), the SNR of an individual map may not be sufficiently large to produce a good mask of source / blank sky, whereas the combination of all data sets do. For those cases, and other examples where you may know, a priori, where to expect emission, a new facility has been added so that the user can define a zero mask externally and then provide it to SMURF.

In this illustrative example, I use one of the OMC-1 maps from 2010 with the s4a array: obs 18,22,23 on 2010216 and obs 27,29,33,34 on 20100218. First, we make a map using the bright_extended configuration (which uses SNR-based zero masking; I also use a lot of down-sampling and big map pixels to speed things up for this example):

makemap ^names_450.txt map_450_be pixsize=8 \
method=iterate \
config='"^/stardev/share/smurf/dimmconfig_bright_extended.lis"'

The resulting image looks reasonably good, but there are obvious negative bowls around the bright extended sources. However, the SNR of the final image is quite good, so I create a new mask based by thresholding the SNR map (4-sigma), and then smoothing it with a 5-arcsec FWHM Gaussian to spread out from where the sources are slightly:

makesnr map_450_be map_450_be_snr

thresh map_450_be_snr mask_450 thrlo=4 newlo=0 \
thrhi=4 newhi=1

gausmooth mask_450 mask_450_sm fwhm=5

thresh mask_450_sm zero_mask_450 thrlo=0.05 newlo=bad \
thrhi=0.05 newhi=1

I can then feed this new "zero_mask" back into makemap and re-run the data:

makemap ^names_450.txt map_450_zm pixsize=8 \
method=iterate \
config='"^/stardev/share/smurf/dimmconfig_bright_extended.lis,itermap=1,ast.zero_mask=1,ast.zero_snr=0"' \
ref=zero_mask_450

The key things to note are: (i) I turned off the SNR thresholding, ast.zero_snr, that is the default for the bright_extended configuration, and (ii) the new zero mask is provided as the REF image (also ensuring that the map and mask are on the same pixel grid), and to use it we set ast.zero_mask=1.

In the following images I show (from top-left): (i) an arm sticking out towards the north of OMC-1 using the default bright_extended reduction; (ii) the same region after using the updated zero-mask; (iii) the first map subtracted from the second map to illustrate the improved response to large-scale structure; (iv) the original combined zero mask based on the SNR of map pixels in individual scans; (v) the new zero mask based on thresholding/smoothing the combined SNR map:

Automatically adding fake sources to real data, part 2

2011-03-31T08:37:00.004-10:00

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.

Recipe Parameters and the Science Archive

2011-03-18T09:33:00.000-10:00

When processing jobs are submitted to the science archive ORAC-DR is configured to load a specialist recipe parameter file associated with the project. This can be used by PIs or surveys to tune processing to their liking. The idea is that PIs or survey teams send us their recipe parameter files and we then make sure that the processing jobs use them.

This week I've updated the recipe parameter system to allow tuning based on the object name as well as the recipe name. In some projects the objects are completely different so a single parameter file was not useful. Hopefully this change will encourage more people to send us parameter files.

To use the object based scheme simply append the object name to the recipe name in the file along with a colon. No spaces in the object name.

  [REDUCE_SCAN:M82]
  PAR1 = A
  PAR2 = B

If there is also a recipe-based entry those parameters will be merged in

  [REDUCE_SCAN]
  PAR1 = Aprime
  PAR3 = C

So in this example if the object name is "M82" PAR1, PAR2 and PAR3 will be set and PAR1 will have a value "A". If the object is no "M82" only PAR1 and PAR3 will be set and PAR1 will have value "Aprime".

This change is currently on the stardev rsync server, or if you have git installed you can update your own ORAC-DR distribution.

Automatically adding fake sources to real data

2011-01-18T08:18:00.002-10:00

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.

Correcting for pointing drifts with PICARD

2011-01-11T08:48:00.011-10:00

If you're lucky enough to have a bright source in your map, you may be able to improve the image fidelity by correcting for pointing offsets between observations. The pipeline now has this capability and it is applied automatically for calibrators, but not for most science targets. I've added a new PICARD recipe, called SCUBA2_REGISTER_IMAGES, which allows you to shift your images to line up at a common reference position before mosaicking. Here's how to do it.

First step, upgrade to the latest Starlink version via rsync.

Step two, create a recipe parameter file (say, params.ini) and add the following lines to it:


[SCUBA2_REGISTER_IMAGES]
REGISTER_IMAGES = 1
REGISTER_RA = HH:MM:SS.S
REGISTER_DEC = +DD:MM:SS.S

where "HH:MM:SS.S" and "+DD:MM:SS.S" are the RA and Dec of the source you wish to use.

Then run


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

and if all goes well, you will end up with a series of output files, one for every input file, which end in "_reg.sdf". These files can then be coadded using MOSAIC_JCMT_IMAGES.

The recipe works by fitting a 2-d profile to the brightest source within a 2x2 arcmin² box centred on the given reference position. The WCS is adjusted to place the reference position at the position reported by the fit.

However, note that this method will only work if the source you're using is present in every image, is compact and clearly resolved from other nearby sources and is the brightest source within the 2x2 arcmin² search box in every input image. It won't work for deep extragalactic maps, nor does it correct for any pointing drifts that occur within an observation. Correcting for those drifts is much more involved, but it is something we'll investigate over the coming months. It's been tested on HL Tau, OMC-1 and AFGL4029 in W5-East.

It's also possible to do the registering and coadding with a single call to MOSAIC_JCMT_IMAGES. In this case you'd want to add the following to a recipe parameter file (in addition to any other coadding parameters):


[MOSAIC_JCMT_IMAGES]
REGISTER_IMAGES = 1
REGISTER_RA = HH:MM:SS.S
REGISTER_DEC = +DD:MM:SS.S

Here's an example using some observations of HL Tau data taken in December 2009. The corrections were as much as 8.5 arcsec in RA (the Dec corrections were less, up to 3 arcsec). The improvements are modest at 850 um (of order 10% increase in peak flux and a few per cent reduction in beam size) but much more significant at 450 um. The peak flux after correcting for pointing drifts is more than 20% higher than before, and the FWHM has decreased by 15% (9.2" to 7.9"). Furthermore, the source elongation (ratio of major and minor axes) is reduced from 16% to less than 3%. Before and after images are shown below.

Left: HL Tau at 450 um, coadd of 9 separate observations from 20091213 and 20091214. Right: same data after registering each image before coadding. Both images are scaled between the same limits.

SCUBA-2 at the ADASS conference

2010-11-29T09:13:00.001-10:00

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.

Adding fake sources to real data

2010-10-27T11:20:00.011-10:00

Until now it has been difficult to characterize the ability of SMURF to reduce maps of known sources. Certainly a basic simulator has been available since before SCUBA-2 began to take real data, but no effort to date has been put into updating the noise model to accurately reflect the real instrument.

An alternative is to simply add signal from synthetic astronomical sources to real time-series. This feature was relatively easy to add and can be quite versatile. The user simply specifies an external 2-dimensional map using the new configuration parameter fakemap to makemap, containing the simulated sources. There are several things to be aware of when using this facility:

the map of fake sources must have the same pixel dimensions as the output map
the units are the same as the output map from makemap (normally pW if working from raw data)
a scale factor can be applied to the map using the fakescale configuration parameter

The following two examples illustrate how this feature might be used.

Response to a point source

In this example we create a beam map using a scan of Uranus (s4a20091214_00015), and add it to the time-series for a scan of a blank field (s4a20100313_00029) to establish the effect of map-making and matched-filtering on point sources. Both of these data sets are public (See SC/19 and this web page).

First, maps of the blank field and Uranus are both reduced:

makemap s4a20100313_00029\*.sdf blank450 method=iterate config=^$STARLINK_DIR/share/smurf/dimmconfig_blank_field.lis

makemap s4a20091214_00015\*.sdf uranus_450 method=iterate config=^$STARLINK_DIR/share/smurf/dimmconfig_bright_compact.lis

Then we cut the central region out of the Uranus map and copy it into a map (using KAPPA:NDFCOPY) with the same pixel dimensions as the blank450.sdf so that we can use it to provide fake signal:

ndfcopy 'uranus_450(3~31,1~31)' fakeblank450 like=blank450

We create a modified dimmconfig that instructs makemap to add in this additional source of signal, but applying a scale factor such that the amplitude of Uranus is 0.01 pW:

^/stardev/share/smurf/dimmconfig_blank_field.lis

fakemap = fakeblank450

fakescale = 0.0471865

We run makemap using this new config:

makemap s4a20100313_00029\*.sdf blank450_fake method=iterate config=^dimmconfig.lis

Finally, we produce match-filtered maps for both the original map, and the map with the fake source added in:

picard SCUBA2_MATCHED_FILTER blank450

picard SCUBA2_MATCHED_FILTER blank450_fake

The following image shows all four maps. The left column contains the raw output of makemap for the original map, and the map with the source added in. The right column shows the corresponding match-filtered images.

The input source had a known amplitude of 0.01 pW. In the makemap output, however, we can see that the filtering during map-making has reduced the amplitude by about 12% to 0.0088. The peak in the match-filtered image, however, is not attenuated quite as severely, at 0.0096. In addition to the peak response, this test also clearly shows the effective PSF for the map: there are faint negative streaks along the scan directions caused by the high-pass filtering. This example suggests that upward corrections to the FCF of order 5 to 10% should be applied to maps of points sources using dimmconfig_blank_field.lis (although this factor should be established on a case-by-case basis).

Response to extended structures

It is known that maps produced with SMURF do not have useful information on scales larger than roughly the array footprint due to the common-mode rejection, and high-pass filtering. Using the fakemap facility, we can now add sources with complicated structures to SCUBA-2 data to see this effect quantitatively.

In this example we use an 850um scan of Orion (s8d20100216_00022), and to it we add some diffuse structure from a publicly available map of Vela at 350um from the Balloon-borne Large-Aperture Submillimeter Telescope:

Similar to the previous example, we first reduce the map of Orion (using dimmconfig_bright_extended.lis),

and then we reduce it again using the BLAST image as a fakemap (positioning it and padding it to the correct size using KAPPA:NDFCOPY) and scaled to a similar signal range as the actual structure in Orion:

Comparing the upper-left region of this map where the Vela data were added with the original BLAST data, it is clear that all of the compact structures are faithfully identified, but the extended structure is not. Since we know what those underlying structures are, we can simply difference the input and output maps to see what we are missing in the SMURF reduction:

For reference, the upper-centre green blob is about 3 arcmin across, of order the array footprint. The smaller-scale dark features at the right are simply an artifact of the real Orion 850um emission that extends into the patch of the map where we inserted the Vela data!

Save time and memory: down-sampling your data

2010-10-26T11:59:00.008-10:00

I've finally gotten around to adding the ability to "down-sample" the data.

SCUBA-2 data are normally sampled at about 200 Hz (actually closer to 180 Hz). This specification was based on the desire to scan up to 600 arcsec/sec while still fully-sampling the 450um beam. With these rates, a single sample covers an angular scale on the sky of 600/200 = 3 arcsec, which corresponds to approximately the 450 FWHM/3.

For most of the observations taken as part of SRO, however, typical scan rates were in the range 100--200 arcsec/sec, meaning that SCUBA-2 samples faster than necessary. We can then, in theory, re-sample the data into longer samples without losing any useful information. SMURF can now do this by specifying the smallest scale you wish to sample in the configuration file:

450.downsampscale = 2

850.downsampscale = 4

These values are measured in arcsec, and by setting to these values you will match the default pixel resolution used by SMURF in each band. It will then determine internally, based on the previous sample rate and mean slew velocity, how much to downsample the data.

The following example is the default reduction of an Orion map at 450um (20100219_00039) using the dimmconfig_bright_extended.lis configuration file:

And then again, except setting 450.downsampscale = 2:

As you can see there is very little difference between the two images.

As makemap is running you will notice:

smf_grp_related: will down-sample from 174.6 Hz to 60.0 Hz

Indicating that the data are being compressed by a factor of about 3. In terms of execution time on my 8-core machine it drops from 18 s to 11 s. The reason it is not a full factor of 3 is because there is a large initial start time as the original data (before down-sampling) are loaded in. Finally, the estimated memory usage drops from 457 MiB to 251 MiB (again, not a full factor of 3 due to other static memory that needs to be loaded).

This test uses a very small data set. The improvements are more pronounced for long scans in which the total memory usage is dominated by the time-series data, and not other statically allocated buffers (such as the map etc.). There will also be larger improvements at 850 um than at 450 um (because it is lower resolution), and for slower scan speeds.

This blog

2010-10-11T13:19:00.000-10:00

You may have noticed that the posting rate has picked up quite a bit in the last few months, as we've been trying to keep everybody up to date with externally visible developments to the JCMT Science Archive.

If you are a reader for this blog, please take a moment to leave a comment and let us know if you would like more/fewer posts, whether the level of detail is too little/too much, the kind of topics you'd like covered, or anything else that you would like to share.

S2SRO products - new re-reduction

2010-10-11T13:00:00.000-10:00

The S2SRO products have been regenerated once again. In this batch you will find

Much improved map for extended emission targets (which now are specially handled)
Improved step-finder
S/N masking for bright sources to prevent negative bowling
Point source observations now reduced with the science recipe rather than the summit "quickie" recipe.

As usual you will find your SCUBA-2 products in the proprietary processed data search - just put in your project code and choose multi-night to get your total product.

(Note that this batch does not include the map-based despiking - that will be in the next run).

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

2010-10-04T12:14:00.014-10:00

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

2010-09-30T14:29:00.007-10:00

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

2010-09-21T16:35:00.001-10:00

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

2010-09-13T09:04:00.013-10:00

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

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

2010-09-09T13:20:00.003-10:00

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').

Processing date now visible in JSA web UI

2010-09-08T13:09:00.000-10:00

By popular request, you now have access to the processing date of your product when you do a JSA search - it's the field at the bottom of the "Program Constraints" section called, unsurprisingly, "Processing Date". Like all the other date fields, you can enter a range to search, or if you just leave it blank but tick the box, you will get the extra column of information in your results page.

Many thanks to Dustin and Sharon for rolling this out so quickly.

SCUBA-2 products update - fresh batch

2010-09-07T13:29:00.000-10:00

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.

JSA FAQ: examining data reduction parameters

2010-09-01T15:10:00.002-10:00

If you know what hislist and provshow do, you can skip this entry.

Say you have downloaded your most excellent project product from the JSA and you are interested to see what went in it. Here's what you do (assuming, of course, that you have already sourced your Starlink startup files).

Convert the FITS file back to NDF (we'll call the result example for brevity):

; convert

; fits2ndf jcmts20100228_00040_850_reduced001_pro_000.fits example

So we now have our example.sdf

To see the history of the file, run the KAPPA hislist (history list) command:

; kappa

; hislist example

In the voluminous output, you will see the full parameters of any command ran on this file, for example, the makemap parameters that were used to generate this SCUBA-2 product:

In order to see what observations went into this product, use the KAPPA provshow (provenance show) command:

; provshow example

which will show you the entire provenance tree for your file.

Thanks to the NDF history and provenance mechanism, every operation on the data files is automatically recorded, so you have a complete record of what was done to the data.

About the JSA Products

2010-08-11T11:37:00.000-10:00

Judging from my inbox, there is quite a lot of enthusiasm about the S2SRO products available from the JSA. We sure love this project, and it's great to see it hitting its stride and being useful to our external users. Still, I just want to throw in a couple of cautions so that people don't get caught out.

We do not have official data releases, as most people understand them. By this I mean that the data is not vetted by anybody at the JAC - we simply don't have the effort for that. Since right now we are still developing, we do trawl for obvious problems and bugs, but there's no scientific oversight of what the processing churns out, and the data is immediately exposed for download. So of course you may take the products, and a lot of them do seem to be publication quality, but you should still work through the reduction cookbook and make sure you understand what was done to your data. The result you download shouldn't be different from what you would get if you run our latest software at home with the recommended parameters. The idea is to get you the best data we can give you now, rather than a perfect version of your data later. The processed products can also help you prioritise which datasets to spend most time on.

As our data pipeline matures, or as we fix new bugs, we do re-reduce the data - and every time we do this the new product replaces the old one. So if you intend to post-process and/or publish using downloaded products, make sure you retain the version of the product that you used in case you need to reproduce your work later.

We do have a plan to allow users to upload their own versions of products into the JSA, but this is still a way off.

SMURF Update (August 5th 2010)

2010-08-05T10:00:00.000-10:00

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.

New version of the SCUBA-2 cookbook

2010-08-04T22:54:00.001-10:00

Version 1.1 of the SMURF SCUBA-2 Data Reduction Cookbook (or SC19, as it is fondly known), is out. We recommend that everybody who works with SCUBA-2 data read through it.

In order to be able to follow the cookbook, you will need to update your local starlink release to the latest version.

You can always get to SC19 from the sidebar of this blog.

JSA FAQ: Product grouping types

2010-08-04T22:32:00.000-10:00

When you search the CADC archive for processed data (either public or proprietary), you will have the option selecting any of four different types of product. These are:

Simple: This is the most processed state of a single observation
Night: This is the product resulting from all observations taken in one night on the same field
Multi-night: This is the product resulting from all observations taken in one field, even if they were taken on different nights. We sometimes call this the "project" co-add, because once observing for that project is finished, it represents all the (groupable) data taken for the project in that field
Public: This is a product consisting of all public observations of a particular field, even if they were taken for multiple projects. At this time, we have not generated any products of this type.

We sometimes refer to the simple product as the "obs" product (after the filename suffix that you will get when you download it). We also refer to products 2-4 as "aggregate" products, because they consist of more than one observation. You should always get an obs product for each observation in your project; aggregate products are made excluding any observations marked as bad.

Here's where you can find this option in the JCMT Science Archive:

Bad, bad, BAD observation!

2010-08-04T21:24:00.002-10:00

This post is an explanation of where we are at the moment with handling quality in the JSA, and what the medium and long-term plan is. Since the topic is bad observations, this should be of particular interest to S2SRO users with very early SCUBA-2 data, since happily, ACSIS observing doesn't result in many of those!

What happens right now

At this time, when we process observations in batches to create night and multi-night (project) co-adds, the system does not use any observation that has been marked as bad in the OMP obslog (observations are good by default). In the case where only one sub-system is bad (for example in the SCUBA-2 case of the 450 being good and the 850 being marked bad) both observations are excluded. Questionable and rejected observations are included in the co-add.

People who are able to set an observation's status to BAD include JAC staff, the active observer at the telescope and the data owners (PIs and co-Is associated with the project). We track who changes the status of an observation, and people are encouraged to leave an explanatory comment as to why they did so.

In the near term

The problem is right now that we do not have enough eyeballs to look closely at the S2SRO data and assess whether every observation is good. In the near term, we know people are inspecting their data and really would like the data owners to take the time to mark an observation as bad in their OMP project pages. Then you can either ask for your products to be re-generated right away, or wait for the next re-processing run (SCUBA-2 data is being reprocessed frequently as the data reduction pipeline improves).

I also think that rejected observations (those that are technically good, but failed to meet a survey's particular QA criteria) should be excluded from the aggregate products - this is an issue we will take up with the surveys.

This is how you can mark your observation as bad when you are not at the JAC:

[In this example, let's say we have already identified that we don't want observation 68 taken on UT 2010-03-06 included in our aggregate products. For the impatient: OMP home page -> Access a project -> Pick your UT date -> Click on Comment]

The real plan

Obviously the current state of affairs is sub-optimal. The two major improvements that are planned in observation quality are:

Allow individual sub-systems to be marked as bad, rather than throwing the baby out with the bathwater. The problem with this is that obslog (which long predates archive processing) only understands the observation level, and so there are significant OMP infrastructure changes that need to be made to allow this.
Develop an interface between ORAC-DR and obslog that will allow the pipeline itself to mark observations as bad (not surprisingly, the most sure-fire to find bad observations is to read what ORAC-DR is telling you).

Both of these are on the cards, but the reality is that they are a lower priority that the main SCUBA-2 work, so it's not possible to promise a timeline for their delivery.

If you zoned out reading the above:

If you find bad observations in your data, take the time to mark them as bad in your project web pages. Watch the video to find out how.

JSA FAQ: Finding your products

2010-08-03T11:34:00.000-10:00

A few folk have had trouble figuring out how to get their proprietary products (processed data). You can do this by making the appropriate JSA query.

Here's how to do this starting from the JCMT home page

And here is how to do this starting from the CADC home page:

Don't forget that you have to use your CADC credentials for this operation, and they have to be associated with your JAC/OMP userid (so that the OMP can tell the CADC system that it is okay for you to access that data). If the above instructions do not work for you, it is likely that this linking of the two accounts has not been done; contact a JCMT staff member to do the deed. You get a CADC userid by applying to their site and picking a username of your choice; your OMP userid was issued to you when you successfully applied for time and typically consists of your last name followed by your first initial.

SMURF update (June 4th 2010)

2010-06-04T16:03:00.000-10:00

We've been making incremental changes to SMURF over the last few weeks so here are a few highlights that have not been covered in other blog entries.

We've updated the extinction correction parameters that scale from CSO to the relevant filter to use the number just calculated by Jessica Dempsey. You can override these values by setting ext.taurelation.filtname in your map-maker config files to the two coefficients "(a,b)" that you want to use (where filtname is the name of the filter). The defaults are listed in $SMURF_DIR/smurf_extinction.def. We have also added a slight calibration tweak to WVM-derived values to correct them to the CSO scale.
We have changed the calling interface for the sc2clean command so that it now takes a map-maker config file as input rather than individual command line parameters. This should make it easier to compare what the map-maker is doing to what sc2clean is doing.
Two new (experimental) models have been added to the iterative map-maker. They are slow and they are experimental:

SMO will do a boxcar smooth to each time series to calculate the low frequency variation. This might be more reliable than using the FLT model if the only aim is to remove low frequencies.
PLN will fit and remove a plane from each time slice.

There have been some enhancements to the step finding code so that it not only finds the steps more accurately but can also find correlated steps that occur at the same place for all bolometers.
The map-maker has a new ast.zero_notlast parameter that can be set to true when used in conjunction with ast.zero_lowhits to disable the final iteration from being forced to zero.

All these changes are available in the 64-bit linux stardev rsync location and also in the most recent Mac Snow Leopard tar ball.

Extinction correction factors for SCUBA-2

2010-06-04T13:26:00.000-10:00

Analysis of the SCUBA-2 skydips and heater-tracking data from the S2SRO data has allowed calculation of the opacity factors for the SCUBA-2 450μm and 850μm filters to be determined.

Some background: the Archibald et al (2002) paper describes how the CSO(225GHz) tau to SCUBA opacity terms were determined for the different SCUBA filters. It was assumed for commissioning and S2SRO that the new SCUBA-2 filters were sufficiently similar to the wide-band SCUBA filters that these terms could be used for extinction correction. For reference the SCUBA corrections were:

Tau(450μm) = 26.2 * (Tau(225GHz) - 0.014)

and

Tau(850μm) = 4.02 * (Tau(225GHz) - 0.001)

The JCMT Water vapour radiometer (WVM) now is calibrated to provide a higher-frequency opacity value which has been scaled to the CSO(225GHz) tau. The WVM (not the CSO 225GHz tipper) data was used for this analysis.

The new filter opacities as determined by the skydip data are as follows:

Tau(450μm) = 19.04 * (Tau(225GHz) - 0.018)

and

Tau(850μm) = 5.36 * (Tau(225GHz) - 0.006)

A follow-up post to this will show analysis of the difference applying the new corrections can make to data combined from multiple observations taken in differing extinction conditions.

It is worth noting that if an individual science map and corresponding calibrator observation is already reduced with the old factors (and your source and calibrator are at about the same airmass and if the tau did not change appreciably), any errors in extinction correction should be cancelled out in the calibration.

Applying FCFs to calibrate your data

2010-06-04T13:03:00.000-10:00

Calculating SCUBA-2 Flux Conversion Factors (FCF's)

Currently SCUBA-2 reduction software: the pipeline and the PICARD recipes produce three separate FCF values. Details of the PICARD recipes can be found on Andy's PICARD page.

For calibration from point sources the FCFs and NEFD's have been calculated as follows:

The PICARD recipe SCUBA2_FCFNEFD takes the reduced map, crops it and runs background removal (and surface fitting parameters are changable in the parameter file).
It then runs the Kappa beamfit program on the specified point source. Calibrators such as CRL618, HLTAU, Uranus and Mars are already hard-coded into the recipe. If it is not, then you can add a line to your parameter file with the known flux: FLUX_450 = 0.050 or FLUX_850=0.005 for example. Beamfit will calculate the peak flux, the integrated flux over a requested aperture (30 arcsec radius default), and the FWHM etc.
It then uses the above to calculate three FCF terms described below.

FCF (arcsec)

FCF(arcsec) = Total known flux (Jy) / [Measured integrated flux (pW) * (pixsize²)]

which will produce an FCF in Jy/arcsec²/pW.

This FCF(arcsec) is the number to multiply your map by when you wish to use the calibrated map to do aperture photometry.

FCF(beam)

FCF(beam) = Peak flux (Jy/beam) / [Measured peak flux (pW)]

producing an FCF in units of Jy/beam/pW.

The Measured peak flux here is derived from the Gaussian fit applied by beamfit. The peak value is susceptible to pointing and focus errors, and we have found this number to be somewhat unreliable, particularly at 450μm. This FCF(beam) is the number to multiply your map by when you wish to measure absolute peak fluxes of discrete sources.

To overcome the problems encountered as a result of the peak errors, a third FCF method has been derived, where the FCF(arcsec) is taken and modeled with a Gaussian beam with a FWHM equivalent to that of the JCMT beam at each wavelength. The resulting FCF calculates a 'equivalent peak' FCF from the integrated value assuming that the point source is a perfect Gaussian.

FCF (beamequiv)

FCF(beamequiv) = Total flux (Jy) x 1.133 x FWHM² / [Measured integrated flux (pW) * pixsize²]

or more conveniently:

FCF(beamequiv) = FCF(arcsec) x 1.133 x FWHM²

where FWHM is 7.5'' and 14.0'' for the 450μm and 850μm respectively. This produces an FCF in units of Jy/beam/pW.

This FCF(beamequiv) and FCF(beam) should agree with each other, however this is often not the case when the source is distorted for the reasons mentioned above. FCF(beamequiv) has been found to provide more consistent results and it is advised to use this value when possible, in the same way as FCF(beam).

Methodology for calibrating your data:

So you have a reduced map for a given date. Each night of S2SRO should have at least one if not more calibrator observations that were taken during the night. A website with this list is currently in the works and I'll add it to the blog when it is completed. For now it is relatively easy to search the night's observations for these. Primary calibrators were Uranus and Mars, and the secondary calibrators are listed on the SCUBA secondary calibrators page. In addition to these, Arp 220, V883 Ori, Alpha Ori and TW Hydrae were tested as calibrators. Their flux properties were investigated with SCUBA (see SCUBA potential calibrators) .

Run your selected calibration observation through the mapmaker using the same dimmconfig as your science data used.
Use PICARD's recipe SCUBA2_FCFNEFD on your reduced calibration observation. This will produce information to the screen and a logfile log.fcfnefd with the three FCFs as mentioned above, and an NEFD for the observation. PICARD by default uses fixed FCF's to calculate the NEFD. (450um: 250 and 850um: 750). If you wish to get an NEFD using the FCF calculated for the calibrator add USEFCF=1 to your parameter file.
Take your selected FCF and multiply your map by it using KAPPA cmult.

Things become slightly more complicated if you wish to use PICARD's matched filter recipe to enhance faint point sources. Again see Andy's PICARD page (link above) for details on the matched filter recipe. If you are normalising the matched filter peak, you will need to run this filter over your calibrator image with the same parameters you used for your science map.

Note: You cannot, at this point, use the FCF(beamequiv) number to calibrate your match-filtered data. This number will now be (usually) disproportionate and just wrong. The FCF(beam) value however, should be preserved by this method.

The peak is truly preserved by this method so two numbers, the FCF(beamequiv) pre-match-filter and the FCF(beam) post-match filter should be close to the same and either of these values can be used to calibrate your match-filtered science map.

It is also worth noting (though perhaps obvious) that after running the match-filter script in peak-normalisation mode, only the peak flux values (and not the integrated sum over an aperture) will be correct. The reverse is true if using sum-normalisation.

Dark Squids and Magnetic Field Pickup

2010-05-28T10:02:00.000-10:00

Over the last couple of weeks I have been experimenting with some data that seem to exhibit the effects of magnetic field pickup in the detectors. The SCUBA-2 readout system amplifies small changing magnetic fields caused by changing currents in the detectors (that are produced by changing optical and thermal loads on the bolometers) using arrays of SQUIDs. Although there is substantial shielding, some external changing magnetic fields may still leak into the raw data - such as motion of the telescope through the Earth's magnetic field, or fields produced by motors in the dome etc. The readouts are multiplexed, and the signals from columns of detectors pass through a common amplifier chain. For this reason, each column has a bare SQUID with no bolometer, a so-called "dark SQUID" (DKS henceforth), to track this non-thermal/optical noise term. The time-series for these data are stored in raw SCUBA-2 data files in the ".more.scuba2.dksquid" extension. They can be plotted, for example, using KAPPA linplot (this is for the DKS in column 1 of the following raw 30s data file):

$linplot 's4a20100223_00061_0002.more.scuba2.dksquid(1,)'

Just like regular bolometer data, DKS can have spikes (like in this example). Generally speaking the DKS haven't shown a lot of structure so we weren't concentrating on using them to clean the bolometer data.

However, on some of the wider-area data sets the DKS signals can be substantially larger, and seem to introduce a significant noise term. One additional feature of the magnetic field pickup signal that makes it easy to distinguish from other noise sources is that its sign is essentially random, and independent of the detector response to optical/thermal loads, i.e. the signal could be rising in one bolometer, and falling in the next. The following movie shows a data cube from the s8d array after the common-mode signal has been removed. The residual signal oscillates in time, with two clear detector populations that are out of phase with eachother:

These data were taken as part of a 0.5 deg diameter scan, and it is clear comparing the dark squid signals with the change in telescope azimuth that they are highly correlated. Combined with the arbitrary sign of the residual noise, the culprit is almost certainly magnetic field pickup.

A dark squid signal (red dots), and the signal smoothed with 200 sample boxcar (white line) over 100s:

Azimuth offset from the map centre over the same period of time:

The next plot shows the dark squids for each working column over ~25 min. with arbitrary vertical offsets for clarity. Note that while similar, each column has a slightly different response:

The iterative map-maker has a partially-implemented model, called "DKS", for iteratively fitting and removing the DKS signals from working detectors along each column. The relevant settings for the configuration file are:

modelorder = (dks,com,gai,ext,ast,noi)
dks.boxcar = 100

Here we have added DKS cleaning as the first step in the iterative model fitting to the data (note also that we have dropped FLT for now, more on that below), and the second option smooths the DKS data with a 100 sample boxcar prior to fitting and removing it from each bolometer.

The next two images show a map solution using only COM to remove low-frequency noise, followed by a solution that uses both DKS and COM:

While not perfect, it is still clear that the DKS removal has removed a lot of the noise that caused streaking in the first image. Remember, there is no additional high-pass filtering as FLT has been turned off in this example.

If "dks" is specified in "exportndf", the fitted DKS model will be written to a file with a suffix "_dks.sdf". However, this file is not particularly easy to view and make any sense of because the dark squids and the fit coefficients are all bundled together in one large array. For this reason, we've added a new SMURF task called "sc2expandmodel" that will evaluate the model and produce a full 3d data cube with the same dimensions/units as the residual data, e.g.

$sc2expandmodel s4a20100228_00016_0002_con_dks expanded_dks

will produce an easily-viewable file called "expanded_dks". All detectors along a column will have the same shape, but the scale/offset will be different for each detector. Examining this file and comparing to the "_res" output from makemap is useful for getting an idea of how strong this magnetic field pickup is for a given data set.

All of this code is still experimental. It currently has several problems:

boundaries with padding/apodization will introduce steps, so padstart/padend/apod should probably be set to zero, meaning that it is not optimal for combining with FLT at this stage.
the DKS are not pre-processed in the same way as the bolometer data to remove spikes / steps etc. This capability will be added shortly.
Some DKS are clearly dead (e.g. DKS 1 in the multi-DKS plot above). Currently we don't have an automatic system for flagging/removing them. Note also that if they are bad, it will be impossible to DKS-clean working detectors along that column (so we need to decide if we will simply have to ignore them when making a map?)

Once we have worked through some of these issues, we will update the default dimmconfig files in SMURF accordingly.

SCUBA-2 filter information

2010-05-18T09:46:00.000-10:00

At our last SCUBA-2/SRO telecon, people asked to see the measured profiles for the SCUBA-2 filters.

Dan has put up some plots here.

PICARD web page

2010-05-17T08:55:00.000-10:00

I've created a new static location for all things PICARD:
http://www.oracdr.org/oracdr/PICARD

The page includes an introduction to PICARD, a list of all available recipes (with more detailed documentation) and a few hints, tips and potential gotchas. I'll keep this page up-to-date, adding new blog entries and/or sending messages to the scuba2dr mailing list as necessary.

Post-processing SCUBA-2 data with PICARD

2010-05-11T13:25:00.001-10:00

Processing raw SCUBA-2 data is done with SMURF or the pipeline (ORAC-DR). What happens after that depends on the user and their level of familiarity with particular software packages. Fortunately the SCUBA-2 software team is here to help out and come up with a series of standardized tools for performing a number of basic post-processing tasks.

Introduction to PICARD

Our tool of choice is PICARD which makes use of the existing ORAC-DR infrastructure as well as our existing knowledge of writing primitives and recipes for the SCUBA-2 pipeline. PICARD is run from the command line as follows (I'll use % as the prompt):

% picard [options] [RECIPE_NAME] [list of files to process]

For example,

% picard -log sf -recpars mypar.lis CROP_JCMT_IMAGES myfiles*.sdf

The most commonly used options are -log and -recpars (the full list of available options can be seen by running picard -h). Both of these options take additional arguments.

The "-log" option controls where the messages from PICARD are printed: "-log sf" will write messages to the terminal window and to a file called .picard_PID.log (where PID is the process ID for picard) in the output directory. To avoid creating the .picard_PID.log files, just specify "-log s".

The "-recpars" option allows the user to pass in a text file containing parameters which can be used in the given recipe. The permitted parameters are listed with the various recipes below. The format of this text file is a list of `parameter = value' entries, with the recipe name given in square brackets:

[RECIPE_NAME]
PARAM1 = VALUE1
PARAM2 = VALUE2

PICARD writes its output files to the current directory (unless the environment variable ORAC_DATA_OUT is defined in which case that location will be used).

There are currently four recipes which may be of interest:

MOSAIC_JCMT_IMAGES
CROP_JCMT_IMAGES
REMOVE_BACKGROUND
SCUBA2_MATCHED_FILTER

These recipes and their parameters are described in more detail below. More recipes will be added as the need arises and as we gain more experience in analyzing SCUBA-2 data. Interested users should update their Starlink installations to get access to these recipes.

MOSAIC_JCMT_IMAGES

Coadd the given files into a single map, taking into account the EXP_TIME and WEIGHTS NDF components. The images are combined using variance weighting and the output variance is derived from the input variances. Currently the recipe uses the KAPPA wcsmosaic task for coadding the images.

The same pixel-spreading method (and any associated parameters) is used for the data and the EXP_TIME and WEIGHTS.

Creates a single output file based on the name of the last file in the list, and with a suffix "_mos" (e.g. mylastfile_mos.sdf).

Available recipe parameters:
[MOSAIC_JCMT_IMAGES]
WCSMOSAIC_METHOD = wcsmosaic pixel-spreading method: see wcsmosaic documentation for available options (default is "nearest")
WCSMOSAIC_PARAMS = additional parameters which may be required for the chosen method

CROP_JCMT_IMAGES

Crop images to the map size in the data header (as specified in the Observing Tool), though this size can be overridden using the recipe parameters below.

Creates an output file for each input file, with the suffix "_crop".

Available recipe parameters:
[CROP_JCMT_IMAGES]
MAP_WIDTH = map width in arcsec
MAP_HEIGHT = map height in arcsec

REMOVE_BACKGROUND

Fit and remove large-scale background variations from images using either KAPPA fitsurface or CUPID findback. See the Starlink documentation on both tasks for more information on the parameters shown below. The option exists to mask out a circular region centred on the source before removing the background.

Be aware that the subtraction of the background will add noise proportional to the RMS deviation between the image and the background fit.

Creates an output file for each input file with the suffix "_back".

Available recipe parameters:
[REMOVE_BACKGROUND]
MASK_SOURCE = flag to mask out a circular region on the source before fitting a background (1 = mask out source; 0 = do not mask out source - the default)
APERTURE_RADIUS = radius of aperture (in arcsec) for masking out source (otherwise 30 arcsec)
BACKGROUND_FITMETHOD = the method for fitting the background, either fitsurface (default) or findback
FITSURFACE_FITTYPE = fittype parameter for fitsurface: polynomial (default) or spline
FITSURFACE_FITPAR = parameters for fit, up to 2 numbers corresponding to NXPAR/NYPAR for fitsurface or KNOTS for spline fit
FINDBACK_BOX = size of box in pixels used by findback for smoothing the image

SCUBA2_MATCHED_FILTER

Apply a matched filter to the data to improve point-source detectability. The images and PSFs are smoothed with a broad Gaussian (default is 30 arcsec but can be varied using the recipe parameter below) and subtracted from the originals. The images are convolved with the modified PSFs. The PSF created by the recipe is a Gaussian with FWHM equal to the JCMT beamsize at the appropriate wavelength (i.e. 7.5 or 14 arcsec).

Creates an output file for each input file with the suffix "_mf", and a PSF file "_psf" if the PSF was not specified as a recipe parameter.

Available recipe parameters:
[SCUBA2_MATCHED_FILTER]
PSF_MATCHFILTER = name of a PSF file (NDF format, will be used for all images)
PSF_NORM = switch to determine whether a PSF is normalized to a peak of unity ("peak" - the default) or a sum of unity ("sum")
SMOOTH_FWHM = FWHM in arcsec of Gaussian to smooth image (and PSF)

Typos in config files

2010-05-07T15:37:00.000-10:00

This week I've made a couple of minor changes to the way the iterative map-maker handles configuration files. On too many occasions I've made a change to a config parameter and seen no effect because of a typo or else I've left a parameter unspecified and had no idea what that really meant.

We've now fixed both these problems:

Defaults are now read from a file called $SMURF_DIR/smurf_makemap.def
When user-supplied parameters are read in they are merged with the defaults and if there is any keyword used that is not present in the default file you will get an error immediately.
When the history is written out (use KAPPA HISLIST to see it) it contains the values that were actually used rather than just the values that were supplied by you. These are much more useful when trying to work out exactly what is different between two reductions.

I did have to change the way config files work in one small way. Previously if you wanted your config file to include overrides based on the wavelength you would have to say "flt_850.something" or "flt_450.something". To allow the above changes to work we've had to reorganize that slightly so you now say "450.flt.something" and can now override all the normal parameters in this way for a particular wavelength.

How Should I Mosaic My SCUBA-2 Data (Redux)?

2010-04-29T11:56:00.001-10:00

Mosaicking your SCUBA-2 with wcsmosaic or makemos is fine if all you're interested in is the data. But what if you want to know the exposure time per pixel in your map to determine if the noise is reasonable? That's easy for each map written by makemap: the exposure time image is stored in the .MORE.SMURF.EXP_TIME NDF extension (which you can view in Gaia). But wcsmosaic doesn't know about this extension, so the EXP_TIME data in a mosaic is only that of the first file.

There are two ways to get properly mosaicked data:

Use the ORAC-DR pipeline and process from raw;
Use the PICARD recipe MOSAIC_JCMT_IMAGES on processed data.

I'll write a post about using the pipeline in the near future but for now I'll highlight the second method.

PICARD comes for free with ORAC-DR (which comes for free with Starlink). It's basically a tool for processing and analyzing reduced data which takes advantage of the same infrastructure used by the pipeline.

At the command line type:

picard -log s MOSAIC_JCMT_IMAGES myfiles*.sdf

It'll tell you what it's doing and exit. This will create a file called "mylastfile_mos.sdf" which has the correct EXP_TIME extension (where "mylastfile" is the name of the last file in the list given to MOSAIC_JCMT_IMAGES). Try it and see.

Note that all the files must be of the same source.

At the moment, it only supports wcsmosaic as the mosaicking task but will support makemos in the future. Keep rsync'ing and one day it'll just be there...

[This post was updated on 20100510 to reflect the change in the recipe name and the file-naming behaviour]

How Should I Mosaic My SCUBA-2 Data?

2010-04-29T10:53:00.000-10:00

I've been asked this a number of times so I thought I'd write something down. Currently the SMURF iterative map-maker calculates models for each chunk of data independently (a chunk is a continuous sequence so many times that will be a single observation but it might be smaller than an observation if a flatfield is inserted mid-way through). It then combines them using a weighted coadd to give you the final map. This means that there is nothing to be gained in throwing all your observations at MAKEMAP apart from it taking much longer to make a map before you can see the result.

In terms of flexibility it is better to make a map from each observation separately. You can then look at each map to see how it looks before combining it. You may also want to remove any low frequency structure at this point as well. To combine these maps into a single map you have two choices using Starlink software:

Use KAPPA WCSMOSAIC with an interpolation or rebinning scheme of your choice. Make sure that you set the VARIANCE parameter to true to enable variance weighting of your mosaic.
If you have already ensured that the maps are made on the same grid (you can use the REF argument in MAKEMAP to ensure this) then you can try CCDPACK MAKEMOS. This can be used to do median or clipped mean stacking and does no interpolation or smoothing. If your images are not aligned on the same grid it is probably better to remake them so that they are, but if that is difficult you can use the KAPPA WCSALIGN command to match up the grids first. Use the GENVAR and USEVAR parameters in MAKEMOS.

In the future it may be that the map-maker will be able to do a better job handling all the data at once itself instead of external mosaicking but at the moment this is not the case.

If you want to keep a track of exposure times you need to track that separately. We will cover that in a follow up post.

Large scale structure divergence

2010-04-23T12:55:00.000-10:00

First, my apologies for the lousy image layout! Next post will be better...

As many of you have noticed, large-scale noise structures in the maps appear to grow with the number of iterations (clearly not the desired behaviour!). This blog entry describes the problem in more detail with an example, and discusses a partial fix.

The map on the left was produced from about 6 min. of scans across a blank-field at 450um using the dimmonfig_faint.lis configuration file, but turning up the number of iterations to 100. As you can see there is a strong gradient across the map -- the more iterations I added, the more the gradient grew.

In the next plot I show the common-mode signal from this map solution (the average signal in all the detectors as a function of time). Mostly we see the fridge temperature oscillations (about every 30s), plus a more slowly varying component (maybe sky?)

Next, I set ast.zero_lowhits = 0.5 in the config file. What this option does is identify all of the pixels in the output map with less than 50% of the mean number of hits, and then sets those pixels to zero after each iteration (i.e. a zero boundary condition). In the next two plots I show the map, and the difference between the common-mode estimated from the two solutions:

What we can can see is that the gradient has been drastically reduced. So what's happening? The common-mode difference plot shows that the map without edge constraints has a large additional periodic signal (about 10% of the main common mode signal) with a frequency that is a factor of a few larger than the 30s fridge oscillations. In fact, this signal is nicely correlated with the higher-frequency scan pattern.

Basically this shows that, left to its own devices, the solution is perfectly happy to pump large-scale structure into the map. It then still manages a flat residual because it simply compensates for this gradient by putting negative signal of the correct amplitude into the common-mode. In other words, the largest scales and the common-mode are degenerate model parameters.

However, I have noticed that even with this edge constraint turned on, there is other strange structure in the map, like this convex curvature. This fact, combined with the substantially worse noise at 850um, shows us that the model is insufficient for accurately describing the data. The edge constraint helps, but that curvature indicates some kind of tension (like the model trying to fit some sort of gradient across the array and failing). I suspect that, lacking any other conditions on the model, the large-scale noise is just growing without bound.

We are presently looking at ways to improve the model for the data and/or introduce other reasonable constraints to help convergence to something sensible! For now, it is worth turning ast.zero_lowhits on if you don't believe you have any strong emission features around the edge of the map (tricky for maps of structure in the galactic plane).

Detection and Correction of DC steps

2010-04-06T02:50:00.000-10:00

The new makemap DC step detection and correction routine (smf_fix_steps) processes each bolometer independently. The processing is divided into two parts; first, potential steps are detected, and then second, each potential step is checked to see if it satisfies various conditions. Once a final set of steps has been determined, the original data stream is modified to correct for these steps. The location of the steps themselves are flagged in the associated quality array.

1) Detection of Candidate steps:

Candidate steps are detected by comparing the difference between the median value in two boxes on either side of each time slice. In the figure below, the time slice with index t0 is being tested by comparing the median of the values in the two red boxes. The number of times slices in each of these median boxes is given by config parameter DCMEDIANWIDTH (default 40 - a small value is used to reduce the smoothing effect which can make steps harder to find). A gap is allowed to exist between these two boxes (controlled by the private parameter dcmediangap, currently set to 3 samples). Each time slice is processed in this way to produce a curve like the one below, which is indicative of the gradient of the bolometer time series. The effect of the median box is to blur sudden jumps so that high gradients are seen over a width comparable to DCMEDIANWIDTH.This curve has a high positive value around the location of the steep upward step in the above curve. Use of a median minimises the effects of spikes, etc.

Unusually large (positive or negative) values are identified in this curve. The threshold is set as a multiple of the RMS value in the curve, with the multiple being given by the DCTHRESH config parameter.
This thresholding results in a miscellaneous collection of samples being selected, which will often, but not always, be associated with a DC step. A genuine step may not occur instantaneously - that is, the rise in value may take several time slices to complete. In addition, the rise may not be monotonic. This means that a single step may be broken up into pieces because of the thresholding. The next step in the process is to merge these pieces together into contiguous blocks. Gaps between selected samples that are shorter that dcminstepgap (a private config parameter currently set to 50), are filled in (i.e. the low gradient samples in the gap are added to the selection of high gradient samples).
After filling in the gaps, any blocks of contiguous selected samples that are shorter than dcminstepwidth (=0.8*DCMEDIANWIDTH) or longer than dcmaxstepwidth (=1.8*DCMEDIANWIDTH) are rejected. The lower width limit is imposed because an instantaneous step will be blurred by the median filtering described above and will have some larger width.
A further restriction is that for a block to be accepted, it must contain predominantly either positive gradient values or negative gradient values. Blocks that contain roughly similar numbers of both are rejected. A block is rejected if the difference between the number of positive and negative values within the block is less than dcminsignratio (=0.8) times the total number of values in the block.
Yet a further restriction is that for a block to be accepted, the median value at the start and end of the block - after filling - must still differ by more than the above threshold value.

This process leaves us with a collection of contiguous blocks of selected samples. Each block is a candidate step with a starting time and an ending time (i.e. each block potentially contains the rise or fall of a step).

2) Testing of Candidate Steps:

Each block of contiguous selected samples found above is now tested to see if it is a real step. This is done as follows:

A least squares linear fit is performed to the sample values just before the step. This fit is performed on a box containing the number of samples specified by the DCFITBOX config parameter (default 400), ending at the lower edge of the step. This fit does two sigma-clipping iterations to minimise the effects of spikes etc on the line. The number of sigma at which to clip is specified by dcthresh2 (private parameter, currently set to 7). This clipping is kept light to avoid loosing the effect of other more gradual undulations in the line. The RMS deviation of the data (minus spikes) about this line is recorded.
A line is fitted to the data just after the step in the same way. The following figure illustrates these two fits - the solid red lines are the fits themselves, and the dotted red lines mark the start and end of the step, as found by the median filtering process described above.
These two lines are used to derive two estimates of the data value at the centre of the step (i.e. mid way between the dashed lines in the above figure). If the difference between these two estimates is larger than a threshold value, the step is accepted as genuine. The threshold value is DCTHRESH times the larger of the two RMS residuals found whilst fitting the two red lines.
If the step is accepted, the difference found in step 3) is subtracted from all subsequent samples following the end of step (i.e. all samples to the right of the right hand red dotted line in the above figure). Earlier sample values are left unchanged. The samples that fall between the two dotted red lines are flagged with the SMF__Q_JUMP quality.
Once all candidate steps have been tested, a final correction is made to all samples in order to ensure that the mean value in the entire data stream is unchanged.

There are a lot of parameters in this process, only some of which are currently public (the ones named in upper case above). Finding appropriate values for these could take some time. Also, the algorithm seems to be fairly fragile, and can miss some visually obvious steps, whilst accepting as steps things which do not really look like steps visually (e.g. steep slopes). Further tweaking will probably be necessary.

A problem occurs where the DC offset gradually rises (or falls) and then suddenly flips back to its original value. The above algorithm often fails to spot the gradual rise, but successfully spots the rapid fall. Correcting for the fall but not the rise can lead to a systematic offset for the remaining data. Hopefully, the blocking recently introduced in the common mode estimation will allow such offsets to be detected and corrected. We shall see.

SMURF and quality flagging

2010-03-25T12:34:00.000-10:00

As promised in one of the recent operations telecons, a blog post about quality flagging, and how all the new quality flagging statistics can help you understand what's going on.

First some background. The raw data is a cube, containing 2 spatial dimensions, and 1 time dimension. There is a single quality array, with these same dimensions, that is used by the map-maker to keep track of various problem areas that are identified both before and during the iterations. Each element of this array is an unsigned byte, and each bit is a flag. You can list the meaning of each bit using Kappa showqual on an exported "_res" model component (with "qua" also specified in exportndf in the dimmconfig file):


BADDA           (bit 1) - "Set iff a sample is flagged by the DA"
BADBOL          (bit 2) - "Set iff all data from bolo to be ignored"
SPIKE           (bit 3) - "Set iff a spike is detected"
DCJUMP          (bit 4) - "Set iff a DC jump is present"
PAD             (bit 5) - "Set iff data are padding"
APOD            (bit 6) - "Set iff data are apodized/boundary"
STAT            (bit 7) - "Set iff telescope was stationary"
COM             (bit 8) - "Set iff data common-mode rejected"

Note that "BADDA" used to be called "BADSAM", and "COM" is also relatively new (common-mode rejection was using BADDA previously, which got confusing).

When you run the iterative map-maker you will now get some information about the quality. When the first iteration starts you will see something like this:


--- Size of the entire data array ------------------------------------------
bolos  : 1280
tslices: bnd:2000(0.2 min), map:150000(12.5 min), tot:152000(12.7 min)
Total samples: 194560000
--- Quality flagging statistics --------------------------------------------
 BADDA:   68400000 (35.16%),         450 bolos
BADBOL:   69312000 (35.62%),         456 bolos
 SPIKE:          0 ( 0.00%),
DCJUMP:          0 ( 0.00%),
   PAD:    2560000 ( 1.32%),        2000 tslices
  APOD:          0 ( 0.00%),           0 tslices
  STAT:    2560000 ( 1.32%),        2000 tslices
   COM:          0 ( 0.00%),
Total samples available for map:  123600000, 64.38% of max

The first part of this output is just telling you about the dimensions of the entire chunk that it is currently processing. With a single sub-array, the number of bolos is 1280 (ignoring for a moment how many are actually working), and "tslices" is the number of time slices (200 Hz samples). This second number is broken down: "bnd" is the number of samples in the boundary region -- these are the parts of the data that are either flagged as padding or apodization. "map" is all of the other samples -- this is the theoretical maximum amount of data that could go into the final map if all the bolometers were working, and there were no spikes, steps or other glitches. "tot" is simply the sum of the two. The "Total samples" line is then product of the number of bolometers with the total number of time slices in the array.

The second part of the report gives a breakdown of the quality flagging bits in the data array, after loading in the data and applying the flatfield, but before any of the pre-conditioning steps, or the first iteration of the model components. In the quality report the first column of numbers gives the total number of samples with those quality flags, followed by percentages of the total data array size in brackets. Since both BADDA and BADBOL apply to entire detectors, the third columns gives the total number of samples that have been flagged, divided by the number of time slices -- giving the numbers of bolometers. "PAD" and "STAT" flags are already set since padding and apodization of the data happens as they are loaded. Since all bolometers are flagged at a given time slice, the third column of numbers gives the number of samples divided by the number of bolometers, giving a number of time slices.

After each iteration, a similar report will be provided:


--- Quality flagging statistics --------------------------------------------
 BADDA:   68400000 (35.16%),         450 bolos  ,change          0 (+0.00%)
BADBOL:   80104000 (41.17%),         527 bolos  ,change   10792000 (+15.57%)
 SPIKE:          0 ( 0.00%),                    ,change          0 (+0.00%)
DCJUMP:    2113837 ( 1.09%),                    ,change    2113837 (+inf%)
   PAD:    2560000 ( 1.32%),        2000 tslices,change          0 (+0.00%)
  APOD:    5120000 ( 2.63%),        4000 tslices,change    5120000 (+inf%)
  STAT:    2560000 ( 1.32%),        2000 tslices,change          0 (+0.00%)
   COM:   10640000 ( 5.47%),                    ,change   10640000 (+inf%)
Total samples available for map:  107990669, 57.79% of max
     Change from last iteration:  -15609331, -12.63% of previous

At the end of the first iteration, this will report quality flags that were set both during pre-processing, and the first calculation of the model components. While some of the quality bits are static throughout the iterations (BADDA,PAD,APOD), the others can change each time. The fourth columns of numbers give the changes both as an absolute number of samples since the last iteration, and as a percentage of the number from the previous iteration (+infinity if the previous iteration didn't have any flagged). Here we see that the number of bad bolometers (BADBOL) has increased from 456 to 527 (15%) from the initial list of good bolometers. These extra detectors are almost exclusively rejected during calculation of the common-mode signal (they have a shape that is a significant outlier from all the other detectors). For this reason the increase in COM (in samples) is nearly identical to BADBOL. No spikes were flagged during this iteration (SPIKE), but about 1% of the data was flagged as having DC steps (in addition to the steps themselves, a small region around each step is also flagged). Finally, "STAT" identifies parts of the data where the telescope was not moving (threshold speed given by "flagstat" in the dimmconfig file). For these data it only flags the padding at the start and end of the data! (this is clearly not useful behaviour, so I will fix it to ignore these parts of the data in the future).

At the end of the individual quality stats is a summary for the entire array. "Total samples available for map" is the number of samples with no quality bits set, and it is also expressed as a percentage of the theoretical maximum discussed above for the "Size of the entire data array" at the start. The final line shows how this number has changed since the previous iteration.

In summary: watch the number of samples (and percentage) of the data going into the final map. If we typically have about 50% working bolometers, we should expect a number about that large going into the map at the end! If that number falls precipitously (or goes to 0), this report will help you to identify the problem. The usual culprits are DC step flagging, and common-mode rejection.

SMURF cookbook

2010-03-18T11:24:00.000-10:00

For the benefit of those not at the telecon this morning (or those that zoned out...)

The question was "is there a DR cookbook". The answer is yes, but it is about two months out of date. It should serve as an introduction, with the material in this blog providing an update. Of course, we will bring the cookbook up to date as things stabilise.

I have added a link to the SMURF cookbook in the "Useful links" sidebar section to the right. If you have a current Starlink installation you can also type "showme sc19" for the cookbook and "showme sun258" for the SMURF manual.

Flatfielding updates

2010-03-17T12:39:00.000-10:00

Summary: Flatfield ramps work really well and SMURF can now automatically handle them in the map-maker.

SCUBA-2 bolometers need to be calibrated to understand how they respond to varying signal coming from the sky and astronomical object. The original plan was to calibrate in the dark (shutter closed). The sequence goes something like:

Select a reference heater value, take a dark frame
Choose a new heater setting, take a dark frame
Take a dark frame at the reference heater value
Choose a different heater setting, take a dark frame
Take a dark frame at the reference heater value

and continue until you have covered a reasonable range of heater settings. As the heater is changed the bolometers read out a different current. Any drifts in the instrument are compensated by averaging the surrounding reference frames and subtracting. This means that you end up with a curve that goes through zero power at the reference heater value. In order to convert this to a flatfield you either fit a polynomial as a function of measured current (so that you can look up the power) or else use "TABLE" mode and do a linear interpolation between measurements either side of the measured current. The gradient of the curve (how the bolometer responds to changes in power) is the "responsivity" and is measured in amps per watt. The responsivity image can be calculated using the SMURF calcflat command.

When you open the shutter the idea is that you "heater track" to the sky. This involves you adjusting the power to the heater such that the sky power detected by the bolometer results in the same current being measured by the bolometer as it measured in the dark. We do this by looking at the signal from a set of tracking bolometers and assume that those bolometers are representative of the others on the array. In reality what happens is that about 80% of the bolometers do more or less read the same signal before and after opening the shutter but the other 20% are in a completely different place. This would not be a problem if the responsivity didn't change for those 20% but unfortunately it does. We have verified this by doing finely spaced pong maps on Mars covering a 6x6 arcmin area. This takes about 15 minutes but gives us a beam map of every single bolometer. Analysing the Mars images showed that the bolometers with the lowest responsivity also measure a very low integrated flux for Mars and so the calibration does change when the shutter is opened.

The solution for this was to change flatfielding to work on the sky rather than in the dark. This works in just the same way as previously, using reference sky measurements to compensate for drift, and the top plot in the figure shows a sky flatfield that is working pretty much perfectly. Finely-spaced maps of Mars confirm that all the bolometers are calibrated to within 10% with no drop off for the low responsivity bolometers.

At this point things were looking good but we still had the issue that the sky flat takes a few minutes and really has to be done every time you do a new setup and probably at least once an hour. They also are very dependent on observing conditions as could be seen on 20100310 and a few days before hand where the sky was terribly unstable despite brilliantly low opacity (0.03 CSO tau). The middle plot below shows a sky flat on 20100310 and it is immediately obvious that the sky is varying very fast and varying the power over a much larger range than the heater is adjusting for. This flatfield failed to calibrate any bolometers at all and we had to resort to dark flatfields to get a baseline calibration (with the associated worries described above).

We had known this was going to be an issue so in the early part of the year we had been modifying the acquisition system to do fast flatfield ramps. Rather than setting the heater, doing an observation, changing the heater, doing an observation we can now change the heater value at 200 Hz (currently we take 3 measurements at each setting though). On 20100223 we enabled sky flat field ramps at the start and end of every single mapping observation and a few days later we added it to focus, pointings and sky noise observations. The bottom plot shows the flatfield ramp for the observation that immediately followed the discrete sky flatfield shown in the middle plot. There is an issue with the very last ramp but the flatfielding software in SMURF had no problem calculating a flatfield for 850 bolometers (SMURF does compensate for drift in the reference heater values). The flatfield ramps are going to help enormously with calibration.

Actually using these flatfields in the map-maker took some work but yesterday I committed changes to SMURF so that flatfield ramps will be calculated and used when flatfielding data in the map-maker (and other SMURF commands). All you need to do is give all the files from an observation to SMURF and it will sort everything out.

I have updated the /stardev and /star rsync server in Hilo (64-bit and 32-bit). There is also a new nightly build available for OSX Snow Leopard 64bit in the usual place.

One final caveat, we have not yet calibrated the resistance of each bolometer relative to the nominal 2 ohms. We have taken data by looking at a blackbody source which should give us a way of tweaking the resistances. When this happens the flatfielding will change slightly and maps will need to be remade (although how critical that is will depend on how much we tweak the bolometers).

JLS/ACSIS DR telecon – 20100304

2010-03-09T09:31:00.000-10:00

Attendance: GAF, JH, CW, PR, JDF, TJ, FE, HST, ACC

Actions from last meeting:

JAC to make a more compact and readable QA report format and make this log available to observers/co-Is following nightly reduction via the OMP (as a downloadable file).
- update: there had been correspondence between JennyH and BradC and although there had been some work done, this item had not been completed before BradC left – ONGOING
JLS teams to provide feedback on what should go in parameter file (once JAC have made initial list). – ONGOING
BradC to work through SLS pipeline and QA requirements as provided to him (and available on SLS wiki)
- update: following BradC's departure, this will have to be picked up by his replacement (hopefully by next month) – ONGOING
For JLS teams to provide JAC with images/data/log of spikes when they come across them in their data.
- update: JAC had received some spike examples from SLS. These have been forwarded to DaveB but given his current priorities we may not see any movement on this for a few months yet (summer maybe?) – ONGOING
ChristineW to pass on NGS QA and moment map making scripts to JAC. – DONE
- ACTION: TimJ to look into whether BradC did code this into the pipeline?
AntonioC to organise the production of pipeline documentation
- update: some minor organisation of documentation but this has fallen down the list of priorities since the start of SCUBA-2 commissioning – ONGOING/STALLED
AntonioC to email coords asking them to provide email contacts of people who are reducing ACSIS data and are willing to act as ACSIS DR contacts. – DONE
JLS coords to provide email contacts of people who are reducing ACSIS data and are willing to act as survey DR contacts (send to FrossieE). – DONE

News from JAC

The software group took a big hit to its available effort when BradC left for new employment in February 2010. We are working on buying in some effort from an experienced ex-Starlink programmer to pick up on many of the tasks that Brad left. This person's effort should start becoming available to us from April 2010.
since the commissioning of SCUBA-2 started and the subsequent onset of SCUBA-2 Shared Risks Observing, JAC has been inundated in SCUBA-2 data
we can now re-reduce data at the JSA across different UT nights to create so-called "project products". This is a very important step and these products are essentially the Advanced Data Products (ADPs) that the JSA will serve.
FrossieE continues to submit data reduction jobs to the JSA but is experiencing problems as makecube has been timing out (in project mode) when cubes are very, very large

Report from GBS

a meeting of the GBS was help in Leiden in February 2010; discussions were held on the QA at that meeting
there are still missing QA steps in the pipeline reduction. For instance, bad scans are propagating through to final product, (bad scans = one or combination of high rms, many bad pixels, high variation across spectra)
- ACTION JAC: this is a high priority for the ACSIS pipeline at the JSA and JAC will look into this when effort materialises

Report from NGS

no issues reported
observing on the ACSIS portion of the survey has completed.
the team is reducing the data themselves via student effort; they will then want to upload the final products to the archive

Report from SLS

the team is focussing on getting all data reduced to then attack QA issue consistently
- ACTION: GaryF to direct FrossieE to SLS wiki where QA requirements are published
PaulR and HollyT have enough data now to compare pipeline reduced data with hand reduced data
- ACTION : PaulR to forward DR documentation to JAC – DONE

Next meeting May 6th.

Data processing at CADC

2010-03-01T18:14:00.000-10:00

Right now we are having trouble reducing SCUBA-2 data for the JCMT Science Archive (JSA).Briefly the problems seem to be:

The JSA wrapper is reporting success even in some cases of failure
Some maps take so long to make that MAKEMAP is timing out during long maps (this will also affect processing on private computers)
There is a bizarre NFS problem on the processing nodes that causes required perl modules to "disappear" when the pipeline is looking for them.

We're working through fixing these, so until then availability for products is a bit patchy. The good news is that processing throughput is great, so catching up with the backlog does not seem to be a problem.

Monitoring changes to the DR software

2010-02-26T09:59:00.000-10:00

Yesterday I gave instructions on how to retrieve the newest version of the Starlink software but it may not always be clear to people what changes are being made to the software to decide whether you want to get an update. The source code repositories for the Starlink software and ORAC-DR have RSS feeds that you can monitor in your standard news feed reader (e.g. Google Reader). Click on the RSS icon in the URL bar for the Starlink repository and the ORAC-DR repository.

SCUBA-2 DR status report

2010-02-25T13:38:00.000-10:00

We have been making some good progress with the data reduction software. Since hawaiki there have been updates to the flatfielding and an improvement in the code that detects working bolometers by comparing the common-mode signal. The down side of all this work is that for people getting their data from shared-risks observing they need to be using the cutting-edge version of SMURF and not the hawaiki version. The switch of flatfielding technique actually means that hawaiki will not generate a reasonable map at all using hawaiki SMURF.

To get the latest version of SMURF there are a number of options:

Use rsync to get a version for Centos5.4 (aka Scientific Linux 5.4 or RedHat Enterprise Linux 5.4) for either 32- or 64-bit linux. We attempt to keep these builds up to date. To see just how new you can run "$STARLINK_DIR/Perl/bin/starversion" and look at the date.
I can periodically make available 64-bit OSX Snow Leopard starlink releases. Check this directory periodically.
Download the source code and build it yourself (instructions on the Starlink home page)

As of 20100223 we have modified data acquisition to include a fast flatfield "ramp" at the start and end of every observation. These provide a direct calibration and should be more accurate than doing standalone flatfield observations as we have done previously (and continue to do). SMURF can not support these ramps at the moment but I am working on the issue. When the new code is working we will re-process the data from 20100223 at CADC.

Starlink release: Hawaiki (Deneb)

2010-01-20T20:31:00.000-10:00

Hawaiki has shipped!

http://starlink.jach.hawaii.edu/starlink/Hawaiki

Checkout the first (it sure won't be the last) version of the SCUBA-2 data reduction cookbook:

http://starlink.jach.hawaii.edu/docs/sc19.htx/sc19.html

Watching the data reduce

2009-11-05T14:58:00.000-10:00

Thanks to the folks at CADC, Dustin Jenkins in particular, JAC now has a really nice interface that allows us to monitor the jobs submitted to their Grid Engine, look for faults and browse the thumbnails of the products in order to spot problems - or just sit back and admire the results :-)

JLS/ACSIS DR telecon

2009-10-13T13:17:00.000-10:00

Attendance: ACC, BEC, TJ, FE, JVB, RP, CW, JdiF
Date: 15 Oct, 2009

Minutes:
1. Review of actions from previous meeting:

i. JAC will provide information on how to rsync the starlink releases to get latest patches/fixes. Information will also include for which operating system these patches/fixes are available. – DONE
JAC to make a more compact and readable QA report format and make this log available to observers/co-Is following nightly reduction via the OMP (as a downloadable file). – ONGOING
For SLS to provide JAC (ie Brad) with list of statistics and requirements for their QA, and also what they want for their reduction recipes to do. – ONGOING. Material has been received from SLS (also available on SLS wiki) and Brad is working through it.
For JLS teams to provide JAC with images/data/log of spikes when they come across them in their data. – ONGOING
ACC to organise the production of pipeline documentation. – STARTED. ONGOING.
ACC to poll for a date and time for next telecon and make these meeting notes available – DONE

2. News from JAC

Nightly reductions are now being carried out at CADC.
There have been tweaks to raster maps production such that the pipeline trims off the tassled edges of maps. Note that this is only for the maps and not the cubes.
We can now regrid to specific wcs coords (e.g galactic).
The infrastructure for controlling the pipeline via a parameter system and a config file is in place. Now need to decide which parameters are to be accessible before implementing the system. JAC will come up with an initial list and then invite feedback from users.
There will be two relevant newsletter articles in Autumn edition of JCMT Newsletter. One on JSA and another on the ORAC-DR pipeline.

ACTION : JLS teams to provide feedback on what should go in parameter file (once JAC have made initial list).

3. Updates from survey teams

GBS : working towards getting a consistent reduction on the data set. Need QA format finalised first.
NGS : data taking almost complete. Data is currently run through manual QA - flagging of bad data and receptors is critical and it’s subtlety requires human interaction. As the results of this is documented and the recorded state of the receptors isn’t going to change, we should be able to transfer that knowledge to an automated system in the future when we have a triggered re-reduction. Reduced products generated by NGS can be uploaded to the JSA. Improvements in moments maps: based on utilising SNR and noise maps. Good for data taken in different weather conditions.

ACTION : CW to pass on NGS QA and moment map making scripts to JAC.

iii. SLS : biggest issue is manpower. GAF has sent information to us to kick off the SLS pipeline QA process. Meeting @ JAC with Paul Ruffle tomorrow to discuss SLS DR issues.

4. Designated DR contacts
JAC would like to have email contacts with individuals in survey teams who are regularly looking at data and have the time and inclination to correspond and work with JAC between meetings. It is important that these individuals do have the time available (e.g. students, post-docs) to run tests and feedback ideas and improvements. Although these telecons are proving useful, progress is slow if we wait for monthly meetings to get feedback on more minor (yet potentially critical) issues. Everybody agreed that this is a good idea.

ACTION : ACC to email coords asking them to provide email contacts of people who are reducing ACSIS data and are willing to act as ACSIS DR contacts.

ACTION : JLS coords to provide email contacts of people who are reducing ACSIS data and are willing to act as survey DR contacts (send to FE).

Automated advanced processing at CADC

2009-10-02T12:14:00.000-10:00

As previously mentioned, the ORAC-DR data reduction pipeline could be run at CADC to generate basic nightly products. This processing used to be run at JCMT on a nightly basis, with the products transferred to CADC.

This processing is now being run at CADC on their processing system. Processing requests are made at 0900 HST every day, and nightly products will be available sometime after that (depending on the amount of processing needed -- scans take longer to reduce than jiggles).

In the near future, effort will be undertaken to process the backlog of ACSIS data.

SCUBA-2: baby steps

2009-09-28T13:10:00.000-10:00

As most of you will know, SCUBA-2 is now on the telescope being put through its paces. Given the high data rate of the instrument, the number of people working on it (in Canada and the UK as well as here at the JAC) and the upcoming early science call, we are racing to make raw data downloadable from the JSA by the beginning of November.

This is a challenge, both because the instrument systems themselves are in flux, and because everybody is so busy with commissioning. Still, thanks to the great efforts of Sharon and her team at CADC, we have coaxed some data into the JSA and searched for it using the test interface. It might not seem very exciting, but this has given most of the infrastructure a good workout, so it is a promising sign that we can meet our target.

ORAC-DR and Starlink on Twitter!

2009-09-04T09:35:00.000-10:00

Much to Tim and Frossie's chagrin, I've created two Twitter accounts for ORAC-DR and Starlink. I haven't completely sorted out what will be done with them at this time, but for now I'll probably use them to disseminate short tips and tricks for using ORAC-DR and Starlink software.

Follow them at http://twitter.com/oracdr and http://twitter.com/starlinksoft!

Thumbnails in search results

2009-09-03T08:32:00.000-10:00

As mentioned previously, the pipeline generates little thumbnails based on the representative image and the representative spectrum products. Now, CADC can show these thumbnails in search results. This hopefully will allow people to quickly identify search results of interest prior to downloading. Clicking on the thumbnails will launch the full preview of the representative image or spectrum as before.

If you have recent data, check it out. We hope to re-reduce the backlog at some point in order to generate products and thumbnails for older data too.

JLS DR telecon - 1st meeting

2009-09-02T08:02:00.000-10:00

Attendance: A. Chrysostomou, R. Tilanus, T. Jenness, R. Plume, M. van der Wiel, J. Di Francesco, G. Fuller, B. Cavanagh, H. Thomas, D. Johnstone, H. Roberts, D. Nutter, J. Hatchell, F. Economou

- initial discussion on whether we will have a SCUBA-2 pipeline ready. There will be something in place for shared risks but basic. More development will have to wait until we have all arrays in place as it is not worth sinking any effort into this at this time.

- some people are having issues getting the pipeline installed and the fact that there is a lack of documentation. If people/institutes are having issues installing (any) Starlink software, then please inform the JAC (stardev@jach.hawaii.edu) providing the relevant details.

ACTION 1: JAC will provide information on how to rsync the starlink releases to get latest patches/fixes. Information will also include for which operating system these patches/fixes are available.

DONE(!): Instructions are available on the starlink web site (http://starlink.jach.hawaii.edu/).
To download the most recent release go to: http://starlink.jach.hawaii.edu/starlink/Releases
To keep up to date with the latest fixes and patches go to:
http://starlink.jach.hawaii.edu/starlink/rsyncStarlink

- GAF requested for more statistics to be made available from the QA. GAF will follow up with specific request to Brad (see Action 3 below)

- it was clarified that the summit pipeline (during normal night-time observing) only runs basic QA on calibrations. After the end of observing, all data taken that night is re-reduced by a “nightly pipeline” which executes the full QA and advanced processing. The reduced data products which result from this are shipped to CADC and can be downloaded with (or without) the raw data in the normal way.

ACTION 2a: JAC to make QA log available to observers/co-Is following nightly reduction via the OMP (as a downloadable file).
ACTION 2b: JAC to make a more compact and readable QA report format.

ACTION 3: For SLS to provide JAC (ie Brad) with list of statistics and requirements for their QA, and also what they want for their reduction recipes to do.

- JH raised some existing issues from the GBS: flatfielding (striping) of early HARP data; some bad baselines are not being picked up by QA; although not as prevalent as in older data, spikes are not trapped by the QA; an investigation is needed on how the gridding should best be done

+++ the flatfielding problem is on Brad’s worklist

+++ we need more feedback from the teams on which bad baselines are not been filtered out

+++ de-spiking data is not a problem that JAC has been able to tackle as yet. Part of the issue is that these do not seem to be as prevalent in data any more and observers (PI as well as JLS) are not reporting the issue any longer. GAF reported that spikes are still present but at a small level, which is an issue for SLS who are looking for weak, narrow lines.

ACTION 4: For JLS teams to provide JAC with images/data/log of spikes when they come across them in their data.

- RPT raised a few issues from the NGLS:

+++ need ability to baseline fit both wide and narrow lines in same data set

+++ need ability to restrict e.g. moments analysis to known velocity range.

+++ QA generally fails for (at least) early NGLS data. Will need to investigate this more but need an easy means to switch off in recipes. This is easy in the main recipe, but less so in the advanced interative part.

- there is a blog available for data reduction and pipeline activities (you’re problably looking at it right now!): http://pipelinesandarchives.blogspot.com/

- the issue of making the pipeline more controllable through a config file to set parameters was discussed. TJ announced that he is developing infrastructure so that the pipeline can be parameterised

- ACC received several emails prior to the meeting. A common theme was the lack of documentation explaining what the pipeline does to data, and how to use the pipeline. JH repeated this concern at the meeting.

ACTION 5: ACC took an action following the close of the meeting to organise the production of pipeline documentation. These will probably take the form of a detailed account of what the pipeline does, and a separate cookbook which explains how to run the pipeline with the different options available.

ACTION 6: ACC to poll for a date and time for next telecon and make these meeting notes available.

Starlink Software Collection - Nanahope (Pollux) version released

2009-07-27T14:30:00.000-10:00

The Nanahope version of the Starlink Software Collection was just released and can be downloaded from here. Highlights include:

GAIA can now visualise 2-D and 3-D clumps created by the CUPID package.

GAIA now has full support for the Virtual Observatory and has been modified to support the SAMP protocol to enable it to communicate with other VO tools.

Automated provenance propagation can now track HISTORY information in addition to provenance. The PROVSHOW command can now list the history of all of the parents in the processing history. HISLIST (and NDF history propagation) has not been changed and still only examines the history of a single path through the processing.

The software can now be built with gfortran 4.4.

Database Replication outage

2009-07-13T11:12:00.000-10:00

Over the weekend we lost database replication to CADC. Until we can sync up the tables we are unable to transfer any raw data to CADC (since CADC only accept files that their systems know about) so data will not be retrievable from the weekend. I'll post an update when transfers are enabled again.

UPDATE: Replication server crashed on Friday night. It's now back up and the tables have been synced with CADC. Transfers have been restarted.

Hierarchical history for NDFs

2009-07-09T22:50:00.000-10:00

The recording of processing history has been part of the NDF library for many years. When an application uses one or more input NDFs to create an output NDF, the NDF library creates a record of the application and its parameter values, and stores this record in the output NDF. It also copies all the history information from the "primary" (usually the first) input NDF into the output NDF.

Whilst it was recognised at the time that it would be nice to copy history from all input NDFs, the exponential growth of history information this could cause was seen to be prohibitive. But 16 years is a long time and we typically now have far greater computing resources. So we've taken the plunge and changed things so that history from all input NDFs is copied into the output NDF. However, to preserve backward compatibility, the new facilities are provided by the provenance routines in the NDG library - the NDF library itself remains unchanged.

This means that applications such as KAPPA:HISLIST, etc, that use the NDF library directly to manipulate history information are unchanged. Instead, the extended history information is stored in the PROVENANCE extension of each NDF, and can be examined using the KAPPA PROVSHOW command. Since there can be quite a lot of history information, it is not shown by default - set the new HISTORY parameter to "YES" when running PROVSHOW to change this default behaviour. Needless to say, NDFs created before these changes were made will not contain any extended history.

A common use for this extended history will be finding the value used for a particular parameter when a selected ancestor was created. We're toying with the idea of a GUI that would make this sort of thing easier by allowing an NDF's "family tree" to be navigated and searched, but for the moment the best thing is probably to use grep on the output of PROVSHOW.

GAIA goes all clumpy

2009-07-08T06:34:00.000-10:00

In the next release GAIA will display CUPID catalogues and masks so you can inspect your clumps in all their detail. This all works in 2 and 3D, which you can see in more detail at on the GAIA support site.

CADC network outage

2009-07-07T22:19:00.000-10:00

From John Ouellette at CADC:

The CADC will be undergoing extensive maintenance from July 18th 0800 PDT to July 19th 1800 PDT. All CADC services, including user access and etransfer, will be unavailable during this period.

During the outage, users will be redirected to a web page stating the reason for the outage and, if possible, we will provide a status update during the work.

CUPID now creates STC-S polygons

2009-06-09T22:08:00.000-10:00

I've added an option to CUPID:FINDCLUMPS to allow it to create an STC-S description (either polygonal or elliptical) for each clump it finds, and add them into the output clump catalogue. I've also modified KAPPA:LISTSHOW so that it can display the STC-S shapes over a displayed image. Below is an example. In the first image, the greyscale (and contours) are the data, and the red lines are the polygonal clump outlines. They overlap slightly in some cases because each polygon is only allowed to have up to 15 vertices, and so only approximates the clump pixel mask.

The next image is the corresponding pixel mask (each colour represents the pixels assigned to a single clump).

Pipeline now generates thumbnails

2009-06-03T11:24:00.000-10:00

The data reduction pipeline now automatically generates PNG thumbnails of rimg and rsp files. These thumbnails are generated in three different sizes, 64x64, 256x256, and 1024x1024. Exif information is also written to these thumbnails, embedding the RA, Dec, source name, orientation, and pixel scale. Only the astro: namespace is currently used (see this ROE page for more information).

Here's a (rather boring) example from last night:

In the near future they will be sent to CADC for automatic ingest.

Pipeline running on CADC grid engine

2009-05-12T08:27:00.001-10:00

Last week I (Tim J) visited CADC to work on integrating ORAC-DR into the CADC grid engine system. This involved making sure that the pipeline wrapper interfaced properly with CADC and that the data retrieval and data capture routines were given the correct inputs.

On Wednesday 6th May we were successful in running four jobs in parallel on the compute cluster. This is a terrific result and leads the way to being able to do the night processing on CADC in short order and then to follow up with processing of project coadds. In the next few weeks I will be working on the code that will query the data archive and submit job requests to grid engine.

It also means that in principal survey teams could request that jobs be submitted to grid engine to make use of the processing nodes.

CADC network outage

2009-05-11T10:29:00.000-10:00

CADC and therefore JCMT data retrievals will be off the air on Wednesday afternoon (HST) due to scheduled network maintenance. If you are having trouble getting through, just try again later.

Near minutely raw data movement

2008-12-19T14:58:00.000-10:00

Raw observation data is being put both in jcmt database & CADC staging area not much later when it arrives.

A new program[0] endlessly checks (currently) about every minute to see if there are any observations to process. If there are, the database is fed, followed by symbolic link creations for CADC's consumption. This should help avoid massive data transfers to CADC twice a day. Note that previously involved programs will keep running concurrently until everybody involved is satisfied that raw data is being entered and/or moved as desired.

All this started yesterday slightly wet, cloudy Hawaiian evening.

[0] enterdata-cadc-copy.pl is a wrapper around JSA::EnterData & JSA::CADC_Copy modules, which were respectively generated from jcmtenterdata.pl & cadcopy.pl.

QA-enabled pipeline released in Hilo

2008-12-04T16:30:00.000-10:00

ORAC-DR has been updated in Hilo to include quality-assurance testing. Based on a number of QA tests, observations are given a pass/questionable/fail status. QA is automatically done on all science observations, and survey-specific QA parameters can be given.

This version will eventually be released to the summit, where it will give telescope operators feedback on which surveys are suitable to do, along with enhancing the JCMT Science Archive pipeline.

SCUBA-2 DR pipeline

2008-06-27T13:40:00.000-10:00

A belated announcement that the SCUBA-2 data reduction pipeline passed its "lab acceptance" earlier this month. Full report at http://docs.jach.hawaii.edu/JCMT/SC2/SOF/PM210/04/sc2_sof_pm210_04.pdf

Initial results of "better" ORAC-DR reduction

2008-04-01T15:52:00.000-10:00

As previously mentioned, ORAC-DR is improving how ACSIS data are reduced. To show the progress between "summit" and "better":

integrated intensity map, group coadd: summit better

integrated intensity map, single observation: summit better

intensity-weighted velocity map, group coadd: summit better

A few notes:

By "summit" pipeline I mean the pipeline currently running at the summit. This pipeline will be replaced by an "improved" pipeline pending JCMT support scientist approval. The "improved" pipeline will not be run at CADC, they will run the "better" pipeline that created the "better" images linked above.

The "group" summit integrated intensity map is not generated by the pipeline, it's created by manually running wcsmosaic to mosaic together the individual baselined cubes (the _reduced CADC products), then collapsing over the entire frequency range. This is how the summit pipeline would create those files, though.

Ditto for the group summit velocity map, except the pipeline wouldn't even create those in the first place, as it doesn't know which velocity ranges to collapse over to get a proper velocity map. This example is just done by naively collapsing over the entire velocity range. The "better" pipeline automatically finds these regions and creates velocity maps -- not only for the coadded group cube, but also for individual observation cubes.

The difference between the "better" pipeline (which is what will be running at CADC) and the "improved" pipeline (which is what will be running at the summit) is very small for this given dataset.

CUPID ClumpFind and backgrounds

2008-03-26T23:39:00.000-10:00

Jenny Hatchell has been comparing the CUPID implementation of the ClumpFind algorithm with the IDL implementation by Jonathan Williams. The IDL version differs in one or two significant respects from the original algorithm published in ApJ, and so CUPID provide a switch that selects either the published algorithm or the IDL algorithm. If the IDL algorithm is selected, Jenny finds that the IDL and CUPID implementations allocate exactly the same pixels to each clump. Good news. And more good news is that the CUPID implementation is at least an order of magnitude faster than the IDL implementation.

However, Jenny noted that the clump sizes reported by CUPID were not the same as those reported by IDL. Both implementations use the RMS displacement of each pixel centre from the clump centroid as the clump size, where each pixel is weighted by the corresponding pixel data value. So in principle they should produce the same values. The difference turns out to be caused by the fact that CUPID removes a background level from each clump before using the pixel values to weight the displacements. IDL , on the other hand, uses the full pixel values without subtracting any background. Thus, increasing the background level under a clump will produce no change in the clump sizes reported by CUPID. IDL, however, will report larger clump sizes due to the greater relative emphasis put on the outer edges of the clump.

So should a background be subtracted or not? Having the reported clump size depend on the background level seems an undesirable feature to me. But if you want to compare CUPID results with other systems (e.g. the IDL ClumpFind in this case) that do not subtract a background, you CUPID also needs to retain the bacground level to get a meaningful comparison. Consequently, I've added a parameter to CUPID:FINDCLUMPS to select whether or not to subtract the background before calculating clump sizes. The default is for the background to be subtracted unless CUPID is emulating the IDL algorithm (as indicated by the ClumpFind.IDLAlg configuration parameter).

If the background is retained in CUPID, Jenny found that the CUPID and IDL clump sizes match to within half a percent. So things look OK.

David

Preliminary Wrapper script released

2008-02-29T17:25:00.001-10:00

A first stab at a wrapper script (jsawrapdr) has been released to users in Hilo. It matches the interface specified by CADC.

What it does:

retrieve data files from the supplied list
convert them to NDF if required
determines the correct ORAC-DR instrument name based on the data
checks that PRODUCT information matches for all files
determines whether to run ORAC-DR or PiCARD
converts products back to FITS

What it doesn't do yet:

Provenance is not quite correct. It is possible to refer to a parent that will not be archived.
There is no standardised approach to logging Standard Output and Standard Error
dpCapture does not automatically copy products to the CADC transfer directory.

It's enough to get us started.

CADC data transfers now working again

2008-02-29T13:44:00.000-10:00

Data transfers to CADC are functioning again. We've had real problems reconfiguring replication to CADC and to our standby server (which use different techniques) but now everything seems to be working. Now that CADC have headers from recent observations they will again start accepting our raw data. Transfers have been initiated and are currently complete to 20080215 (there are quite a few files to transfer). Data retrieval requests from users via the OMP will shortly be redirected back to CADC.

ORAC-DR: CADC+batch mode

2008-02-25T16:11:00.000-10:00

A brief mind-dump on the processing steps ORAC-DR will probably take at CADC when run in batch mode:

_cube files created. These are then forgotten about by ORAC-DR but are stored by CADC.
Run initial steps of Remo's script on time-series data. Removes any gross time-series signal through collapsing and rudimentary linear baselining.
Run MAKECUBE using every member observation of a Group, creating tiles.
Run remainder of Remo's script on each tile, which uses a combination of smoothing and CUPID to create baseline region masks and to remove baselines.
Take the baseline region mask from the previous step along with the original input time-series data, and throw them through UNMAKECUBE. This will create time-series masks.
Apply the time-series masks to the original time-series data.
Run MFITTREND with a high-order polynomial (or spline, or whatever) on the masked time-series data. These cubes shouldn't have any signal and should be pure baseline.
Subtract the baselines determined in the previous step from the original time-series data.
Run MAKECUBE on the baselined time-series data for each observation to create the _reduced / _rimg / _rsp files.
If necessary, WCSMOSAIC the _reduced files for each observation to create an "even better" group, which can then be used to determine a better mask and then possibly iterate through the UNMAKECUBE to _reduced generation steps.

The Wrapper

2008-02-25T14:03:00.000-10:00

Background: there is a system called, imaginatively, the wrapper. Its purpose is to wrap the data processing specifics so as to present a generic interface to the CADC data processing system that is under development.

The wrapper is on TimJ's to-do list, and so is at the mercy of his higher priority SCUBA-2 work. In an attempt to push something out to CADC before working on the SCUBA-2 translator, he is writing a prototype with the following functionality:

has a stub dpRetrieve, emulating the system that will eventually fetch the data needed from the CADC database
examining the data to determine whether it is raw or already a product
converts any FITS files to NDF
runs ORAC-DR or PiCARD as appropriate given the above information
converts any NDF products back to CADC-compliant FITS
calls a stub dpCapture (the real dpCapture imports any products into the CADC system)

The main problem that stops this from being more than a prototype is the provenance system. In our NDF based systems provenance is a time series - file A turns into file B which turns into file C.... eventually resulting in file E, the final product. So the provenance looks like this: A, B, C, D, E. In the CADC system, provenance is the nearest parent existing in the archive. So, if only A, B, D and E exist in the database (because C happens to be an intermediate file of no lasting importance) the provenance for E is D, but the provenance for D is B. Therefore, the wrapper has to make sure that at the end of any processing the provenance is correctly fixed to display only parents existing in the CADC archive.

The intended solution for this is for DavidB to commit some NDG patches to allow TimJ/the wrapper to remove C (in the previous example) from the provenance. Also, the wrapper needs to rename A, B, D and E to the CADC naming convention so they can be matched to entries in the archive.

So as not to hold other parts of the project up, the intent is for the prototype to be delivered to CADC in the next few days without this provenance-related functionality, and to come back and fix this when the SCUBA-2 translator work allows.

OMP to CADC connection is down

2008-02-13T10:30:00.001-10:00

At the tail end of last week we had a hardware failure at the summit with our primary OMP database server. We switched to the new Sybase 15 64-bit servers but they have not been configured correctly to replicate the JCMT header table to CADC (full database replication is working to the backup 64-bit server in Hilo). Until we get the CADC replication up and running there will be no transfers of raw data to CADC. This is because CADC validates transfers against it's copy of the header table and rejects observations that are unknown to them.

We hope to have replication running by early next week but in the mean time we have reconfigured OMP data retrievals to serve the raw data files from JAC.

specx2ndf now creates variance

2008-02-11T02:51:00.000-10:00

Forget to mention that I've modified specx2ndf so that it creates a Variance component in the output NDF based on the Tsys value in the specx file. The variance values in the output NDF are constant since each specx file seems to contain only a single Tsys value.

Processing 3D cubes with FFCLEAN

2008-02-04T06:39:00.000-10:00

I've just modified kappa:ffclean so that it can:
1) process 3D cubes. It will do this either by processing the cubes as a set of independent 1D spectra, or as a set of independent 2D images (see new parameter AXES)
2) store the calculated noise level in the output variance array (see new parameter GENVAR)

This was motivated by my experiments with the new smurf:unmakecube command as a means of getting an estimate of the noise level in each residual spectrum.

Creating artifical time series from a sky cube

2008-02-01T06:53:00.001-10:00

I've just added a new command to smurf called UNMAKECUBE. It takes a (ra,dec,spec) sky cube and generates artifical time series by resampling the sky cube at the detector positions in a set of existing reference time series NDFs. It's a sort of inverse to MAKECUBE. It should be useful for investigating baselines, and for iterative map making. I'm currently playing around with it, using data from Christine Wilson.

Time accounting for JLS

2008-01-31T09:58:00.000-10:00

Another issue for JLS observing is with time accounting. How to account in the OMP the time for data that has been deemed unacceptable. These data do not get charged to the surveys so if the data were initially QUESTIONABLE retroactive action will need to be taken to correctly account for that time (notwithstanding the issue of shared calibrations...see later). So that we can track how much time has been REJECTed by each survey, there shall be a special project code (eg. MJLSG00) for each which we will use to charge REJECT observations to.

The idea is that when the obslog flag changes:

1. an email is triggered to ACC and PIs notifying the change

2. if the change is to REJECT then the release date is automatically changed to TODAY (or equivalent)

3. ACC runs up nightrep for the night in question and changes the time accounting accordingly.

4. changes to flags are propagated to CADC

In a situation where calibrations are provided by the observatory this system should work flawlessly (and a tool which takes care of the time accounting automatically would also be feasible). However, in the current situation where calibrations are shared amongst the projects, it is difficult to do the time accounting properly in this scheme as it is not immediately obvious how much calibration time should be taken with the REJECTed observation. It would have to be recalculated.

How to properly deal with BAD/QUESTIONABLE data within the JLS

2008-01-31T09:41:00.000-10:00

The problem stems from what happens to questionable data which remains in some form of limbo until its status is deemed to be GOOD or BAD. Setting data to BAD is an issue in itself as such data are usually BAD because of a fault. However, the Legacy nature of JLS means that some data will be deemed to be unacceptable for the surveys and so should not be processed into Advanced Data Products. These data are not intrinsically BAD and so the plan is for them to be immediately released to the public.

We resolve this issue by having a new quality parameter in obslog - REJECT. Definitions:

GOOD: data is good and is processed by the pipeline

QUESTIONABLE: data may have problems with it - a human needs to look and make a decision on its quality

BAD: data is not good, do not process

REJECT: data does not meet agreed standards for survey

We don’t expect to be using the REJECT flag during normal PI observations. Furthermore, it is expected that with working QA in the survey pipelines the number of REJECT and QUESTIONABLEs will be small (but there will be enough, especially at the beginning as we're bedding the system, that we need to deal with them appropriately).

The following table summarises what happens to data with these flags:

_cube _reduced group proprietary

GOOD Y Y Y Y

BAD Y N N N

QUEST. Y Y N Y

REJECT Y Y N N

ADP charged VO/master product

GOOD Y Y Y

BAD N N N

QUEST. N Y N

REJECT N N Y

(N.B: QUESTIONABLE data should not be combined into the public VO product as those data are in an undefined quantum state and until their wave functions have collapsed into either of GOOD, BAD or REJECT you don’t know what to do with them)

ACSIS DR and hybrid observations

2008-01-30T15:45:00.001-10:00

In the version of the ACSIS pipeline that will ship in the upcoming Starlink "humu" release subsystems from ACSIS are processed independently of each other. This is fine for true multi-subsystem observations such as simultaneous C18O and 13CO modes but is not the correct thing to do when a true hybrid mode is observed. This was not an issue when we relied on the "real-time" DR to do the sub-band merging but since the decision was made to configure 500MHz and 2GHz observations as pseudo multi subsystem observations (albeit with shared reference pixel so that the channels are aligned) it has become obvious that ORAC-DR needed to be modified.

This week I modified the internals of ORAC-DR to allow it to recognize hybrid observations and treat them as a single "frame" (in ORAC-DR speak). It was a little more involved than expected because there was no merging of header information in the pipeline outside of bespoke implementations in UKIRT spectroscopy and SCUBA-2 classes. I moved some header parsing code into a base class and removed the special code from UKIRT (SCUBA-2 is still special). This required some minor changes to all the header translation code but was worth it since all instruments can make use of the header merging.

Next step is to actually use this information to merge the sub-bands. One minor caveat in all this is that ORAC-DR will still not try to merge multi-subsystem observations simply because subsystems overlap. If that is required (for example for the SLS) then we need to be told the requirements.

DB replication to CADC

2008-01-28T12:48:00.000-10:00

Some confusion as to where we are with getting replication from the ASE 15 system to CADC. Phone call with anubhav, timj and isabella tomorrow (Tuesday) 12pm HST.

scubadev

2008-01-28T10:51:00.000-10:00

Scubadev's mysterious lack of speed continues to be a concern, especially since it is a baseline spec for the summit DR computers. Had a chat with Tim where we concluded that after the Starlink release is out of the way (1-2 weeks) we will take it down and turn it into a standalone gentoo box, then a standalone CentOS box. This should allow us to narrow down whether the problem is related to hardware, OS, or the integration in to the network.

OMP ACSIS data retrievals

2008-01-28T10:31:00.001-10:00

Last week we had a strange problem where data could be retrieved from the OMP from November onwards but between Feebruary and October 2007 data retrievals failed because the OMP sent the wrong filenames to CADC. Turned out that the new test database (running Sybase 15) had been loaded with data up to end of October and that was triggering a new logic path through the OMP. Usually, the OMP failed to find any entries in the database and fell back to looking on the data disk for raw data. This always works and always finds the right files. When rows are found in the database there is no need to look on disk (the database is much faster than looking at files) and once the test database was initialised the DB lookups were working properly. The only problem was that the query to the ACSIS database did not return the filename information from the FILES table and therefore the OMP was forced to guess the filename. For data taken since we renumbered the subystem numbers the guess was wrong and CADC were asked to serve files that didn't exist.

I fixed the problem last week and now retrievals work with database and file lookup. I was able to cleanup quite a lot of code in the process and the Astro::FITS::HdrTrans module was made a little cleverer and can now tell the difference between a database result and a header read from a file. Apologies to people who experienced retrieval problems over the past 2 weeks.

ACSIS wrapper script

2008-01-25T13:25:00.000-10:00

Short telecon arranged for Monday 12:30 HST (after regular JSA 12:00 HST CADC team meeting) so that Tim can put forward his ideas (mostly to Sharon) about the ACSIS processing wrapping script. There is some interest in testing this by February 7th, which will be tight.