</aside>

<aside>

Download tutorial4_script_selfcal.py

This script was written for CASA 6.6.6 (ALMA pipeline version)
Source: First Look at Imaging (CASA Guide)
<https://casaguides.nrao.edu/index.php/First_Look_at_Self_Calibration_CASA_6.6.1>
Adapted for ALMA Data Reduction Workshop at Universidad de Chile 
Pietro Curone ([email protected]) 29 October 2025

</aside>

Dataset

TW Hya calibrated continuum measurement set (435 MB, already contained in the twhya_firstlook.tar you can download from the Before the workshop page)

twhya_calibrated.ms.tar (Download)

Your working directory should look like:

Tutorial4_selfcal >> ls
tutorial4_script_selfcal.py
twhya_calibrated.ms.tar

Untar with:

# In bash
tar -xvzf twhya_calibrated.ms.tar

Tutorial4_selfcal >> ls
tutorial4_script_selfcal.py
twhya_calibrated.ms
twhya_calibrated.ms.tar

Open CASA

Open a CASA terminal

Tutorial4_selfcal >> casa

Import useful libraries

# In CASA
import os
import numpy as np

Split out the science target (but no average in frequency or time!)

os.system('rm -rf twhya_science_target.ms')

split(vis='twhya_calibrated.ms',
      outputvis='twhya_science_target.ms',  
      datacolumn='data',
      field='TW Hya')
      
listobs(vis='twhya_science_target.ms',
        listfile='twhya_science_target.ms.txt',
        overwrite=True)

Create the initial model (Step 1)

Refer to Tutorial 2: Continuum Imaging for the imaging description. Here we are going to use the same imaging parameters, except for:

the nsigma threshold that we set higher (6 times the rms noise) to be more conservative and to include only the brightest regions in the model
the savemodel='modelcolumn' parameter, to save the CLEAN model.

Create the mask:

# Create the mask
mask_ra  = '11h01m51.837s'
mask_dec = '-34.42.17.216'
mask_pa  = 0  # mask position angle in degrees (here 0 as we want a circular mask)
mask_inc = 0   # mask inclination in degrees

mask_maj = 1.3  # semimajor axis of mask in arcsec
mask_min = mask_maj*np.cos(mask_inc * np.pi/180) #semiminor axis of mask in arcsec 

mask = f'ellipse[[{mask_ra}, {mask_dec}], [{mask_maj}arcsec, {mask_min}arcsec], {mask_pa}deg]'

CLEAN image p0, with p0 = “round 0 of phase-only selfcal”

imagename = 'twhya_cont_p0'  # p0 = "round 0 of phase-only selfcal"

os.system(f'rm -rf {imagename}.*') # Remove old versions in case you've run this before

tclean(vis='twhya_science_target.ms',
       imagename=imagename,
       specmode='mfs',
       gridder='standard',
       deconvolver='hogbom',
       cell=['0.1arcsec'],
       imsize=250,
       weighting='briggs',
       robust=0.5,
       mask=mask,    # Use the mask we just defined
       niter=1000000,    # Just choose a high number
       nsigma=6,    # CLEAN down to a threshold of 6 times the rms noise
       fastnoise=False,     # Better noise calculation
       interactive=False,
       savemodel='modelcolumn')

Create first set of gain solutions (Step 2)

refant (reference antenna): we do not actually measure absolute phase values; so we fix the phase of one antenna to zero, which gives us relative phases for all the others. We usually pick an antenna near the center of the array with as little flagging as possible. A good way to choose is to check the top 3–5 antennas ranked by the hif_refant task in the weblog.
gaintype: G finds solution for both polarizations independently, while T solves a common gain for both polarizations, since it averages them.
p1 = “round 1 of phase-only selfcal”

gaincal(vis="twhya_science_target.ms",
        caltable="gains_p1.cal",   
        solint="inf", # Notice this choice
        calmode="p",
        refant="DV22",
        gaintype="G")

We found solution:

Per antenna: always
Per polarization: because gaintype='G'
Per spectral window: default behavior unless combine='spw'
Per scan: because solint='inf'

Image solutions

plotms(vis='gains_p1.cal',
       xaxis='time',
       yaxis='phase',
       gridrows=3,
       gridcols=3,
       iteraxis='antenna',
       plotrange=[0,0,-50,50],
       coloraxis='corr',
       titlefont=7,
       xaxisfont=7,
       yaxisfont=7,
       plotfile='twhya_selfcal_p1_scan.png',
       showgui = False)

Screenshot 2025-10-29 at 00.07.16.png

Apply first set of gain solutions (Step 3)

**i**nterp tells CASA how to apply calibration solutions to visibilities at times where no solution was directly calculated, because solint in gaincal only gave solutions at specific intervals.

applycal(vis="twhya_science_target.ms",
         gaintable=["gains_p1.cal"], # The calibration table we just made
         interp="linear")

Split out corrected datacolumn (Step 4)

os.system('rm -rf twhya_science_target_p1.*')

split(vis="twhya_science_target.ms",
      outputvis='twhya_science_target_p1.ms',
      datacolumn='corrected')

New image (Step 5 of this round and step 1 of next round)

This image will show us if the SNR is getting better:

if SNR increases, we will proceed with a new round (p2, then p3, p4, etc.), decreasing the solint in gaincal at each step (e.g., 360s in p2, 120s in p3, 60s in p4, 30s in p5, 18s in p6). We’re now already producing the model for the p2 round, using savemodel='modelcolumn' .
if the SNR stays the same or decreases, we stop here and take the previous step’s image (p0 in this case) as the final one.

imagename = 'twhya_cont_p1'  

os.system(f'rm -rf {imagename}.*') 

tclean(vis='twhya_science_target_p1.ms',  # Notice new p1 ms
       imagename=imagename,
       specmode='mfs',
       gridder='standard',
       deconvolver='hogbom',
       cell=['0.1arcsec'],
       imsize=250,
       weighting='briggs',
       robust=0.5,
       mask=mask,    # Use the same mask as before
       niter=1000000,    # Just choose a high number
       nsigma=6,    # CLEAN down to a threshold of 6 times the rms noise
       fastnoise=False,     
       interactive=False,
       savemodel='modelcolumn')

p0 image:

peak intensity: 288 mJy
rms noise: 8 mJy
SNR = 36

p1 image:

peak intensity: 316 mJy
rms noise: 4.8 mJy
SNR: 113

SNR increased! We could continue with next rounds.

Screenshot 2025-10-29 at 00.17.31.png

Auto selfcal (only works on Linux)

<aside>

Described and stored in this GitHub repository

https://github.com/jjtobin/auto_selfcal

</aside>

It can only work on Linux as it needs the CASA Viewer, unavailable on MacOS.

This code should only be executed within the CASA 6.4+ environment.

Codes

Create a new directory (Tutorial4_autoselfcal) and clone the GitHub repo

Tutorial4_autoselfcal >> git clone <https://github.com/jjtobin/auto_selfcal.git>

Copy the Python script from the repo into your working directory

Tutorial4_autoselfcal >> cp -r auto_selfcal/*py .

Dataset

TW Hya calibrated continuum measurement set (435 MB, already contained in the twhya_firstlook.tar you can download from the Before the workshop page)

Logistics

Hands-on tutorials

Lectures

Dataset

Open CASA

Split out the science target (but no average in frequency or time!)

Create the initial model (Step 1)

Create first set of gain solutions (Step 2)

Apply first set of gain solutions (Step 3)

Split out corrected datacolumn (Step 4)

New image (Step 5 of this round and step 1 of next round)

Auto selfcal (only works on Linux)

Codes

Dataset