SCUBA-2 at the ADASS conference

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

Adding fake sources to real data

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

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

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

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

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

Map-based Despiking

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

ast.mapspike = 10

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

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



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

ast.mapspike = 0

spikethresh = 10

spikebox = 50

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



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

Flagmaps and Shortmaps

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

flagmap = (BADBOL,SPIKE,DCJUMP)

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



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

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

shortmap = -1

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



Choosing frequency-domain filters based on angular scale

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

filt_edgehigh = 0.3

flt.filt_edgehigh = 0.3

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

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

filt_edge_largescale = 600.

450.filt_edge_smallscale = 2.

850.filt_edge_smallscale = 4.

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

850.flt.filt_edge_largescale=300

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

flt

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

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


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

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

flagslow = 30

450.flagfast = 600

850.flagfast = 980

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

Suppressing large-scale noise with zero-masks

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

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

The first is to define a circular masked region:

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

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

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

ast.zero_snr = SNR

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

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

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

Signal = ASTronomical + COMmon + FiLTered

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

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

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



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



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



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

SCUBA-2 Pointings in the science archive

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

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

Spatial noise vs noise maps

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

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

Larger image here

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

Larger image here

Larger image here

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

S2SRO data processing notes

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

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

The S2SRO data maps

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

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

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

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


Calibration

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

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

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


Noise Estimates

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

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

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

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

Processing date now visible in JSA web UI

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

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

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

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)

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

Version 1.1 of the SMURF SCUBA-2 Data Reduction Cookbook  (or SC21, 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 SC21 from the sidebar of this blog.

Updated to change reference from SC19 to SC21

JSA FAQ: Product grouping types

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:

1. Simple: This is the most processed state of a single observation
2. Night: This is the product resulting from all observations taken in one night on the same field
3. 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
4. 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!

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:

1. 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.
2. 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

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)

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

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

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:

1. 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).
2. 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.
3. 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) .

1. Run your selected calibration observation through the mapmaker using the same dimmconfig as your science data used.
2. 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. 
3. 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

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

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

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

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

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)?

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?

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:

1. 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.
2. 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

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