Self-Calibration Template CASA 6.1.1

This guide continues with products created in Image the Continuum Template. All imaging parameters are set in the previous guide. You will need calibrated_final_cont.ms to proceed. Commands for this guide can be found in scriptForImaging_template.py available on github.

Self-Calibration on the Continuum (optional)

If you have a high signal-to-noise detection of your source, you can use the source itself to calibrate the phases and potentially the amplitudes of the visibilities as a function of time. This technique is called self-calibration and takes advantage of the fact that interferometer data is over-constrained (we have more equations than we have solutions). The recommended signal to noise depends on the number of antennas. For an array of 40 antennas, you should have a Peak/RMS of at least 45 in your continuum image to attempt self-calibration. You can use the following relationship to find the minimum Peak/RMS for your dataset.

[math]\displaystyle{ \frac{Peak}{RMS} \gt 3 \sqrt{N-3} \sqrt{\frac{t_{int}}{t_{solint}}} }[/math]

where:

Peak is the maximum value in the continuum image.

RMS is the RMS of the non-self-calibrated continuum image.

N is the number of antennas.

t_int is the total on source integration time.

t_solint is the solution interval used to calculate the corrections. In this case, use the scan length.

See the NRAO Live! presentation When, why, and how to do self-calibration for more information on this relationship.

Self-calibration is an iterative process. You start by generating a model of your source using tclean. Then you use this model to determine the gains as a function of time using the gaincal task. Finally you apply these solutions to the data and re-image. These steps are repeated until you are happy with your model or the solution interval is too short to reach the necessary signal-to-noise. In general, you start with phase-only self-calibration and only do amplitude calibration at the end of the self-calibration process if there are amplitude-based gain artifacts in the data (see the Imaging and Deconvolution talk from the Synthesis Imaging Summer School for more detail. Amplitude calibration should be used with caution because it has the potential to change the fluxes of sources in your data.

Self-calibration can either be performed with line or continuum data. Here we demonstrate how to do it with continuum data. The principle is the same for line data, except you form your image using the brightest line emission. For an example of self-calibration with line data, see the VLA high frequency spectral line tutorial for IRC+10216.

---

The first thing you should do when beginning self-calibration is to save the original flags in your data set. This step is necessary because when you apply the calibration to a data set it flags data that is not associated with a solution. Applying a bad calibration table to your data can result in the excessive flagging of your data. Saving your original flags allows you to restore the original state of your data easily. We also recommend you save your flags after each applycal command. This will allow you to go back to a given point in the self-cal process without starting from the beginning. See Self_Calibration_Template#Restart_Self-Calibration if you need to start over or go back a step in self-cal.

# in CASA
flagmanager(vis=contvis,mode='save',versionname='before_selfcal',merge='replace')

Now we can set up the parameters for the self calibration. Begin by setting the parameters for your ms and continuum image name. The ms calibrated_final_cont.ms was created in the previous part of this guide: Image the Continuum Template.

# in CASA
contvis = 'calibrated_final_cont.ms'         
contimagename = 'calibrated_final_cont'

If you run the self-calibration process more than once, you will need to remove any models that could be left in the measurement set. This can can be accomplished by using the delmod task as well as the clearcal task (which will need to be uncommented in the code below).

# in CASA
# delmod(vis=contvis,otf=True,scr=True)
# clearcal(vis=contvis)

You should use the same reference antenna that was used during calibration. For pipeline reductions, this can be found in the stage hif_refant as the first antenna in the list for each execution. For manual reductions, the reference antenna will be stated near the top of *.ms.scriptForCalibration.py. If there are multiple executions, make sure to use an antenna that has good solutions and is present in all executions. The Analysis Utilities task au.commonAntennas can help you find antennas that meet this criteria. If you haven't installed Analysis Utilities, see Obtaining Analysis Utilities for instructions. Tasks such as plotants or listobs/vishead can also be used to find antennas in all executions.

# in CASA
refant = 'DV09' # pick a reference antenna.

In the example below, we combine the signal from all spectral windows to improve the signal-to-noise for our gain solution. When combining the spectral windows, you need to map the solution from the combined spectral window (0) to the individual spectral windows using the spwmap parameter. This parameter is a list where the index of the element in the list indicates the spectral window and the value for that index the window that it is mapped to. For example, if we have three spectral windows in the original data set and use combine='spw' for our gain solution, we set spwmap=[0,0,0] to map spw=0 to spw=0, spw=1 to spw=0, and spw=2 to spw=0.

# in CASA
spwmap = [0,0,0] # mapping self-calibration solutions to individual spectral windows. Generally an array of n zeroes, where n is the number of spectral windows in the data sets.

You then begin the process of shallowly cleaning your continuum data to create an initial model for your data in the model column of your data set. Usually you only should do at most a few hundred iterations on the brightest source(s) in the field.

# in CASA
# shallow clean on the continuum

for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
    rmtables(contimagename + '_p0'+ ext)

tclean(vis=contvis,
       imagename=contimagename + '_p0',
       field=field,
       #phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight=True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize=imsize, 
       cell= cell, 
       weighting=weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       savemodel='modelcolumn',
       usepointing=False)

# Note number of iterations performed.

Before proceeding, make sure that tclean has saved the model column. You should see a message like the following

>>> INFO .... ------ Predict Model ------
>>> INFO ... Saving model column

If you don't see this message, you should repeat the above tclean command with the following modifications. This will populate the model column.

niter=0
calcpsf=False
calcres=False

As you continue with the self-cal, always make sure the model column is saved after running tclean.

Next you take that model and use it to determine the per-antenna phase solutions as a function of time using gaincal. We use the rmtables command here to completely eradicate any previous solution table from CASA's memory. Note that we start here with a fairly long solution interval that is the length of a scan (solint='inf').

# in CASA
# per scan solution
rmtables('pcal1')
gaincal(vis=contvis,
        caltable='pcal1',
        field=field,
        gaintype='T',
        refant=refant,
        calmode='p',
        combine='spw',
        solint='inf',
        minsnr=3.0,
        minblperant=6)

Inspect the logging messages output by gaincal to see how many solutions were expected/attempted/succeeded. If you have a large number of failed solutions, do not proceed further! This usually means that you do not have enough signal-to-noise on your source to proceed with self-calibration. Note that if you apply a calibration table with many failed solutions to the data, it will flag the data associated with these solutions. You may try changing the minsnr and minblperant parameters in the gaincal call above to find more solutions, but this option should be used with caution because these parameters are already low. This is especially relevant for 7-meter datasets, where achieving 6 baselines per antenna may be difficult. If you choose to modify these parameters, make sure to examine the results in detail.

Calibration solve statistics per spw: (expected/attempted/succeeded):
   Spw 0: 235/235/235
   Spw 1: 235/235/235
   Spw 2: 235/235/235

Figure 1: Solutions from the first round of self-cal.

You should also check the solutions (above) that were obtained using plotcal. The solutions should vary smoothly with time and there should not be any large outliers. The reference antenna’s phases should be flat in time.

# in CASA
# Check the solution
plotms(caltable='pcal1',
        xaxis='time',
        yaxis='phase',
        iteration='antenna',
        plotrange=[0,0,-180,180])

If you are satisfied with the solutions, apply them to the ms. Note: this step usually takes a while, so go ahead and get a cup of coffee.

Figure 2: Continuum image after the first round of self-cal.

# in CASA
# apply the calibration to the data for next round of imaging
applycal(vis=contvis,
         field=field,
         spwmap=spwmap,
         gaintable=['pcal1'],
         gainfield='',
         calwt=False,
         flagbackup=False,
         interp='linearperobs')

Save the flags in case you need to go back to this step.

# in CASA
flagmanager(vis=contvis,mode='save',versionname='after_pcal1')

Using our new gain solutions, we can generate an improved model for our source using tclean. During this second tclean iteration, you should clean a bit deeper than previously, but not all the way down to the noise. Remember the goal here is to build up a good model for your source, so you don't want to include things that are not real emission.

# in CASA
# clean deeper
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
    rmtables(contimagename + '_p1'+ ext)

tclean(vis=contvis,
       imagename=contimagename + '_p1',
       field=field,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight=True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize=imsize, 
       cell=cell, 
       weighting=weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       #pbcor=True, #if final image
       savemodel='modelcolumn',
       usepointing=False)
# Note number of iterations performed.

You should inspect the new image in the viewer to assess how well self-calibration worked. In general, the noise should be lower and the signal-to-noise ratio higher than in the original image and there should be no major new artifacts in the image. If the signal-to-noise ratio does not increase, you may not have enough signal-to-noise on your target to run self-calibration. If you don't see an improvement in the image, you may wish to try a new reference antenna by first clearing the calibration with clearcal (see commands above). If this still does not help improve the selfcal images you may proceed to continuum subtraction (if you are working on a line project).

If the image has improved you can proceed with the self-cal and generate a new gain solution table. The next gaincal call uses a shorter solution interval (solint) to generate solutions on a shorter time interval. Note that we always want to generate our model based on the latest gaincal solutions. The calibration tables, however, should be generated by comparing the original visibilities to the model. This prevents bad solutions from propagating through different iterations of self-calibration.

# in CASA
# shorter solution
rmtables('pcal2')
gaincal(vis=contvis,
        field=field,
        caltable='pcal2',
        gaintype='T',
        refant=refant,
        calmode='p',
        combine='spw',
        solint='30.25s', # solint=30.25s gets you five 12m integrations, while solint=50.5s gets you five 7m integration
        minsnr=3.0,
        minblperant=6)

# Check the solution
plotms(caltable='pcal2',
        xaxis='time',
        yaxis='phase',
        iteration='antenna',
        plotrange=[0,0,-180,180])

# apply the calibration to the data for next round of imaging
applycal(vis=contvis,
         spwmap=spwmap,
         field=field,
         gaintable=['pcal2'],
         gainfield='',
         calwt=False,
         flagbackup=False,
         interp='linearperobs')
flagmanager(vis=contvis,mode='save',versionname='after_pcal2')

Now we can repeat the tclean/gaincal/applycal process until we start seeing many failed solutions and only small improvements to the data. Remember to clean a bit deeper each time to improve your model and be cautious about what emission you are including in your tclean mask. You should also decrease the solution interval each time to better model the gain variations with time. Note that it generally only take 3-4 rounds of phase self-calibration to produce a good solution.

# in CASA
# clean deeper
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
    rmtables(contimagename + '_p2'+ ext)

tclean(vis=contvis,
       imagename=contimagename + '_p2',
       field=field,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight=True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize=imsize, 
       cell=cell, 
       weighting=weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       #pbcor=True, #if final image
       savemodel='modelcolumn',
       usepointing=False)
# Note number of iterations performed.

# shorter solution
rmtables('pcal3')
gaincal(vis=contvis,
        field=field,
        caltable='pcal3',
        gaintype='T',
        refant=refant,
        calmode='p',
        combine='spw',
        solint='int',
        minsnr=3.0,
        minblperant=6)

# Check the solution
plotms(caltable='pcal3',
        xaxis='time',
        yaxis='phase',
        iteration='antenna',
        plotrange=[0,0,-180,180])

# apply the calibration to the data for next round of imaging
applycal(vis=contvis,
         spwmap=spwmap,
         field=field,
         gaintable=['pcal3'],
         gainfield='',
         calwt=False,
         flagbackup=False,
         interp='linearperobs')

flagmanager(vis=contvis,mode='save',versionname='after_pcal3')

for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
    rmtables(contimagename + '_p3'+ ext)

tclean(vis=contvis,
       imagename=contimagename + '_p3',
       field=field,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight=True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize=imsize, 
       cell=cell, 
       weighting=weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       #pbcor=True, #if final image
       savemodel='modelcolumn',
       usepointing=False)

# Note number of iterations performed.

Figure 3: Amplitude based residuals after phase self-cal.

Figure 4: Image after amplitude self-cal.

The next section performs an amplitude calibration. This is an optional part of self-calibration. If you are happy with the results of your phase calibration, you can stop here. However, if you see amplitude-based artifacts, you can attempt to improve the situation using amplitude self-calibration. Since amplitude solutions are inherently less constrained than phase solutions, we use a longer solution interval only here. If you have a complex source with lots of extended emission, you may set a uvrange limit on the data to avoid downweighting the large scale emission

While the phase calibration involved iterating with different solution intervals, the amplitude self-calibration will only use an infinite solution interval.

# in CASA
rmtables('apcal')
gaincal(vis=contvis,
        field=field,
        caltable='apcal',
        gaintype='T',
        refant=refant,
        calmode='ap',
        combine='spw',
        solint='inf',
        minsnr=3.0,
        minblperant=6,
#        uvrange='>50m', # may need to use to exclude extended emission
        gaintable='pcal3',
        spwmap=spwmap,
        solnorm=True)

plotms(caltable='apcal',
        xaxis='time',
        yaxis='amp',
        iteration='antenna',
        plotrange=[0,0,0.2,1.8])

applycal(vis=contvis,
         spwmap=[spwmap,spwmap], # select which spws to apply the solutions for each table
         field=field,
         gaintable=['pcal3','apcal'],
         gainfield='',
         calwt=False,
         flagbackup=False,
         interp='linearperobs')

flagmanager(vis=contvis,mode='save',versionname='after_apcal')

# Make amplitude and phase self-calibrated image.
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
    rmtables(contimagename + '_ap'+ ext)


tclean(vis=contvis,
       imagename=contimagename + '_ap',
       field=field,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight=True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize=imsize, 
       cell=cell, 
       weighting=weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       savemodel='modelcolumn',
       pbcor=True, # apply the primary beam correction since this is the last image.
       usepointing=False)

Compare the final image to the initial image and see if the image dynamic range (ratio betwen peak flux and rms noise) has improved and phase- and amplitude-based errors have improved.

Finally, you should split out the results of your self-calibration for safe-keeping.

# in CASA
split(vis=contvis,
outputvis=contvis+'.selfcal', datacolumn='corrected')

Restart Self-Calibration

If you would like to revert to a certain point in the self-cal process, use the following commands to do so.

# in CASA
# uncomment the following to revert to a given point in the iterative process (here after pcal2 has been applied)
# flagmanager(vis=contvis, mode='restore',versionname='after_pcal2')
# clearcal(contvis)
# delmod(contvis,field=field,otf=True)

Then return to the applycal step that applied the desired calibration table. In the example code above, you would return to where the table pcal2 was applied.

If you are unhappy with the self-calibration, use the clearcal and delmod tasks and restore the original flags to return your ms to it’s original pre-self-cal state.

# in CASA
# uncomment the following to revert to pre self-cal ms
# flagmanager(vis=contvis, mode='restore',versionname='before_selfcal')
# clearcal(contvis)
# delmod(contvis,field=field,otf=True)

If you have line data, you will subtract the continuum and apply the self-calibration results to the line data on the next page, Spectral Line Imaging Template.

Return to the Main Page