MERLIN data: simplified procedure

Overview

This tutorial describes the data reduction for MERLIN 5-GHz radio continuum observations of the radio galaxy 3C277.1. This version covers the simplest steps. For more information see www.merlin.ac.uk. Please note that the general principles are common to most interferometer data but the exact route through data reduction and the parameters such as solution intervals depend on the nature of the array, the observer's intentions and the quality of the data, so different data sets would require modified inputs.

These data were taken using 6 antennas, Ca (32 m) and Mk Ta Da Kn De, all 25 m. Their sensitivities vary, with De being by far the least efficient at this frequency, but is is important for North-South coverage. The array contains baselines up to 217 km. A bandpass correction and approximate flux scale have already been applied and the data are averaged to one 14-MHz channel. All four polarizations (RR RR LR RL) are present but we will just use total intensity. Antenna positions are geocentric, terrestrial E, N positive.

The data can be downloaded from 3C277.1.MULTTB which is a multi-source file containing:

3C277.1 (target)
1300+580 (phase-reference source, known to be point-like)
3C286 (primary flux scale calibrator)
OQ208 (point source flux scale calibrator)

Two flux scale calibrators are used because sources like 3C286 with stable, well-known amplitudes are usually extended over several arcsec or more and hence are heavily resolved by MERLIN. On the other hand, point-like, bright sources like OQ208 are usually variable.

Information: If you want to short-circuit any of the steps, or compare your data (see if you can do better!), 3C277.1.FITS is a multi-source file containing target and calibration sources which has already been (mostly) edited. 3C277.1.SPLIT contains the fully-edited target with all calibration from calibration sources applied, ready for imaging and self-calibration. 3C277.1.CAL.SPLIT is also fully self-calibrated and 3C277.1.FITSIMAGE is the final image in full polarization. These can be loaded into CASA using importuvfits for the first 3 files and importfits for the image.

The script below covers:

Loading and inspecting data
Editing
Initial phase-calibration of calibration sources
Setting the flux scale with reference to 3C286
Time-dependent amplitude calibration of phase-reference source
Applying the solutions to the target, imaging and self-calibration
Image visualisation and measurements

This script is designed to be used by copying or typing the input to one task at a time, based on the text with a cyan background like this.

It attempts to explain the least complicated route through the data calibration and imaging processes. It is possible to combine some steps at the expense of more complicated housekeeping, which might be important for very large data sets. You can plot any calibration step or compare its effects on the data, so you can investigate the impact of changing the values of parameters. Note that the default for most tasks is not to overwrite existing data or calibration files, so the task will fail if the directory or file already exisits. The exception is the output of the task clean, which will carry on further iterations of deconvolution on an existing output image.

Further information on any task e.g. importuvfits is obtained by clicking on the link the first time the task is mentioned, or from the CASA cookbook, or by typing

help('importuvfits')

You can review the possible inputs to a task at any time by e.g.

inp importuvfits

If a parameter is set to False, or if an optional parameter takes several arguments, these may only appear if the parameter is set to True or to a value requiring the arguments. Make sure that you review the inputs to each task using inp before running it. In some places you have to insert the parameters yourself similar to previous inputs. This format is convenient for careful inspection of all parameter values before running each task (with some built-in error checking). Other examples in CASA Guides illustrate formats more suitable for scripting.

Import, inspect and flag the data

This assumes that all data are in the present working directory; if not, use the path in front of the fitsfile name.

Start CASA and set variables for convenience:

prefix ='3C277.1R'
msfile = prefix + '.ms'

The msfile is the Measurement Set (MS) used by CASA for visibility data. If you want to delete all data files created here by CASA at any time, type

!rm -r 3C277.1R*

Load the fits data using importuvfits and list the contents using listobs.

default importuvfits
vis=msfile 
fitsfile='3C277.1.MULTTB'
inp

importuvfits()

listobs()

You will get a few warnings in the Logger which you can ignor, as long as listobs gives something like this in the Logger (some columns/messages omitted):

	##### Begin Task: listobs            #####
================================================
MeasurementSet Name:  /home/amsr/scratch/ALMA/CASANOTES/3C277.1/3C277.1R.ms      MS Version 2
===============================================================================
   Observer:      Project:   
   Observation: MERLIN2
   Data records: 81210       Total integration time = 49913 seconds
   Observed from   18-Apr-1995/17:08:02.0   to   19-Apr-1995/06:59:55.0 (UTC)
   listobs::ms::summary
   ObservationID = 0         ArrayID = 0
  Date        Timerange (UTC)          Scan  FldId FieldName    nVis   Int(s)   SpwIds
  18-Apr-1995/17:08:02.0 - 17:09:54.0     1     10 1300+580     225    7.99     [0]
              17:10:37.0 - 17:17:49.0     2     11 3C277.1      825    7.99     [0]
              17:18:36.0 - 17:19:56.0     3     10 1300+580     165    7.99     [0]
              17:20:35.0 - 17:27:55.0     4     11 3C277.1      840    7.99     [0]
              17:28:42.0 - 17:29:54.0     5     10 1300+580     150    7.99     [0]
.....

              21:43:01.0 - 22:04:53.0    56      0 3C286        2475   7.99     [0]
              22:06:00.0 - 22:59:52.0    57      1 OQ208        6075   7.99     [0]
.....
	           (nVis = Total number of time/baseline visibilities per scan)
..... 
	Fields: 4
	  ID   Code Name         RA            Decl           Epoch   nVis   
	  0         3C286        13:31:08.2873 +30.30.32.9590 J2000   2475   
	  1         OQ208        14:07:00.3944 +28.27.14.6880 J2000   6075   
	  10        1300+580     13:02:52.4649 +57.48.37.6175 J2000   11745  
	  11        3C277.1      12:52:26.2859 +56.34.19.4882 J2000   60915  
	   (nVis = Total number of time/baseline visibilities per field) 
	Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)
	  SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)TotBW(kHz)  Ref(MHz)    Corrs           
	  0          1   TOPO  4993.5       14000        14000     4993.5     RR  RL  LR  LL  
	The SOURCE table is empty: see the FIELD table
	Antennas: 6:
	  ID   Name  Station   Diam.    Long.         Lat.         
	  0    1     DEFFORD   0.0  m   +002.08.35.0  +51.54.50.9  
	  1    2     CAMBRIDG  0.0  m   -000.02.19.5  +51.58.50.2  
	  2    3     KNOCKIN   0.0  m   +002.59.44.9  +52.36.18.4  
	  3    4     DARNHALL  0.0  m   +002.32.03.3  +52.58.18.5  
	  4    5     MK2       0.0  m   +002.18.08.9  +53.02.58.7  
	  5    6     TABLEY    0.0  m   +002.26.38.3  +53.06.16.2  
listobs::::casa
	##### End Task: listobs              #####
	##########################################

Each integration is 8 sec (allowing for the annoying formatting of floating point numbers). This is the duration of a single complex visibility, corresponding to the phase and amplitude for one baseline, one polarization. The polarization products present (also termed correlations or corr in CASA) are LL, LR, RL and RR i.e. full polarization. The target:phase reference source are observed on an 8:2 min cycle. 3C286 and OQ208 are observed in the middle of the run. The diameters of all antennas are given as 0.0 m but as the target does not extend far out in the primary beam we can ignor this. Here, there is only one frequency sub-band, termed a spectral window (spw), containing a single 14-MHz channel.

The paths swept out by the baselines, as seen in projection from the sources, are known as uv tracks. Plot these for the phase-reference source using plotms.

default(plotms)
vis=msfile
xaxis='u'
yaxis='v'
field= '1300+580'
spw='0'
correlation='LL,RR'
inp

plotms

At present the conjugate visibilities are not plotted and the units are metres (not wavelength units) so this may not look like a plot which you are used to, but it shows that the maximum projected baseline is about 210 km and the minimum is 10 km. At a wavelength of 6 cm, this means that the resolution of the array is about 50-60 mas and the largest spatial scale which can be imaged is about 1 arcsec.

Using the Axes tab of plotms to plot amplitude as a function of uvdistance ('UVDist_L' is in wavelength units). You can use the Display tab to colour the points by the first antenna in each baseline. This shows a lot of scatter associated with the antenna contributing baselines around 1.5-2 Mlambda, which is Defford, known to be less sensitive at 6cm wavelength. Colouring by antenna 2 shows that most of the other low points come from antenna 6, Tabley (use the locate facility to find this, easiest to get someone to show you how). This is shown in the figure below.

Right-click on any plot to see it full-size.

The plot shows some bad data but the amplitudes are otherwise fairly constant with uv distance, confirming that the phase-reference, 1300+580, source is point-like.

Experiment with plotting the other sources and confirm which sources are extended or pointlike as expected. In Data Selection set spw 0:3~12 and for Averaging use Time 30 (seconds), which known to be suitable to maximise signal-to-noise in raw data, as demonstrated later.

The target, 3C277.1, has more flux on short baselines and complex structure.

Change the axes to plot the phase reference source amplitudes as a function of time. Zoom in on the later times which show bad data. Use the Tools tab to activate Hover and Display. This means that if you select phase-reference data to flag, you can note the same time-range to flag for the target. Use Locate to see which antennas are responsible (the times etc. to flag are given below, so if in a hurry don't spend too much time on this).

Phase-ref amplitudes (averaged every 30 sec). The Locate function shows which antennas are responsible for each of the major patches of bad data. You can flag these interactively or by noting the time ranges and using the flagdata task, see next step. You can also flag the isolated low points (which are unlikely also to affect the target).

Use flagdata to edit the main periods of bad data on the phase-reference and target together. Small falls in amplitude will be taken care of by calibration, but multiplying up very faint data adds too much noise.

default flagdata
vis=msfile
spw='0'
field='1300+580','3C277.1'
antenna='6&*'
timerange='1995/04/19/02:01:00~1995/04/19/03:49:00'

inp
flagdata

antenna='4&*'
timerange='1995/04/19/05:23:00~1995/04/19/05:38:00'
flagdata

antenna='2&*'
timerange='1995/04/19/06:24:00~1995/04/19/06:35:00'
flagdata

Check the messages in the logger and check by re-plotting the data. When visibility data are loaded into CASA, as well as creating a MS, a .flagversions directory is created with a (usually empty) version 'Original'. You can use flagmanager to back up flags (or restore a previous version). If all the flagging so far has worked, back up the flags.

default(flagmanager)
vis=msfile
mode='save'
versionname = msfile+'_calsource.flags'

inp
flagmanager()

Inspect the amplitudes of OQ208 and 3C286 as a function of time and then phases as a function of time. For point sources, with a good position, it should be flat at zero. In the plots below, the phase wanders about but can be traced easily from integration to integration. If the phase looks like pure noise, you could try averaging up but if that still does no good on a calibration source, the time should be flagged. Here, it looks all OK unless you plot the data already flagged due to bad amplitude.

The phases for 1300+580, for all baselines to the most distant antenna (Cambridge), change by up to a full turn in a day, seen in the upper plot. The middle plot shows that the phase only changes by about 10 deg or less in 5 min. The data shown in red in the lowest plot was already flagged due to low amplitudes and also has bad (apparently random) phase.

The phase-reference source is observed for about 90 sec in between each scan on the target. The plots above show that if we choose a 30-sec solution interval, we will achieve a phase accuracy of about 5 degrees or better, and this also allows accurate interpolation onto the target data.

First phase solutions

The origin of phase is arbitrary and in order to achieve consistent calibration, especially for polarization, the phase of one antenna is set to zero and the other phases are worked out relative to that. Plot the antenna distribution using plotants to chose a reference antenna with good data with many short baselines (slower phase-rates so less chance of solutions failing).

vis = msfile plotants

Ignor messages about not finding the position of MERLIN. In some versions of CASA the plot may come out with E and W reversed. I suggest '5' (Mk 2) as the refant.

MERLIN antenna positions

Close the plotter or you may have file lock problems in future. If you do find that plotms or some other task is hanging, or you get messages about waiting for a lock, close any (other) CASA visualisation tools and type

clearstat

The task gaincal is used to derive time-dependent solutions. The default model, of a point source at the phase centre, will be used to use a chi-squared minimisation method to derive the corrections needed per antenna to make the phases as close as possible to zero. This will average all unflagged channels.

default('gaincal')
vis=msfile
caltable=prefix+'_cals.phcal'
field = '3C286,OQ208,1300+580'   
spw='0'
solint = '0.5min'
refant='5'    
calmode='p'
minsnr = 3.0
minblperant= 3

inp
gaincal()

Check the logger and plot the solutions:

default('plotcal')
caltable=prefix+'_cals.phcal'
field = '1300+580'
subplot = 611
iteration='antenna'
xaxis='time'
yaxis='phase'
plotrange=[-1, -1, -180, 180]

inp
plotcal()

Phase reference source phase solutions (the range is [-180,+180] deg)

Check the solutions for the other calibration sources. Even where the phase solutions are going round quite fast, they should not be random. At present in CASA, the only way to check the effect of the solutions on the data is to apply the solution table using href="http://casa.nrao.edu/docs/taskref/applycal.html">applycal (eventually it will be possible to apply tables on the fly whilst running a task). As the MERLIN antennas have different sensitivities we use calwt=False (so that data for the less sensitive antennas are not given too small a weight to be included properly in self-calibration).

default('applycal')
vis = msfile
gaintable= prefix+'_cals.phcal'
field='1300+580,3C286,OQ208'
calwt=False   
parang=False

inp
applycal()

Plot the corrected phase (you can only select this after displaying the data).

default('plotms')
vis=msfile
xaxis='time'
yaxis='phase'
field='1300+580'
spw='0'
correlation='LL,RR'
avgtime='30'

inp
plotms

In the viewer, Axes tab, set the Data Column to 'corrected'. You should see that the phase is within about 5 degrees of zero except for a short period near the middle of the run, as shown below.

Corrected phase for the phase-reference source (all antennas, 30-sec average).

Setting the flux scale and time-dependent amplitude calibration

3C286 is a standard flux and polarization calibration source, with properties tabulated by Baars (1977). It is somewhat resolved even on the shortest MERLIN baselines at 5 GHz, hence a sub-standard flux density (96% of Baars) is inserted in the MODEL data column of the measurement set using setjy.

default('setjy')
vis = msfile
field = '3C286'
fluxdensity=[7.077,0.0]
inp

setjy()

gaincal is now used to derive amplitude solutions using a point model for the primary flux calibration source, 3C286 and the point source OQ208. The plot below shows that 3C286 is heavily resolved and we will select a uv range less than 500,000 wavelengths. Use plotms to reproduce this and check that the data in this range does come from at least 3 antennas. This gives 3 baselines and is the minimum possible for normal phase solutions.

3C286 is heavily resolved on all but the shortest baselines.

The phase and bandpass solution tables must be applied since gaincal uses the Data column. Set the uvrange to the range where 3C286 is almost unresolved.

default gaincal
vis=msfile
caltable=prefix+'_cals.a1cal'
field='3C286,OQ208'
selectdata=True
uvrange='0~500klambda'
solint='2min'
minblperant =2
refant='5'
gaintable=prefix+'_cals.phcal'
calmode='ap'

inp
gaincal()

Check the solutions. The solutions for 3C286 are close to unity but those for OQ208 are not because we do not yet know its correct flux density. Note that even asking for 4 plots only produces 3, confirming that solutions were only found for the 3 antennas providing the shortest baselines.

caltable = prefix+'_cals.a1cal'
xaxis='time'
yaxis='amp'
field=''
antenna = '4,5,6'
plotrange=[-1,-1,-1,-1]
subplot=411
iteration='antenna'
plotcal()

fluxscale is used to compare the solutions for OQ208 with those of 3C286 to find the flux density of the former.

default('fluxscale')
vis=msfile
caltable=prefix+'_cals.a1cal'
fluxtable=prefix+'_OQ208.flux'
reference='3C286'
transfer='OQ208'

inp
fluxscale()

Note the OQ208 flux density in the logger; the scatter should be less than 1%. The exact value will depend on whether you have done much flagging. I got

Flux density for OQ208 in SpW=0 is: 2.43903 +/- 0.00453 (SNR = 537, nAnt= 3)

The fluxtable will contain the proper scaling, but since we have only used a subset of antennas we have to set the flux density of OQ208 (insert your value below) and repeat the amplitude calibration using all antennas:

default('setjy')
vis = msfile
field='OQ208'
fluxdensity=[2.439,0] 

inp
setjy()

default gaincal
vis=msfile
selectdata=False
field='OQ208,1300+580'
solint = '2min'
uvrange=''
antennas = ''
refant = '5'
minsnr=3
minblperant = 4
caltable = prefix+'_cals.a2cal'
gaintable = prefix+'_cals.phcal'

inp
gaincal()

Use plotcal to check the solutions (see previously); assuming that they are OK use fluxscale to find the flux density of the phase reference source on all baselines

vis=msfile
caltable=prefix+'_cals.a2cal'
fluxtable=prefix+'_1300+580.flux'
reference='OQ208'
transfer='1300+580'

inp
fluxscale()

Note the flux density of 1300+580 from the logger, check that it is close to the expected value. I got

Flux density for 1300+580 in SpW=0 is: 0.433907 +/- 0.000773 (SNR = 561, nAnt= 6)

3C277.1R_1300+580.flux now contains correctly-scaled amplitude solutions (check in plotcal if you have time).

Apply the calibration and prepare for imaging

We now have solution tables for phase and amplitude as a function of time, and as a function of channel, and for polarization leakage and polarization angle. Apply these to the target and the phase-reference source using applycal.

default('applycal')
vis = msfile
gaintable= prefix+'_cals.phcal',default('applycal')
vis = msfile
gaintable= prefix + '_cals.phcal', prefix+'_1300+580.flux',
field='1300+580,3C277.1'
calwt=False   
parang=True  

inp
applycal()

This fills the CORRECTED data column with visibilities modified by the various calibration factors. Use plotms to check the data (and compare it with the uncalibrated data) split out the target. Further imaging and self-calibration are quicker and simpler using a single-source file.

default('split')
vis=msfile
datacolumn = 'corrected'
outputvis='3C277.1R.split.ms'
field='3C277.1'

inp
split()

Inspect the calibrated 3C277.1 data using plotxy, which is supposed to be replaced by plotms but the latter does not yet allow multiple sub-plots. In plotxy, use the Next button to see all baselines. Any obviously isolated points can be flagged. Inspect phase also and use the zoom if necessary but remember some errors will be corrected by self-calibration.. After later self-calibration, you might have another look and edit remaining isolated points.

default(plotxy)
vis = '3C277.1R.split.ms'
xaxis='time'
yaxis='amp'
datacolumn='data'
correlation='RR LL'
iteration='baseline'
field = '3C277.1'
subplot=511

plotxy()

Amplitude and phase of the target after applying reference source and bandpass calibration, on all baselines to the reference antenna. The signal-to-noise is good enough to see the source in short time intervals on all baselines so it is possible to use self-calibration.

Self calibration and imaging of the target

Imaging consists of two processes. The calibrated visibility data are Fourier transformed to make a Dirty Map. This is what will appear first in the Viewer during interactive cleaning. It contains many artefacts due to the sparse sampling of the sky by the interferometer. A plot of the Dirty Beam is shown lower down, demonstrating the sidelobes above 20% of the peak which extend over several beamwidths (defined as FWHM, about 50 mas). Cleaning is a process in which the brightest point(s) in the map are identified and a fraction, say 10%, of their flux density is recorded as Clean Components (CC) and the sidelobe response for the dirty beam centred on that position is scaled and subtracted from the map. After a certain number of iterations (known as minor cycles) the accumulated CCs are FT'd and subtracted from the visibility data, the remaining visibility data are FT'd to the map plane and the process continues (major cycles). When there is nothing left much above the noise, the cumulative list of CC are added to the residuals and a final Clean map is made, convolved with the restoring beam.

Image the target using clean. The longest baseline is 217 km and the wavelength is 6 cm. Use this to estimate the interferometer resolution and hence a suitable pixel size to give at least 3 pixels across the beam. A minimum image size of at least 512 pixels is needed for for adequate deconvolution of the rather extended dirty beam produced by the sparse MERLIN array. psfmode='hogbom' provides a large beam patch for deconvolving extended sidelobes.

Note that if interactive=True you have to set a cleaning mask yourself. This approach is suitable for a bright, extended source such as 3C277.1. Start clean as described below. When the Viewer appears, experiment with these steps:

Zoom in so that you get a good view of the two brightest patches of emission.
The flux density under the cursor is reported below the image, identify approximately what colour represents zero and use the Polygon option to draw two roughly ellptical regions just above zero round the two bright patches.
When you are happy with one region, close it either by drawing a line which is continuous or by double-clicking until green edit points appear. Double-click again inside the green shape until it goes white-ish in order to fix it and move on to the next region, if wanted.
If you draw overlapping shapes, they will merge. If you want to remove all or part of a shape, use the Erase button, and draw the area to be removed (don't forget to return to Add mode afterwards).
To start cleaning, click on the green curved arrow and watch the messages in the Logger.
At first, when the Viewer reappears after each 100 iterations (or whatever number you have set), you will see that, as the brighter emission is Cleaned away, a larger area appears including the quasar-like core; extend/add masks to include the newly-visible, fainter emission.
Eventually, the whole image will appear fairly smooth and the target regions are no longer brighter; click the red button to stop cleaning.

default('clean')
vis='3C277.1R.split.ms'
cell=['0.015arcsec', '0.015arcsec']
imsize=[512,512]
psfmode='hogbom'    
imagermode='csclean'
cyclefactor=4
stokes='I'
interactive=True
niter=1000
gain=0.05
field='3C277.1'
imagename='3C277.1R_clean1'
interactive=True
stokes='I'
mask=''

inp
clean()

Don't forget to double-click inside each mask you draw to make it turn white(ish) to activate it. Ignor warnings about not finding position of MERLIN. To save the clean inputs for future, use saveinputs

saveinputs ('clean', '3C277.1.clean1.saved')

Check progress in the logger. In this instance we have let CASA fit the natural synthesised beam, the size is reported in the logger (and also recorded in the image header). The Clean model is left in the MS in the MODEL_DATA column. Use the viewer to load 3C277.1R_clean1.image as a raster map and use the 'data display options' to tweak it. You can access these via the spanner icon or the Data menu option Ajust.

viewer('3C277.1R_clean1.image')

You can also plot the Dirty Beam, '3C277.1R_clean1.psf'.

The dirty beam and its sidelobes.

There are several ways to measure the noise and the peak. If you draw a box and double-click, the rms etc. will be displayed. To set the measurement in a variable (here called rms_target) use imstat, setting the limits to avoid the image.

default('imstat')
imagename='3C277.1R_clean1.image'
box='20,20,492,120'
stokes='I'
noise_target = imstat()
rms_target=noise_target['rms'][0]
print rms_target

I got 0.00087 Jy. Use the menu Data>Load, select 3C277.1R_clean1.image as a contour image and use the 'data display options' to set the base contour level to the rms_noise and the relative contour levels to [-3,3,6,12,24,48,96].

First image of the target, contours at [-0.0025,0.0025, 0.005, 0.01, 0.02, 0.04, 0.08] Jy/beam

You can also overlay the Clean Components (3C277.1R_clean1.model) as a Marker Map using the Data > Open menu (set the x and y increments to 1.) Self-calibrate with gaincal using the model (the Fourier transform of the clean components) left in the MS. It is OK to selfcal in I if V is zero, which is usually the case at all but the highest resolution for extragalactic cm-wave radio continuum.

default('gaincal')
vis='3C277.1R.split.ms'
caltable=prefix+'_self1.phcal'
field='3C277.1'
refant='5'
minblperant=3
minsnr=1
calmode='p'
solint = '2min'
parang=True

inp
gaincal()

Check the logger and inspect the solutions

caltable=prefix+'_self1.phcal'
xaxis='time'
yaxis='phase'
subplot=611
iteration='antenna'
plotrange = [-1, -1, -180, 180]

plotcal()

You expect to see just residual solutions; if they look OK use applycal.

Target phase solutions (the range is [-180,+180] deg).

default('applycal')
vis='3C277.1R.split.ms'
gaintable=prefix+'_self1.phcal'
field='3C277.1'
parang=True 
calwt=False

inp
applycal()

In clean, set the threshold to the measured noise, rms_target Restore the saved inputs (3C277.1.clean1.saved is just a text file, you can view or edit it) using the python command execfile. If you saved your boxes give the file name as the argument of mask. You can change these interactively in the clean viewer. Restore the previous, saved clean inputs.

execfile('3C277.1.clean1.saved')

Inspect and adapt the inputs, using the previous noise level as the minimum clean flux. You should stop cleaning when the minimum clean components have sidelobes fainter than the noise. You can calculate the stopping criteria more accurately by inspecting the sidelobe ratio of the dirty beam and comparing that with the expected final noise level. The mask is set below to use the final mask from the previous round of clean. You can leave this blank to start afresh, or you can modify the input mask as you did before.

inp

imagename = '3C277.1R_clean2'
threshold=str(rms_target)+'Jy'
niter=2000
mask='3C277.1R_clean1.mask'

clean()

saveinputs ('clean', '3C277.1.clean2.saved')

Inspect 3C277.1R_clean2.image as previously in the viewer and check that the noise has gone down.

default('imstat')
imagename='3C277.1R_clean2.image'
box='20,20,492,120'
stokes='I'
noise_target = imstat()
rms_target=noise_target['rms'][0]
print rms_target

In this example we go straight to amplitude self-calibration. If the previous phase solutions showed large changes you might do more phase-only calibration first to improve the model, perhaps trying shorter solution intervals and different boxes. Gaincal only works from the original, not corrected data so apply the previous phase self-cal solutions as the gaintable in gaincal.

default('gaincal')
vis='3C277.1R.split.ms'
caltable=prefix+'_self2.ampcal'
field='3C277.1'
refant='5'
minsnr=1
calmode = 'ap'
minblperant=4
solint='3min'
gaintable=prefix+'_self1.phcal'
parang=True

inp
gaincal()

Inspect the solutions.

caltable=prefix+'_self2.ampcal'
iteration='antenna'
subplot=611
yaxis='amp'
plotrange = [-1, -1, 0.8, 1.2]
plotcal()

The amplitude corrections should be within a few percent of unity. Plot the phase corrections as well, these should be very small residuals.

Target amplitude corrections, range [0.8,1.2]

Apply both 'p' and 'ap' calibration tables with applycal

vis='3C277.1R.split.ms'
gaintable=prefix+'_self1.phcal',prefix+'_self2.ampcal'
calwt=False
parang=True

inp
applycal()

Finally, clean. Remember that you can change the threshold or the number of iterations/cycles, or modify the mask, interactively.

execfile('3C277.1.clean2.saved')
inp

imagename = '3C277.1R_clean3'
mask='3C277.1R_clean2.mask'
threshold=str(rms_target)+'Jy'
niter=5000
clean()

View and check that the noise has gone down:

default('imstat')
imagename='3C277.1R_clean3.image'
box='20,20,492,120'
stokes='I'
noise_target = imstat()
rms_target=noise_target['rms'][0]
print rms_target

You can check the accuracy of your model by plotting the Fourier transform of the Clean components (in the MODEL data column of the MS) against the CORRECTED visibilities:

default('plotxy')
vis='3C277.1R.split.ms'
xaxis='uvdist'
yaxis='amp'
datacolumn='corrected'
selectdata=True
correlation='LL'
field='3C277.1'
timebin='180' # 3 min, same as solution interval
plotcolor='blue'
overplot=true
plotxy()

datacolumn='model'
plotcolor='red'
plotxy()

Target model (red) overplotted on the calibrated visibility data.

Some scatter in the data is due to real structure below the CC cutoff and to unavoidable noise, which is worse on intermediate baselines because these include Defford, the least sensitive antenna at 5 GHz. If the model is a very poor fit, try more editing and/or self-calibration and/or change the clean inputs

3C277.1 after self-calibration, contour levels in logarithmic multiples of 0.4 mJy/beam (3 sigma_rms).

If time, try to improve the image... read the help for clean and try multi-scale clean, experiment with a smaller gain factor or more careful boxing. You could also split out 3C286 and 1300+580 and image these in full polarization as above. 3C286 will have strong sidelobes due to the poor uv coverage but boxing the central peak and just cleaning for about 100 iterations should produce an image good enough to check the polarization angle, which should be 33 deg. Additional techniques in cleaning include:

You can change the number of iterations per cycle, total number of cycles and the threshold interactively. If you click the blue button, the Clean will proceed without further interaction till it reaches the threshold or the maximum number of iterations.
You can modify the colour table or even add a wedge using the 'spanner' tool.
You can open another image as a contour map, e.g. from a previous round of cleaning, and use e.g. the 3-sigma rms contour to guide you in drawing the mask.
The Tools menu lets you load or save a mask in various ways (as an alternative to recycling the mask from the previous cleaning). See the Cookbook for details.

Anita Richards 19 June 2010 revised January 2011