Part 1: CASA Basic EVN Calibration

« Return to homepage
« Return to EVN continuum

Guidance for calibrating EVN data using CASA

This guide uses capabilities of CASA 4.7.2; CASA 5+ does work apart from the Tsys calibration step (and the script isn’t fixed yet, see Issue #1).

It is suitable for good-quality EVN data which does not require the use of ‘rate’ in calibration. It also requires a prepared gaincurve table and helper scripts if you are starting from the fitsidi data. You need about 20 G disk space, although one could manage with less by moving early versions of these files to storage.

This section is organised as follows:

  1. Data and supporting material
  2. How to use this guide
  3. Data loading and inspection
  4. Frequency-related calibration
  5. Apply calibration and split out phase-ref - target pairs

1. Data and supporting material

« back to top

n14c3 is an EVN network monitoring experiment. It contains two target - phase-reference pairs of bright sources (plus a 5th source which we will not use). The full data reduction here needs:

These are all provided on your desktop for the DARA workshops, but the fits IDI files can also be obtained from here along with more information about the observations. The other files are available in the tar file: NME_DARA.tgz

These are the sources:

Number Name Position (J2000) Role
0 J1640+3946 16:40:29.632770 +39.46.46.02836 target
1 3C345 16:42:58.809965 +39.48.36.99402 phase-cal
2 J1849+3024 18:49:20.103406 +30.24.14.23712 target
3 1848+283 18:50:27.589825 +28.25.13.15523 phase-cal
4 2023+336 20:25:10.842114 +33.43.00.21435 not used

J1849+3024 is also used as bandpass calibrator because it is a bright, compact source with good data on all baselines. If you are reducing a new data set, you would check the suitability of the bandpass calibrator during early data reduction.

2. How to use this guide

« back to top

This web page presents inputs for CASA tasks to be entered interactively at a terminal prompt for the calibration of the averaged data and initial target imaging. This is useful to examine your options and see how to check for syntax errors.

To start:

  1. Make a directory called something you will remember, e.g. EVN, and work in it. (Do not do data reduction within the CASA installation).
  2. Copy the data and files you need to this directory. <path***> is the original location of these files; Note in general, *** means something for you to fill in. e.g:
mdir EVN
cd EVN
cp <path***>/n14c3_1_1.IDI? .
cp <path***>/NME_DARA.tgz .
cp <path***>/*py .

tar -zxvf NME_DARA.tgz          # extract the additional material (also the scripts  which will be supplied).

3. Data loading and inspection

« back to top

pwd                   # tells you present working directory
ls
which should show

flagSH1.flagcmd      key.py          n14c3.antab        NME_3C345_skeleton.py  NME_J1849.py
flagSH.flagcmd       n14c3_1_1.IDI1  n14c3.gc           NME_all.py             tsys_spectrum.py
flagTar1Ph1.flagcmd  n14c3_1_1.IDI2  NME_3C345.py       NME_DARA.tgz

a. Load data and correct for Earth rotation

« back to top

# in CASA
os.system('rm -rf n14c3_prefix.ms*')              # make sure no old files are present

default(importfitsidi)                            # initialise task
inp()                                             # list all inputs

fitsidifile=['n14c3_1_1.IDI1','n14c3_1_1.IDI2']   # fits files to convert
vis='n14c3_prefix.ms'                            # output ms file name
constobsid=True                                   # data could be processed together
scanreindexgap_s=15                               # only separate scans if gap >15 sec (or source change)
inp()                                             # check your inputs

importfitsidi()                                   # if happy, execute the task

inp() provide checks that you should normally carry out but the script won’t repeat these for every task.

Note A few warnings about DIAMETER and POLCALA/B and zero/negative scan Nos may appear because CASA is not yet fully adapted for VLBI data but normally can be ignored if the task finishes ok. Check you CASA logger just in case.

Important You can review the inputs to the last time you ran a task by typing:

# In CASA:
ls      # This should show the MS you just created and importfitsidi.last

!more importfitsidi.last   # ! allows you to run any linux shell command
This will show something like:

taskname           = 'importfitsidi'
fitsidifile        =  ['n14c3_1_1.IDI1', 'n14c3_1_1.IDI2']
vis                =  'n14c3__prefix.ms'
constobsid         =  True
scanreindexgap_s    =  15
#importfitsidi(fitsidifile=['n14c3_1_1.IDI1', 'n14c3_1_1.IDI2'],vis="n14c3__prefix.ms",constobsid=True,scanreindexgap_s=15)

The first format is for the interactive input (i.e. what you type directly into CASA), and the second format (without ‘#’) is for use when scripting (automatic calibration! woo!).

# In CASA:
os.system('rm -rf n14c3.ms')
os.system('rm -rf n14c3.ms.flagversions')

default(fixvis)
vis='n14c3_prefix.ms'
outputvis='n14c3.ms'

fixvis()

# In CASA:
os.system('rm -rf n14c3.ms.listobs')

default(listobs)
vis='n14c3.ms'
listfile='n14c3.ms.listobs'
listobs()

Important The listobs output give information on what is contained in your data set i.e. antennas, sources, scans etc. The following section breaks this file down:

b. Fix antenna and Tsys tables

« back to top

These steps may not always be needed as CASA and EVN data formats become more compatible especially with the VLBI tools getting released with CASA v5.3+.

The antab table has to be converted into a subtable within the MS using tsys_spectrum.py. This is then be used to generated a Tsys correction table using gencal. The system temperature is measured every few minutes. This compensates roughly for the different levels of signal from different sources, the effects of elevation and other amplitude fluctuations. Each antenna has a characteristic response in terms of Kelvin of system temperature per Jy of flux, so Tsys is also used to scale the flux density.

# In CASA:
inbase='n14c3'                   # pass base name to script
execfile('tsys_spectrum.py')     # execfile runs any suitable script inside CASA

os.system('rm -rf n14c3.tsys')

default(gencal)  # This reads the Tsys information we just wrote into the MS and derives a calibration table
vis='n14c3.ms'
caltable='n14c3.tsys'
caltype='tsys'

gencal()

Plot the table generated.

# In CASA:
default(plotcal)
caltable='n14c3.tsys'
xaxis='freq'
yaxis='tsys'
subplot=431
iteration='antenna'
plotcal()

alt text

You can see the changes with each source change.

Three antennas do not have Tsys measurements but the .gc (gain curve) table provides a scaled gain-elevation correction which scales the visibility amplitudes, although without any allowance for weather or source contribution

listobs gave the antenna diameters as zero since the information is not written into the fits files for EVN data. Some antennas have off-axis feeds and so the offsets also have to be inserted.

# In CASA
# Copy these arrays of values (no line breaks in each array)
ants=  ['EF','WB','JB','ON','NT','TR','SV','ZC','BD','SH','HH','YS','JD']
diams= [100.0,300.0,75.0,25.0,32.0,32.0,32.0,32.0,32.0,25.0,24.0,40.0,25.0]
axoffs=[[0.013,4.95,0.,2.15,1.831,0.,-0.007,-0.008,-0.004,-0.002,6.692,2.005,0.]
        ,[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]
        ,[0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.]]

# Modify the antenna table
tb.open('n14c3.ms/ANTENNA', nomodify=F)
tb.putcol('DISH_DIAMETER', diams)
tb.putcol('OFFSET',axoffs)
tb.close()

default(listobs)     # Check in listobs
vis='n14c3.ms'
listobs()

This time we don’t write a text file because the parameter listfile was not set, so just look in the logger to ensure that the antenna diameters now appear.

c. Inspect and flag data

« back to top

# In CASA
default(flagdata)
vis='n14c3.ms'
mode='manual'
autocorr=T

flagdata()

# In CASA
default(plotants)
vis='n14c3.ms'

plotants()

alt text

# In CASA
default(plotms)
vis='n14c3.ms'
xaxis='frequency'
yaxis='amp'
field='1848+283'        # The phase-cal/bandpass cal
avgtime='3600'          # Will only average within scans unless additionally told to average scans too
antenna='EF&*'          # Plot all baselines to the largest and most sensitive antenna
correlation='RR,LL'     # Only plot the parallel hands; the cross hands are fainter (and we won't be using them).
coloraxis='antenna2'

plotms()

alt text

  # In CASA
  default(flagmanager)
  vis='n14c3.ms'
  mode='save'
  versionname='preSVandEndChans'

  flagmanager()

  # In CASA
  default(flagdata)
  vis='n14c3.ms'
  mode='manual'
  antenna='SV'

# In CASA
default(plotms)
vis='n14c3.ms'
xaxis='channel'
yaxis='amp'
field='1848+283'        # The phase-cal/bandpass cal
avgtime='3600'          # Will only average within scans unless additionally told to average scans too
antenna='EF&JB'          # Plot all baselines to the largest and most sensitive antenna
correlation='RR,LL'     # Only plot the parallel hands; the cross hands are fainter (and we won't be using them).
iteraxis='spw'

plotms()

# In CASA
default(flagdata)
vis='n14c3.ms'
mode='manual'
spw='0:0~5;29~31,2:0~5;29~31,4:0~5;29~31,6:0~5;29~31,1:0~2;27~31,3:0~2;27~31,5:0~2;27~31,7:0~2;27~31')

flagdata()

# In CASA
default(plotms)
vis='n14c3.ms'
xaxis='time'
yaxis='amp'
field='1848+283'        
spw='0~7:13~20'        # Average a few central channels where the response is stable
avgchannel='8'  
antenna='EF&*'          
correlation='RR,LL'     
coloraxis='baseline'

plotms()


# In CASA
default(flagdata)
vis='n14c3.ms'
mode='quack'
quackinterval = 5     # 5 seconds, or just over 2 intgrations

flagdata()

# In CASA
default(plotms)

vis='n14c3.ms'
xaxis='time'
yaxis='amp'
scan='60~70'            # a few scans including the bad data; plot all fields in this time range
spw='0~7:13~20'       
avgchannel='8'  
antenna='EF&HH'         # just the relevant baseline          
correlation='RR,LL'     
coloraxis='field'       # distinguish sources

plotms()


# In CASA
default(flagdata)
vis='n14c3.ms'
antenna='HH'
mode='manual'
scan='62~65'

flagdata()

# In CASA
default(plotms)

vis='n14c3.ms'
xaxis='time'
yaxis='amp'
field='1848+283'
spw='0~7:13~20'       
avgchannel='8'  
antenna='EF&SH'         # just the relevant baseline          
correlation='RR,LL'     
coloraxis='spw'         # distinguish spw

plotms()

There are many short, bad periods, some only affecting certain spw. The easiest way to cope with this is to write a flag list. Identify antennas/spw/times, adding 1s to the start and finish of each flagging period. They are all shorter than one scan so the same flags aren’t applied to the target but later, we will look for similar bad data on the target. Take a look at the flag list file: flagSH.flagcmd which is in your current directory and cross-check that these entered values are affected.

If you are happy, then we use this file to automatically flag those channels specified in the file:

# In CASA
default(flagdata)
vis='n14c3.ms'
mode='list'
inpfile='flagSH.flagcmd'

flagdata()

# In CASA
default(plotms)
vis='n14c3.ms'
xaxis='time'
yaxis='amp'
field='1848+283'
spw='0~7:13~20'
avgchannel='8'
antenna='EF&*'
coloraxis='corr'
correlation='RR,LL'

plotms()

a. Delay calibration

« back to top

The signal from each VLBI antenna is timestamped in order to allow synchronisation of all antennas when they are correlated. Small clock errors cause a frequency-dependent phase shift in the correlated signal, i.e. a slope of the phase across the band.

# In CASA
default(plotms)
vis='n14c3.ms'
xaxis='freq'
yaxis='phase'
field='1848+283'
avgtime='3600s'
antenna='EF&*'
coloraxis='corr'
iteraxis='baseline'
correlation='RR,LL'

plotms()

The phases have a random wiggle (which will be corrected later) as well as a systematic slope, which is more or less the same for each scan.

# In CASA
default(plotms)
vis='n14c3.ms'
xaxis='scan'
yaxis='amp'
field='1848+283'
spw='0~7:13~20'
avgchannel='8'
avgtime='3600'
antenna='EF&*'
coloraxis='antenna2'
correlation='RR,LL'

plotms()

Important: Delay calibration derives the gradient of phase against frequency, which can divided into the complex data to remove the systematic slope.

# In CASA
default(gaincal)
vis='n14c3.ms'
caltable='n14c3_bpcal.K'  # File name for output calibration table
field='1848+283'
scan='38'
solint='150s'             # longer that the scan, to make sure just one solution
refant='EF'               # The most sensitive antenna, with many shorter baselines, as the zero of phase.
gaintype='K'              # Solve for delay
parang=T                  # Apply the parallactic angle correction for the feed rotation of Alt-Az telescopes

gaincal()

default(plotcal)
caltable='n14c3_bpcal.K'
xaxis='freq'
yaxis='delay'
iteration='antenna'
subplot=431

plotcal()

The solutions are smooth in frequency and don’t have massive values (30ns is ok for EVN observations). Note that the delay solutions are proportional to the baseline length (see the plotants image for reference)

b. Pre-bandpass time-dependent phase correction

« back to top

We need now need to correct channel-by-channel errors but in order to get high enough signal to noise, we need to average the data in time. However, in order to do that, we need to average in time. The delay corrections just derived are good enough to allow preliminary channel averaging to derive time-dependent solutions before doing this.

# In CASA
default(plotms)
vis='n14c3.ms'
xaxis='time'
yaxis='phase'
field='1848+283'
spw='0~7:13~20'
avgchannel='8'
antenna='EF&*'
coloraxis='spw'
iteraxis='baseline'
correlation='RR'

plotms()

And zooming into a single scan reveals:

Some notes on these plots:

If we are happy, then we use this solution interval and the task gaincal to derive the time dependent phase solutions. Note that gaincal compares the phase on each baseline with the zero phase for a point source model and uses a chi-squared minimisation method to decompose this into a per-antenna correction.

# In CASA

os.system('rm -rf '+inbase+'n14c3_bpcal.p0')
default (gaincal)
vis='n14c3.ms'
caltable='n14c3_bpcal.p0'             
field='1848+283'
solint='30s'          
refant='EF'
gaintype='G'
calmode='p'                     # Solve for phase only as a function of time
gaintable=['n14c3_bpcal.K']     # Apply the delay corrections
interp='nearest'   
parang=T

gaincal()

Two solutions fail for 2 spectral windows in one solution interval, which is OK to ignore because that is just a small fraction of the time on this source.

# In CASA

default(plotcal)
caltable='n14c3_bpcal.p0'
xaxis='time'
yaxis='phase'
iteration='antenna'
subplot=431

plotcal()

For each spw, in a different color, the corrections have a systematic shape. These solutions are also smooth in time which means (normally) that these are ok!

c. Bandpass correction

« back to top

We use the bandpass calibration source 1848+283 to derive channel-by-channel corrections for phase and amplitude. We apply the phase vs. time calibration table ()n14c3_bpcal.p0) to allow the data to be averaged in time, and we apply the delay and Tsys/gain curve correction tables to provide preliminary per-spw corrections to the slope of the phase and the overall amplitude fluctuations with time.

Since we have not provided a different model, bandpass uses a default 1 Jy model at the phase centre; hence we normalise the amplitude corrections, so there are different corrections as needed per channel but they are normalised by the overall product (over all time on this source but in each spw).

#In CASA

os.system('rm -rf n14c3.B')
default(bandpass)
vis='n14c3.ms'
caltable='n14c3.B'           # output bandpass calibration table
field='1848+283'
solint='inf'
combine='scan'               # average all the data in time
refant='EF'    
solnorm=T                    # Keep the product of corrections unity so the flux scale is not changed overall
gaintable=['n14c3_bpcal.K','n14c3_bpcal.p0','n14c3.tsys','n14c3.gc']
interp='nearest'
parang=T

bandpass()

The solutions fail for the end channels we have flagged but the rest are OK.

# In CASA
default(plotcal)

caltable='n14c3.B'
xaxis='freq'
yaxis='phase'
iteration='antenna'
subplot=431

plotcal()

xaxis='freq'
yaxis='amp'
plotcal()

The phase corrections are mostly small, larger close to the remaining end channels, and the amplitude corrections are close to one.

Apply calibration and split out phase-ref - target pairs

Apply the initial calibration

« back to top

The frequency-dependent calibration i.e. delay, Tsys/gain curve and bandpass, is stable in time over the few hours of these observations and we apply these table to all sources.

Note: we don’t apply the time-dependent phase corrections for the bpcal, because these are not applicable to other sources observed at different times (usually a separate bandpass calibrator is observed before or after the main observations; these are unusual because it is also a phase-cal, but it is convenient to ignore this role here).

# In CASA
default(applycal)
vis='n14c3.ms'
field=''                      # Apply to all sources
gaintable=['n14c3.tsys','n14c3.gc','n14c3_bpcal.K','n14c3.B']
interp='nearest'
parang=T
applymode='calonly'           # Interpolate or extrapolate over data with missing solutions (don't flag)

applycal()

default(plotms)
vis='n14c3.ms'
xaxis='freq'
yaxis='amp'
field='1848+283'
ydatacolumn='corrected'   # Applying calibration creates a new 'corrected' data column
correlation='RR'
coloraxis='baseline'
avgtime='3600'
antenna='EF&*'

plotms()

yaxis='phase'
plotms()

xaxis='time'
avgtime=''
avgchannel='32'
plotms()

field=''
coloraxis='field'
yaxis='amp'
plotms()

Amplitude and phase are reasonably flat as a function of frequency but different for different antennas and the phase still has a big variation with time, along with some amplitude discrepancies.

Split out each pair of sources

« back to top

For convenience, split out each phase-reference - target pair separately. The ‘corrected’ column will just be ‘data’ in the newMSs.

We use the task split to do this:

# In CASA
os.system('rm -rf J1640+3946_3C345.ms')
os.system('rm -rf J1640+3946_3C345.ms.flagversions')

default(split)
vis='n14c3.ms'
outputvis='J1640+3946_3C345.ms'
field='J1640+3946,3C345'
datacolumn='corrected'

split()

os.system('rm -rf 1848+283_J1849+3024.ms')
os.system('rm -rf 1848+283_J1849+3024.ms.flagversions')
outputvis='1848+283_J1849+3024.ms'
field='1848+283,J1849+3024'
datacolumn='corrected'

split()

Ignore any ‘file left unfinished’ messages if the lists are OK)

os.system('rm -rf J1640+3946_3C345.ms.listobs')
default(listobs)
vis='J1640+3946_3C345.ms'
listfile='J1640+3946_3C345.ms.listobs'
listobs()

os.system('rm -rf 1848+283_J1849+3024.ms.listobs')
vis='1848+283_J1849+3024.ms'
listfile='1848+283_J1849+3024.ms.listobs'
listobs()

Check the listobs files and go on to the next part CASA_1848+283_J1849+3024. Congratulations this is the end of part 1. This is a scripted example of phase-referencing, imaging and self-calibration.