Self-Calibration Template CASA 6.1.1: Difference between revisions

From CASA Guides
Jump to navigationJump to search
(47 intermediate revisions by 3 users not shown)
Line 4: Line 4:
== Self-Calibration on the Continuum (optional) ==
== 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 S/N of at least 45 in your continuum image to attempt self-calibration.
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.
<center>
<equation id="SNR">
<math>\ V_{XX} = (I + Q_{\psi}) + U_{\psi}(d_{Xj}^{*} + d_{Xi})</math>


<math> S/R </math>
<math>
\frac{Peak}{RMS} > 3 \sqrt{N-3} \sqrt{\frac{t_{int}}{t_{solint}}}
</math>


</equation>
<xr id="SNR" nolink/>. Use this equation to determine the minimum SNR in order to attempt self-calibration.
</center>
<math>Peak/RMS</math>


<math>\frac{Peak}{RMS} > 3 * {{math|{{radical|N-3}}}} * {{math|{{radical|\frac{Tint}{Tscan}}}}}</math>
where:


Self-calibration is an iterative process. You start by generating a model of your source using {{clean}}. 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 you run out of 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.
''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}}. 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 [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].
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].
For a more detailed explanation of self-calibration, please 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].


----
----


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


<source lang="python">
<source lang="python">
Line 41: Line 47:
</source>
</source>


If you run the self-calibration process more than once, you will need to remove any models that could be left in the image header; which can be accomplished by using the {{delmod}} task as well as the {{clearcal}} task (which will need to be uncommented in the code below).
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).


<source lang="python">
<source lang="python">
Line 49: Line 55:
</source>
</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'''. 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}} or {{listobs}}/{{vishead}} can also be used to find antennas in all executions.  
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}} or {{listobs}}/{{vishead}} can also be used to find antennas in all executions.  


<source lang="python">
<source lang="python">
Line 59: Line 65:


<source lang="python">
<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.
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>
</source>
Line 69: Line 76:
# shallow clean on the continuum
# shallow clean on the continuum


for ext in ['.flux','.image','.mask','.model','.pbcor','.psf','.residual','.flux.pbcoverage']:
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
     rmtables(contimagename + '_p0'+ ext)
     rmtables(contimagename + '_p0'+ ext)
   
 
clean(vis=contvis,
tclean(vis=contvis,
      imagename=contimagename + '_p0',
      imagename=contimagename + '_p0',
      field=field,
      field=field,
#     phasecenter=phasecenter, # uncomment if mosaic.     
      #phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
      mode='mfs',
      # mosweight = True, # uncomment if mosaic
      psfmode='clark',
      specmode='mfs',
      imsize = imsize,
      deconvolver='hogbom',
      cell= cell,
      # Uncomment the below to image with nterms>1.
      weighting = weighting,
      #deconvolver='mtmfs',
      robust=robust,
      #nterms=2,
      niter=niter,
      imsize = imsize,  
      threshold=threshold,
      cell= cell,  
      interactive=True,
      weighting = weighting,  
      usescratch=True, # needed for 4.3 and 4.4 (and maybe 4.5)
      robust=robust,
      imagermode=imagermode)
      niter=niter,  
      threshold=threshold,  
      interactive=True,
      gridder=gridder,
      savemodel='modelcolumn',
      usepointing=False)


# Note number of iterations performed.
# Note number of iterations performed.


</source>
</source>
Before proceeding, make sure that {{tclean}} 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}} 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 runnint {{tclean}}.




Line 96: Line 121:


<source lang="python">
<source lang="python">
# in CASA
# per scan solution
# per scan solution
rmtables('pcal1')
rmtables('pcal1')
Line 122: Line 148:
</figure>
</figure>


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


<source lang="python">
<source lang="python">
In CASA
# in CASA
# Check the solution
# Check the solution
plotcal(caltable='pcal1',
plotcal(caltable='pcal1',
Line 137: Line 164:


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.
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 id="Pcal1.png">
<figure id="Pcal1.png">
[[File:Pcal1.png|thumb|<caption>Continuum image after the first round of self-cal.</caption>]]
[[File:Pcal1.png|thumb|Figure 2: Continuum image after the first round of self-cal.]]
</figure>
</figure>


Line 155: Line 181:
</source>
</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 {{clean}}. During this second clean 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.
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.


<source lang="python">
<source lang="python">
# in CASA
# clean deeper
# clean deeper
for ext in ['.flux','.image','.mask','.model','.pbcor','.psf','.residual','.flux.pbcoverage']:
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
     rmtables(contimagename + '_p1'+ ext)
     rmtables(contimagename + '_p1'+ ext)


clean(vis=contvis,
tclean(vis=contvis,
      field=field,
      imagename=contimagename + '_p1',
#     phasecenter=phasecenter, # uncomment if mosaic.     
      field=field,
      imagename=contimagename + '_p1',
      # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
      mode='mfs',
      # mosweight = True, # uncomment if mosaic
      psfmode='clark',
      specmode='mfs',
      imsize = imsize,
      deconvolver='hogbom',
      cell= cell,
      # Uncomment the below to image with nterms>1.
      weighting = weighting,
      #deconvolver='mtmfs',
      robust=robust,
      #nterms=2,
      niter=niter,
      imsize = imsize,  
      threshold=threshold,
      cell= cell,  
      interactive=True,
      weighting = weighting,  
      usescratch=True, # needed for 4.3 and 4.4 (and maybe 4.5)
      robust=robust,
      imagermode=imagermode)
      niter=niter,  
 
      threshold=threshold,  
      interactive=True,
      gridder=gridder,
      #pbcor = True, #if final image
      savemodel='modelcolumn',
      usepointing=False)
# Note number of iterations performed.
# Note number of iterations performed.
</source>
</source>


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, clear the calibration with clearcal (see commands above) and proceed to [[Image_Line#Continuum_Subtraction_for_Line_Emission | continuum subtraction]] (if you are working on a line project).
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 [[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}} 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.
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.
Line 219: Line 256:
         flagbackup=False,
         flagbackup=False,
         interp='linearperobs')
         interp='linearperobs')
flagmanager(vis=contvis,mode='save',versionname='after_pcal2')
</source>
</source>




Now we can repeat the {{clean}}/{{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 clean 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.
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.


<source lang="python">
<source lang="python">
# in CASA
# clean deeper
# clean deeper
for ext in ['.flux','.image','.mask','.model','.pbcor','.psf','.residual','.flux.pbcoverage']:
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
     rmtables(contimagename + '_p2'+ ext)
     rmtables(contimagename + '_p2'+ ext)


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


Line 279: Line 323:
         interp='linearperobs')
         interp='linearperobs')


for ext in ['.flux','.image','.mask','.model','.pbcor','.psf','.residual','.flux.pbcoverage']:
flagmanager(vis=contvis,mode='save',versionname='after_pcal3')
 
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt']:
     rmtables(contimagename + '_p3'+ ext)
     rmtables(contimagename + '_p3'+ ext)


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


# Note number of iterations performed.
# Note number of iterations performed.
</source>
</source>
<figure id="Imaging-tutorial-selfcal-3.png">
[[File:Imaging-tutorial-selfcal-3.png|thumb|Figure 3: Amplitude based residuals after phase self-cal.]]
</figure>
<figure id="Apcal.png">
<figure id="Apcal.png">
[[File:Apcal.png|thumb|<caption>Image after amplitude self-cal.</caption>]]
[[File:Apcal.png|thumb|Figure 4: Image after amplitude self-cal.]]
</figure>
</figure>


Line 342: Line 397:
         flagbackup=False,
         flagbackup=False,
         interp='linearperobs')
         interp='linearperobs')
flagmanager(vis=contvis,mode='save',versionname='after_apcal')


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


clean(vis=contvis,
 
      imagename=contimagename + '_ap',
tclean(vis=contvis,
      field=field,
      imagename=contimagename + '_ap',
#     phasecenter=phasecenter, # uncomment if mosaic.     
      field=field,
      mode='mfs',
      # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
      psfmode='clark',
      # mosweight = True, # uncomment if mosaic
      imsize = imsize,
      specmode='mfs',
      cell= cell,
      deconvolver='hogbom',
      weighting = weighting,
      # Uncomment the below to image with nterms>1.
      robust=robust,
      #deconvolver='mtmfs',
      niter=niter,
      #nterms=2,
      threshold=threshold,
      imsize = imsize,  
      interactive=True,
      cell= cell,  
      usescratch=True, # needed for 4.3 and 4.4 (and maybe 4.5)
      weighting = weighting,  
      imagermode=imagermode)
      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>
</source>
Line 374: Line 438:
</source>
</source>


If you are unhappy with the self-calibration, use the {{clearcal}} task and restore the original flags to return your ms to it’s original pre-self-cal state.
== 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}} 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.


<source lang="python">
<source lang="python">
Line 383: Line 460:
# delmod(contvis,field=field,otf=True)
# delmod(contvis,field=field,otf=True)
</source>
</source>
[xxx AAK there should be also a note on the best way to restart the self-cal process. This is somewhat complicated because it requires keeping track of the data, model, and corrected data columns as well as the flags. the easiest way to restart might be to put a split after the apply cal step in each round of self cal and do a copy of the ms before the first stepk but I need to think more about this.]


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]]'''.
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]]'''
'''[[Guide_NA_ImagingTemplate | Return to the Main Page]]'''

Revision as of 15:09, 11 December 2019

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.

tint is the total on source integration time.

tsolint 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,field=field,otf=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 runnint 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 id="Plotcal_image.png">

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

</figure>

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
plotcal(caltable='pcal1',
        xaxis='time',
        yaxis='phase',
        timerange='',
        iteration='antenna',
        subplot=421,
        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 id="Pcal1.png">

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

</figure>

# 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
plotcal(caltable='pcal2',
        xaxis='time',
        yaxis='phase',
        timerange='',
        iteration='antenna',
        subplot=421,
        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
plotcal(caltable='pcal3',
        xaxis='time',
        yaxis='phase',
        timerange='',
        iteration='antenna',
        subplot=421,
        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 id="Imaging-tutorial-selfcal-3.png">

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

</figure> <figure id="Apcal.png">

Figure 4: Image after amplitude self-cal.

</figure>

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)

plotcal(caltable='apcal',
        xaxis='time',
        yaxis='amp',
        timerange='',
        iteration='antenna',
        subplot=421,
        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