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.