View source for Self-Calibration Template CASA 6.6.1

This guide continues with products created in '''[[Image_Continuum | 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 [https://github.com/aakepley/ALMAImagingScript github].

This guide features CARTA, the “Cube Analysis and Rendering Tool for Astronomy,” which is the new NRAO visualization tool for images and cubes. The CASA viewer (imview) has not been maintained for a few years and will be removed from future versions of CASA. We strongly recommend using CARTA, as it provides a much more efficient, stable, and feature rich user experience. A comparison of the CASA viewer and CARTA, as well as instructions on how to use CARTA at NRAO, is provided in the CARTA section of the CASA docs. 

== 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>
\frac{Peak}{RMS} > 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<sub>int</sub>'' is the total on source integration time.
 
''t<sub>solint</sub>'' is the solution interval used to calculate the corrections. In this case, use the scan length.

See the NRAO Live! presentation [https://science.nrao.edu/facilities/alma/naasc-workshops/nrao-cd-wm16/Selfcal_Madison.pdf 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_6.6.1}}. Then you use this model to determine the gains as a function of time using the {{gaincal_6.6.1}} 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 [https://science.nrao.edu/science/meetings/2016/15th-synthesis-imaging-workshop/documents/wilner_vla16.pdf 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 [https://casaguides.nrao.edu/index.php?title=VLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216#UV_Continuum_Subtraction_and_Setting_Up_for_Self-Calibration 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_6.6.1}} 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.

<source lang="python">
# in CASA
flagmanager(vis=contvis,mode='save',versionname='before_selfcal',merge='replace')
</source>

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_Continuum | Image the Continuum Template]]'''.

<source lang="python">
# in CASA
contvis = 'calibrated_final_cont.ms'         
contimagename = 'calibrated_final_cont'
</source>

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_6.6.1}} task as well as the {{clearcal_6.6.1}} task (which will need to be uncommented in the code below).

<source lang="python">
# in CASA
# delmod(vis=contvis,otf=True,scr=True)
# clearcal(vis=contvis)
</source>

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 [https://safe.nrao.edu/wiki/bin/view/ALMA/CommonAntennas au.commonAntennas] can help you find antennas that meet this criteria. If you haven't installed Analysis Utilities, see [https://casaguides.nrao.edu/index.php?title=Analysis_Utilities Obtaining Analysis Utilities] for instructions. Tasks such as {{plotants_6.6.1}} or {{listobs_6.6.1}}/{{vishead_6.6.1}} can also be used to find antennas in all executions. 

<source lang="python">
# in CASA
refant = 'DV09' # pick a reference antenna.
</source>

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.

<source lang="python">
# 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.
</source>

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.


<source lang="python">
# 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.

</source>

Before proceeding, make sure that {{tclean_6.6.1}} has saved the model column. You should see a message like the following 
<pre style="background-color: #E0FFFF;">
>>> INFO .... ------ Predict Model ------
>>> INFO ... Saving model column
</pre>
If you don't see this message, you should repeat the above {{tclean_6.6.1}} command with the following modifications. This will populate the model column.
<pre style="background-color: #E0FFFF;">
niter=0
calcpsf=False
calcres=False
</pre>
As you continue with the self-cal, always make sure the model column is saved after running {{tclean_6.6.1}}.


Next you take that model and use it to determine the per-antenna phase solutions as a function of time using {{gaincal_6.6.1}}. We use the {{rmtables_6.6.1}} 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').

<source lang="python">
# 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)
</source>

Inspect the logging messages output by {{gaincal_6.6.1}} 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_6.6.1}} 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.

<pre style="background-color: #fffacd;">
Calibration solve statistics per spw: (expected/attempted/succeeded):
   Spw 0: 235/235/235
   Spw 1: 235/235/235
   Spw 2: 235/235/235
</pre>

[[File:Plotcal_image.png|thumb|Figure 1: Solutions from the first round of self-cal.]]

You should also check the solutions (above) that were obtained using {{plotms_6.6.1}}. 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. 

<source lang="python">
# in CASA
# Check the solution
plotms(vis='pcal1',
            xaxis='time',
            yaxis='phase',
            iteraxis='antenna',
            plotrange=[0,0,-180,180],
            gridrows=3,
            gridcols=3)
</source>

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.

[[File:Pcal1_6.5.4.png|thumb|Figure 2: Continuum image after the first round of self-cal.]]

<source lang="python">
# 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')
</source>

Save the flags in case you need to go back to this step.
<source lang="python">
# in CASA
flagmanager(vis=contvis,mode='save',versionname='after_pcal1')
</source>

Using our new gain solutions, we can generate an improved model for our source using {{tclean_6.6.1}}. 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.

<source lang="python">
# 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.
</source>

You should inspect the new image in CARTA to assess how well self-calibration worked. If using NRAO machines, you can open a new terminal tab, cd to the working directory, then type: 

<source lang="bash>
# in bash
carta --no_browser
</source>

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 [[Image_Line#Continuum_Subtraction_for_Line_Emission | 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_6.6.1}} 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.

<source lang="python">
# 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(vis='pcal2',
            xaxis='time',
            yaxis='phase',
            iteraxis='antenna',
            plotrange=[0,0,-180,180],
            gridrows=3,
            gridcols=3)

# 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')
</source>


Now we can repeat the {{tclean_6.6.1}}/{{gaincal_6.6.1}}/{{applycal_6.6.1}} 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.

<source lang="python">
# 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(vis='pcal3',
            xaxis='time',
            yaxis='phase',
            iteraxis='antenna',
            plotrange=[0,0,-180,180],
            gridrows=3,
            gridcols=3)

# 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.
</source>
[[File:Imaging-tutorial-selfcal-3.png|thumb|Figure 3: Amplitude based residuals after phase self-cal.]]
[[File:Apcal_6.5.4.png|thumb|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.

<source lang="python">
# 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(vis='apcal',
            xaxis='time',
            yaxis='phase',
            iteraxis='antenna',
            plotrange=[0,0,-180,180],
            gridrows=3,
            gridcols=3)

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)
</source>
	
Using CARTA, compare the final image to the initial image and see if the image dynamic range (ratio between 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. 

<source lang="python">
# in CASA
split(vis=contvis,
outputvis=contvis+'.selfcal', datacolumn='corrected')
</source>

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

<source lang="python">
# 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)
</source>
Then return to the {{applycal_6.6.1}} 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_6.6.1}} and {{delmod_6.6.1}} tasks and restore the original flags to return your ms to it’s original pre-self-cal state.

<source lang="python">
# 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)
</source>

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

'''[[Guide_NA_ImagingTemplate | Return to the Main Page]]'''