Chapter 6
Tailoring Your Reduction

 6.1 Adding and amending parameters
 6.2 Writing out models & intermediate maps
 6.3 Large-scale filtering
 6.4 Fitting COM for each sub-array
 6.5 Flagging bad data
 6.6 Using external masks
 6.7 Skyloop
 6.8 Troubleshooting

6.1 Adding and amending parameters

The configuration file dimmconfig_jsa_generic.lis is a good configuration file to use as a first pass for reducing any data, and is the default configuration used by the pipeline’s REDUCE_SCAN recipe.

You can create your own personalised configuration file from scratch or copy one of the provided ones to your local directory and edit it.

The first line of each specialised configuration file is often the path to another configuration file — known as the parent configuration file. The new configuration file inherits all the parameter values defined by the parent configuration file, and then goes on to specify additional parameter settings which may supplement or over-ride those defined in the parent configuration file. Parameters that are not set to a specific value by either configuration file assume the default values listed in Appendix H1.

Often you will want to use one of the pre-defined configuration files included within Starlink tree as the parent configuration. Remember that any parameters appearing in your configuration file automatically override the values supplied by the parent file. As an example, consider the following text file “myconf”:

  % cat myconf
  
  #  This is an example configuration file.
  ^$STARLINK_DIR/share/smurf/dimmconfig_jsa_generic.lis
  
  numiter=-100
  maptol=0.01
  

The first thing to note is that blank lines and comment lines (i.e. lines beginning with the hash character “#”) are ignored and can be used to document your configuration file. Next note that this example file uses $STARLINK_DIR/share/smurf/dimmconfig_jsa_generic.lis as its parent (the path to the parent file must be preceded by an up-caret (̂) character). Thus, all the parameter values set by $STARLINK_DIR/share/smurf/dimmconfig_jsa_generic.lis are first read, before adding in other parameter settings. In this case, values are assigned to maptol and numiter, over-riding the values provided by the parent file.

You can use more than one parent file if required (each one a separate line, and preceded by an up-caret). Each specified parent will be read in turn, with parameter settings read from later ones having priority over those read from earlier ones.

You can also add or amend parameters by listing them directly on the command line. They are appended to the configuration file name as a comma separated list as shown in the example below. Be sure to include all the necessary quotation marks.

  % makemap in=’s8*.sdf’ out=850map \
         config=’"^dimmconfig_jsa_generic.lis,numiter=-50,exportndf=(flt,noi),itermap=1"’

For full details of all the possible ways of specifying groups of parameter values, see Section “Specifying Groups of Objects” in “Starlink User Note 95.

What parameters can be changed?
Appendix H lists the parameters that are more likely to be of interest to you when creating your own configuration files. You can also change any of the other more esoteric parameters not included in that list — see Appendix SUN/258 within Starlink User Note 258 for a full list – but we do not advise this.

Tip:
If some feature is switched on either by default or within the parent configuration file, it can be switched off if required by assigning a suitable value to the corresponding parameter within your own configuration file. The value needed to do this will be given in the parameter description in Appendix H — for instance it may be <undef>, or 0, or some other special value.

Note: any parameter can be made wavelength dependent by adding the prefix 450. or 850., e.g. flt_edge_largescale applies to both 450 μm and 850 μm whilst 450.flt_edge_largescale applies to 450 μm only. Be aware that if both are specified, unqualified values (no prefix) take priority over qualified values.

6.2 Writing out models & intermediate maps

itermaps
Setting the parameter itermap = 1 writes out the map containing the astronomical signal after each iteration. Setting itermap = 2 adds the QUALITY component. These can be visually inspected with

  % gaia 850map.more.smurf.itermaps

to help determine an appropriate number of iterations. Alternatively, you can view several itermaps simultaneously side-by-side (e.g. Figure 5.1) using Kappa as described in Section 5.3.2.

Viewing itermaps is useful when a fixed number of iterations have been requested (i.e. a positive value for numiter) and the map solution diverges before they have completed. See also Section 5.3.2 for how to view these maps whilst makemap is still running.

shortmaps
If the parameter shortmaps is non-zero, a map is made from every group of adjacent time-slices (as specified by the parameter). These are stored as an NDF extension and can be viewed Gaia.


Viewing ITERMAPs

  % stackframes map.more.smurf.itermaps \
  sort=false map_itermaps

Stack the individual itermaps into a single cube (Kappa PASTE can also be used).

pict

The output map map_itermaps can be opened with Gaia. The data used in this example is the Galactic map reduced in Section 7.2. The Spectral plot window shows the value for a single pixel and the three chunks are easily identified. You can select the Animation tab in the Display image sections window and click Start to loop through the itermaps for each iteration. The ‘movie’ will appear in the main Gaia window.

pict pict pict

These windows show the itermaps map at 1, 10, and 30 iterations. A specific iteration can be selected using the Index of plane slider on the Display image sections window.


Figure 6.1: Example using the Smurf command stackframes and Gaia to view the ‘itermap’ for each iteration.


  % gaia 850map.more.smurf.shortmaps

You can view the shortmaps and itermaps more conveniently by stacking them into a single cube using the Smurf command stackframes. This cube can then be viewed as a ‘movie’ with Gaia, using the animation option to loop through the itermaps. See Figure 6.1 for instructions.

exportmodel
This parameter has been discussed in Section 9.8 and allowed you to see the model that was fit for each component specified by the modelorder parameter.

6.3 Large-scale filtering

Some of the most important parameters to experiment with are the filtering options. By default, no filtering is applied during the pre-processing stage2 (see filt_edge_largescale and filt_edge_smallscale). A high-pass filter is used during the iterative stage if moelorder includes FLT3, and selected a suitable value for the filter size is crucial for maps containing extended emission.

The maximum spatial scale of structure that can be recovered by the map-maker is determined by the scanning speed and frequency cut applied to the data:

speed[arcsec/s] frequency cut[Hz] = scale size[arcsec] (6.1)

The default (i.e. if no configuration is supplied) filter sizes in arc-seconds are:

450.flt.filt_edge_largescale = 200
850.flt.filt_edge_largescale = 480.

To make your life easier, these parameters allow you to specify the filter limits in terms of spatial scale in arc-seconds—in this case 480 arcsec at 850 μm and 200 arcsec at 450 μm. For example, at 850 μm, recovering scales of 480 arcsec at a scan speed of 600 arcsec/sec (default for a 1 degree pong) corresponds to a frequency of 1.25 Hz.

Choosing a high-pass filter is especially important for the recovery of extended emission. The dimmconfig_bright_extended configuration file sets flt.filt_edge_largescale = 480 for both 450 μm and 850 μm. Be aware that increasing filter sizes decreases the flatness of your background. A compromise must be made between extended structure and the flatness of your map. See Figure 6.2 for an illustration of the effect of flt.filt_edge_largescale on your map.

The scanning speeds are fixed for a given observing mode; you can find out the speed at which your data were taken from the SCAN_VEL keyword in the FITS header (see Section 9.2).

Flattening the background
There is an option to reduce the noise in your background introduced by setting a high value for the large-scale filter. The parameter flt.filt_edge_largescale_last filters the regions outside your AST mask (see Section 3.5.1) on a shorter scale for the last iteration only, thereby producing a much flatter background. Note that this is probably a bad idea if you intend to co-add several observations, as it removes potentially real structure in the background regions that could otherwise be recovered by co-adding several observations. In general, use of flt.filt_edge_largescale_last should be seen as a cosmetic enhancement, since it results in differing filter sizes being used inside and outside the AST mask.

Also note that when using flt.filt_edge_largescale_last the variances stored in the final map are from the penultimate iteration in order to avoid using the artificially reduced variances created on the last iteration.


pict   pict
Figure 6.2: Highlighting the effects of high-pass filtering on your map. (Left) Map made with 850.flt.filt_edge_largescale300. (Right) Map made with 850.flt.filt_edge_largescale1000. All other configuration parameters remain the same.


6.4 Fitting COM for each sub-array

A useful option to improve the flatness of your maps is to fit the COM model independently for each sub-array. This is particularly effective if you find you have one sub-array noisier than the others.

This comes with the warning however that you will lose information on any scales larger than the area covered by a single sub-array. It is therefore not recommended if you have very large-scale extended structure.

To initialise this option set com.perarray = 1.

6.5 Flagging bad data

Bolometers that have higher noise levels are down-weighted when forming the map. By default, a separate noise estimate is created for each 15 second section4 of each bolometer time-stream, but an alternative length can be specified using parameter noi.box_size.

Dividing each bolometer up into sections helps to remove ‘scuffs’ or other noise artefacts you might see in your error map due to a sub-array (or arrays) temporarily jumping to a higher noise state. As the section length tends to 1 second we find some of the source signal being down-weighted. Higher than this 15 seconds and the map-maker becomes less sensitive to the higher noise states.

Other parameters you may want to try include flagfast and flagslow. You may find that setting flagfast to less than the default of 1000 arcsec/sec will help reduce the effect of any ‘smearing’ of sources (and of noise) in maps, while setting flagslow greater than the default of 300 arcsec/sec helps to flatten the edges of maps. To determine reasonable values for your data-set you should do jcmtstate2cat and view the scan speed using Topcat. See Section 9.3 for details.

6.6 Using external masks

For an introduction to the purpose and effects of masking, see Section 3.5.

As an S/N mask is redetermined after each iteration it changes with the map which can sometimes cause convergence problems. The mask will also depend on the amount of data going into the map and the pixel size. A fixed externally supplied mask can get around these problems.

The sequence below is a summary of the procedure for generating and supplying an external mask. In this example the mask is generated from the map produced by an initial run through the map-maker. Alternatively maps from other observatories can be used.

These steps are followed in the example in Section 7.2.

Step 1 Generate a map covering your region. This may be by simply running the map-maker on your data as shown below.
  % makemap in=’s8*.sdf’ out=850map \
            config=’"^dimmconfig_bright_extended.lis"’

The alternative is to access a map from a different data-set or even a different telescope, e.g. a map downloaded from the Herschel Science Archive. For instructions on converting from FITS to NDF see Appendix G.

Step 2 Make a signal-to-noise map using the Kappa command makesnr.
  % makesnr 850map 850map_snr
Step 3 Threshold this S/N map to set everything below 3σ to 0 and everything above to 1.
  % thresh 850map_snr 850map_mask thrlo=3 newlo=0 thrhi=3 newhi=1

This generates a mask which has an unrealistic hard 3σ cut-off. Step 4 is performed to smooth the the edges of your mask.

Step 4 Smooth the thresholded map with a Gaussian filter of FWHM of 5 pixels (= 20 arcsec). Then it is again thresholded, this time keeping everything above 5 % of the 0 level as the mask and setting the rest to bad.
  % gausmooth 850map_mask 850map_mask_sm fwhm=5
  % thresh 850map_mask_sm 850map_mask_zm thrlo=0.05 newlo=bad \
    thrhi=0.05 newhi=1
Step 5 Finally the map is re-made with this mask supplied as an external file. Notice that the extra parameters required to pick up this external mask are being appended to the configuration file on the command line rather than editing the file itself.
  % makemap in=’s8*.sdf’ out=850map_zm ref=850map_mask_zm \
            config=’"^dimmconfig_bright_extended.lis,ast.zero_mask=1,\
                     ast.zero_snr=0"’

6.7 Skyloop


pict

Figure 6.3: Illustration of the skyloop approach to map-making compared with the standard map-maker.


Traditionally, the map-maker divides a non-contiguous sequence of time series data into chunks. It processes each chuck independently before co-adding them as a final step in the the reduction—see Figure 6.3.

This means for each chunk the map-maker has to start from scratch determining the AST model, and the benefit of long integration times spent building up the signal is lost. Configuration files that use signal-to-noise masks especially suffer from this approach as the signal-to-noise in each individual chunk can remain low and fainter extended structure is not recovered.

The skyloop command is a script that runs makemap multiple times performing just a single iteration on each occasion. It starts by performing a single iteration of makemap from which a co-added map is generated. This map is then supplied as an initial estimate of the sky for the next invocation of makemap. On this next invocation, the initial sky estimate is subtracted from the cleaned time-series data and the COM, GAI, FLT, EXT models are subtracted. This produces a new model of the sky (from the current iteration) to which the sky estimate (from the previous iteration) is then added. In this way the signal from all of the chunks is built up over the iterations and is all included in the final map estimate when convergence is reached.

Be aware that skyloop uses a lot of disk space. Setting environment variable STAR_TEMP to a suitable location before you start will prevent skyloop from crashing if you run out of temporary storage space.

  % setenv STAR_TEMP /some/directory/with/a/lot/of/space

skyloop can then be called in a way very similar to makemap, with a configuration file specified on the command line.

  % skyloop in=^myfiles.lis out=map_skyloop config=^dimmconfig_bright_extended.lis

6.8 Troubleshooting




PROBLEM

POSSIBLE SOLUTION



I have blobs in my map that look like big thumbprints.

Try adding dimmconfig_fix_blobs into your configuration file (see Section 3.8). Note that this sets com.sig_limit = 5 which is somewhat conservative. You can experiment by lowering this, but we would not recommend lower than 2 because you may then lose too many samples.

In addition, check you are not using flt.notfirst = 1 as this can make blobs worse.



I want to recover more extended structure.

There is a trade off between extended emission and noise in your map. If you are willing to accept more low frequency noise you can increase the filter scale with flt.filt_edge_largescale. The default is 480 arcsec but you could try 600 arcsec. To reduce the increased background noise you can set flt.filt_edge_largescale_last = 200. This sets the background filtering to 200 arcsec for the final iteration only, though you can go as low as you want with it. Note, this should be seen as a mainly cosmetic effect as it causes the map to use different filter sizes in different regions, making interpretation of the map difficult.



I want a flatter background.

Try com.perarray = 1, although be aware this will lose structure on scales larger than a sub-array. If you are chasing extended emission see the point above. For a more uniform background set flt.filt_edge_largescale_last to a small value to get harsh filtering on your final iteration.



I have linear striations in my map making my background look scratchy.

Try setting com.corr_abstol = 0.8 [default=0.2]. This rejects more bolometers with deviant common-mode signals. However, as more bolometers are removed there are fewer data available for your final map, resulting in higher noise.



My map will not converge.

Try adding dimmconfig_fix_convergence into your configuration file (see Section 3.8). This prevents the masks from changing after ten iterations.




1These default values are defined in the file $SMURF_DIR/smurf_makemap.def.

2However, the dimmconfig_blank_field configuration over-rides this default.

3The dimmconfig_blank_field configuration omits FLT from modelorder.

415 seconds is half a sub-scan.