Self-Calibration Template CASA 6.1.1: Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
m (Updated plotms commands)
(101 intermediate revisions by 5 users not shown)
Line 1: Line 1:
This guide continues with products created in '''[[Image_Continuum]]'''. All imaging parameters are set in the previous guide. You will need '''calibrated_final_cont.ms''' to proceed.
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].
 
== Self-Calibration on the Continuum (optional) ==
== Self-Calibration on the Continuum (optional) ==


This step of Imaging is optional and can be attempted in cases where the expected rms is not reached on the continuum or the line data. Self calibration solutions can only be determined when a source exhibits a strong continuum emission, preferably near the phase center. With this being said, self calibration should not be performed on data sets with poor uv coverage but can be performed on “partial” data sets (like ACA _TC datasets). It is not recommended on projects that are non-detections or contain weak sources.  
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}}. 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].


For a more detailed explanation of self-calibration, please see Ed Fomalont and Jeff Mangum’s [https://staff.nrao.edu/wiki/pub/NAASC/DataReductionDiscussion/SelfcalBestPractices.pdf Best Practices for Self Calibration of Spectral Line Image Cubes] '''Editing Note: this is an internal to NRAO link'''
----


Begin by declaring your ms and continuum image name.  
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">
# in CASA
# in CASA
contvis = 'calibrated_final_cont.ms'        
flagmanager(vis=contvis,mode='save',versionname='before_selfcal',merge='replace')
contimagename = 'calibrated_final_cont'
</source>
</source>


Next, choose a reference antenna that is in the array by using the tasks plotants or listobs/vishead. For data sets with multiple executions, you will want to choose an antenna that is present in all of the executions, which can be found by using the task au.commonAntennas(). You should use the same reference antenna that was used during calibration.  
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">
<source lang="python">
# in CASA
# in CASA
refant = 'DV09' # pick a reference antenna.
contvis = 'calibrated_final_cont.ms'       
contimagename = 'calibrated_final_cont'
</source>
</source>


You will need to map the self-calibration solutions to individual spectral windows. This is generally an array of n zeros, where n is the number of spectral windows for the ms. This is the case because spwcombine is used, which combines all solutions to a single spectral window numbered 0. By setting an array of n zeros we are applying the solution to each spectral window within the ms.
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">
# in CASA
# in CASA
spwmap = [0,0,0] # change to an array of n zeros
# delmod(vis=contvis,otf=True,scr=True)
# clearcal(vis=contvis)
</source>
</source>


Save all of the initial flags in case the final self calibration solution still doesn’t meet your expectations. The task applycal() will flag data that does not have solutions.  
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">
# in CASA
# in CASA
flagmanager(vis=contvis,mode='save',versionname='before_selfcal',merge='replace')
refant = 'DV09' # pick a reference antenna.
</source>
</source>


In the case of running self-calibration 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 script).
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">
<source lang="python">
# in CASA
# in CASA
# delmod(vis=contvis,field=field,otf=True)
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.
# clearcal(vis=contvis)
</source>
</source>


You then begin the process of shallowly cleaning your continuum. Once the parameters for cleaning are set, and clean is run interactively, the next step is to calibrate antenna based gains as a function of time (gaincal) and apply solutions for next cleaning. With each additional clean, the time iterations will get shorter. At this point you should clean until satisfied with the end result.  
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">
<source lang="python">
Line 49: 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>
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 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').
<source lang="python">
# in CASA
# per scan solution
# per scan solution
rmtables('pcal1')
rmtables('pcal1')
Line 84: Line 136:
</source>
</source>


You may see messages about failed solutions when running gaincal. Any data that does not have solutions in the calibration table will be flagged when the calibration is applied to the ms. If too much data is flagged, the sensitivity will decrease. Look for messages in the logger about the success in of the calculations. The image above shows the expected/attempted/succeeded numbers in the logger window. It is not a problem to have a few failed solutions, but be cautious if it seems many are failing. Note that you will only see successful calculations for spw 0.
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.
 
<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>


The minimum snr (minsnr) and minimum baseline per antenna (minblperant) parameters can be decreased to increase the number of solutions found.
[[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 plotcal().  
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',
plotms(caltable='pcal1',
         xaxis='time',
         xaxis='time',
         yaxis='phase',
         yaxis='phase',
        timerange='',
         iteration='antenna',
         iteration='antenna',
        subplot=421,
         plotrange=[0,0,-180,180])
         plotrange=[0,0,-180,180])
</source>
</source>


If you are satisfied with the solutions, apply them to the ms and clean the new visibilities.
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.png|thumb|Figure 2: Continuum image after the first round of self-cal.]]


<source lang="python">
<source lang="python">
Line 115: Line 175:
         flagbackup=False,
         flagbackup=False,
         interp='linearperobs')
         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}}. 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
# 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>


Inspect the new image in the viewer to see if the first iteration of self-cal worked. Find and compare the new RMS, beam size, and dynamic range of the image to the calibrated_final_cont_p0.image. If the RMS noise increased, do not proceed with the self calibration. You may try changing the minsnr and minblperant parameters in the gaincal call above to find more solutions. If you are unable to improve the image using self-calibration, undo the self-calibration and skip to the continuum subtraction section of this guide.
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 RMS noise and dynamic range improved, proceed with self-calibration below. The next gaincal() call uses a shorter solution interval (solint). Continue on with the self-calibration as long as the image statistics continue to improve.
 
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.


<source lang="python">
<source lang="python">
Line 158: Line 235:


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


Line 175: Line 250:
         flagbackup=False,
         flagbackup=False,
         interp='linearperobs')
         interp='linearperobs')
flagmanager(vis=contvis,mode='save',versionname='after_pcal2')
</source>
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">
# 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 212: Line 299:


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


Line 230: Line 315:
         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>
[[File:Imaging-tutorial-selfcal-3.png|thumb|Figure 3: Amplitude based residuals after phase self-cal.]]
[[File:Apcal.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


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


<source lang="python">
<source lang="python">
Line 272: Line 369:
         solnorm=True)
         solnorm=True)


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


Line 288: Line 383:
         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>
Compare all of the images produced during self-calibration. Only apply solutions that continue to improve the image statistics.  
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.


The results of the self-cal should be split out via the split command and saved as a *.splitcal file.
Finally, you should split out the results of your self-calibration for safe-keeping.  


<source lang="python">
<source lang="python">
Line 319: Line 424:
</source>
</source>


If you are unhappy with the self-calibration, use the clearcal task to return your ms to it’s original pre-self-cal state. In addition to the flagmanager task.  
== 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 328: Line 446:
# delmod(contvis,field=field,otf=True)
# delmod(contvis,field=field,otf=True)
</source>
</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]]'''

Revision as of 08:49, 24 August 2021

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