VLBA Basic Spectral Line Tutorial DRAFT
This CASA Guide is for Version 6.6.3 of CASA. If you are using a later version of CASA and this is the most recent available guide, then you should be able to use most, if not all, of this tutorial.
Overview
This CASA Guide describes the procedure for calibrating a phase-referenced VLBA observation of methanol masers in W3OH. The data were taken specifically for this tutorial. The observation made use of the DDC observing personality, using dual polarization with one spectral windows in each polarization. The spectral window has a bandwidth of 32 MHz divided into 256 spectral channels centered at 6668 MHz for the methanol maser line.
This tutorial will focus on calibrating the data and creating a continuum (Stokes I) image of the phase reference calibrator and images of the masers in the W3OH.
Those who have already worked through the VLBA Basic Phase-referencing tutorial will recognize that many of the calibration steps in this tutorial are very similar (sometimes identical). In general, the initial calibration (data ingestion through amplitude calibration) of spectral line data are the same. The major differences between continuum and spectral line projects are in the target imaging and specific spectral line analysis steps.
Please note that CASA should be used with caution when calibrating VLBA data. At the current time, CASA should only be used to calibrate relatively simple VLBA observations (basic continuum, simple phase-referenced, etc.). In particular, CASA is currently not recommended for calibrating the following types of VLBA observations:
- Polarimetric — Calibration of resolved polarized sources is not yet available in CASA.
- Astrometric — Earth orientation parameter (EOP) and pulse-cal tone corrections are not yet available in CASA.
- Spectral line projects requiring fringe-rate mapping — Fringe-rate mapping is not yet available in CASA.
- Low Frequency (<4 GHz) — Total electron content (TEC) corrections have not yet been verified to work correctly for VLBI observations in CASA. Also, note that many C-band (4 to 8 GHz) VLBA observations will benefit from TEC corrections.
- VLBA+Y27 or other VLBA+ arrays that require the use of antab files — CASA does not currently support the use of antab files for ingesting calibration data (system temperatures, gain curves, etc.)
- VLBA data correlated prior to December 2009 — CASA cannot load data correlated with the old VLBA hardware correlator. Your data must have been correlated with the DiFX correlator.
If your observation involves any of the above, you should use AIPS to calibrate your data.
A note on VLBA+Y1 observations: It is currently possible to calibrate VLBA+Y1 (VLBA and single VLA antenna) in CASA. Any data taken after 2022 July 06 should work with minimal extra steps (namely, checking that the Y1 antenna diameter is correct). VLBA+Y1 data observed prior to 2022 July 06 may not include the system temperature and/or gain curve values need for calibration. For these data sets, refer to VLBA Scientific Memo 39 for details on how to prepare the data for calibration in either CASA or AIPS.
How to Use This CASA Guide
There are a number of possible ways to run CASA, described in more detail in Getting Started in CASA. In brief, there are at least three different ways to use CASA:
- Interactively examining task inputs. In this mode, one types taskname to load the task, inp to examine the inputs, and go once those inputs have been set to your satisfaction. Allowed inputs are shown in blue and bad inputs are colored red. The input parameters themselves are changed one by one, e.g., selectdata=True. Screenshots of the inputs to various tasks used in the data reduction below are provided, to illustrate which parameters need to be set. More detailed help can be obtained on any task by typing help taskname. Once a task is run, the set of inputs are stored and can be retrieved via tget taskname; subsequent runs will overwrite the previous tget file.
- Pseudo-interactively via task function calls. In this case, all of the desired inputs to a task are provided at once on the CASA command line. This tutorial is made up of such calls, which were developed by looking at the inputs for each task and deciding what needed to be changed from default values. For task function calls, only parameters that you want to be different from their defaults need to be set.
- Non-interactively via a script. A series of task function calls can be combined together into a script, and run from within CASA via execfile('scriptname.py'). This and other CASA Tutorial Guides have been designed to be extracted into a script via the script extractor by using the method described at the Extracting_scripts_from_these_tutorials page. Should you use the script generated by the script extractor for this CASA Guide, be aware that it will require some small amount of interaction related to the plotting, occasionally suggesting that you close the graphics window and hitting return in the terminal to proceed. It is in fact unnecessary to close the graphics windows (it is suggested that you do so purely to keep your desktop uncluttered).
If you are a relative novice or just new to CASA, it is strongly recommended to work through this tutorial by cutting and pasting the task function calls provided below after you have read all the associated explanations. Work at your own pace, look at the inputs to the tasks to see what other options exist, and read the help files. Later, when you are more comfortable, you might try to extract the script, modify it for your purposes, and begin to reduce other data.
NOTE: It is not recommended to use scripts to generate "science-ready" VLBA data products. Nearly all VLBA observations will require some amount of hands-on calibration. Calibration and imaging scripts should only be used as "first pass" attempts at calibration, which can be useful to determine if the observation resulted in a detection and to identify problems in the data.
Obtaining the Data
This Guide is intended to cover the entire process one would follow for calibrating their own VLBA observation. Therefore, we will start with a FITS-IDI file rather than a Measurement Set. The FITS-IDI file for this Guide is: TL016B.idifits (right-click and select "Save Link As..."; file size: 1.8 GB). Additionally, you will need to download the zoom band file containing the high spectral resolution data for the masers.
Alternatively, you may download the data from the NRAO Archive. In the archive search inputs, enter "TL016" in the Project Code field. Expand the results for TL016B and look for the file VLBA_TL016B_tl016b_BIN0_SRC0_1_220307T210851.idifits. Select the file, click "Download", and follow the instructions. Once you have downloaded the FITS-IDI file from the archive, it will be helpful to change the filename to "TL016G.idifits".
The Observing Log
Before diving into the calibration, it is always a good idea to look over the observing log to check for notes from the operators that can inform us about potential issues with the data (missing stations, bad weather, etc.). These logs are always emailed to the PI's of an observation, but you can also access them later from the NRAO's vlbiobs fileserver. To locate an observing log, first find the directory for observing month and the last two digits of the year (for the observation used in this Guide, that directory is jan24 for January 2024). Once inside the proper month+year directory, look for the project code (in this case, tl016d). Look for a file named <project code>log.vbla (tl016dlog.vlba).
The observing log for this particular observation looked like this:
VERY LONG BASELINE ARRAY OBSERVING LOG
--------------------------------------
Project: TL016G
Observer: Linford, Justin
Project type: TEST
Obs filename: tl016d.vex
Date/Day: 2024JAN07/007
Ants Scheduled: SC HN NL FD LA PT KP OV BR MK
=UT-Time===Comment===============================================MF#===%AD==AMD=
Operator: Trevor Seddon
2359 Begin
2359-0416 *BR out, FRM stuck in SX position wo-19696 FRM 100 257
2359-0416 *KP out, elevation up limit and brake interlock fault WEA 100 257
that are likely caused by ice tripping the switches
2359-0112 *HN out to avoid snow WEA 100 73
0000 Date is 2024JAN08/008.
0030 %OV windy
0228 New Operator: Paul Padilla
0416 End.
Downtime Summary:
Total downtime : 587 min
Percentage downtime of observing: 22.8%
Average downtime per hour : 13.7 min
Total scheduled observing time (# Antennas): 2570 min (10)
Notes:
* = Entries where data was affected.
% = Entries where data may have been affected.
& = Entries where the site tech was called out.
WEA = Weather entries.
MF# = Maintenance form or major downtime category associated with a problem.
%AD = The percentage of an antenna affected by a problem.
AMD = Total antenna-minutes downtime for a problem.
Tsys = System Temperature (TP/SP x Tcal/2)
ACU = Antenna Control Unit
FRM = Focus/Rotation Mount
RFI = Radio Frequency Interference
CB = Circuit Breaker
Creating the Measurement Set
Before beginning our data reduction, we must start CASA. If you have not used CASA before, some helpful tips are available on the Getting Started in CASA page. Remember to start CASA in the directory containing the observation data.
Once you have CASA up and running, it is time to get the data into a format that CASA can use. Unlike VLA data, a VLBA observation is only available as a FITS-IDI file and cannot be downloaded as a CASA Measurement Set. So, the first step in calibrating a VLBA observation with CASA is to create a Measurement Set from the FITS-IDI file. To do this, we will use the task importfitsidi:
# In CASA
importfitsidi(fitsidifile='TL016G.idifits', vis='tl016g.ms', constobsid=True, scanreindexgap_s=6)
The scanreindexgap_s parameter is used to reconstruct scan boundaries in those cases where sources do not change between scans. In general, it is good to set scanreindexgap_s to some non-zero number to help CASA properly organize the scan list. The recommended value is 15 seconds, but shorter values may work as well (although you probably don't want to go much shorter than about 5 seconds). If you find that the resulting MS contains too few scans, run importfitsidi again with scanreindexgap_s set to a smaller number. If your MS has too many scans, especially multiple scans on the same source when you think it should just be one scan, run importfitsidi again with scanreindexgap_s set to a larger number.
NOTE: You may see warnings in the CASA window like "Table not yet implemented" and "Keyword has wrong data type". This is normal behavior for VLBA data and due to the fact that CASA does not yet handle the pulse-cal tone data in the FTIS-IDI file. Just ignore them. However, do not ignore any SEVERE warnings, as these can indicate there are serious problems with the data.
Inspecting the Data
Now that we have a Measurement Set, we need to look things over. This is where we will start to identify possible problems, but also find reliable antennas to use as the reference antenna.
The Observation Summary
It will be useful later to have the basic information about the observation. The task listobs will return a list of all the scans, the sources observed, which stations were used, and the frequency setup. It is possible to run listobs in two ways: printing information in the CASA logger, or saving the information to a file.
To simply display the information in the CASA logger:
#In CASA
listobs(vis='tl016g.ms')
You should see the listobs output in the CASA logger window:
OUTPUT TO FILL IN LATER
NOTE: You can also assign the listobs output to a python dictionary (e.g., "obs_dict") by typing "obs_dict = listobs(vis='tl016g.ms')".
It is usually useful to have a copy of the listobs output in a file that you can refer to later. To save the listobs output to a file named "tl016b_listobs.txt", run listobs again with listfile='tl016b_listobs.txt'.
#In CASA
listobs(vis='tl016g.ms', listfile='tl016g_listobs.txt')
A few things that are useful to note from the listobs output:
- While the antenna number, field number, and spw number all start at 0, the scan number starts at 1.
- The fringe finder and bandpass calibrator is J0102+5824, with 4 scans of about 4 minutes each.
- There is also a back-up fringe finder and bandpass calibrator, J0217+7349.
- The phase reference calibrator is J0306+6243. Each phase reference scan is about 1 minute long.
Identifying a Good Reference Antenna
VLBA calibration makes use of a reference antenna: one antenna that all other antennas can be referenced to when generating solutions. Usually, you will want to use one of the more central stations: FD, KP, LA, or PT.
To determine which is the best to use as the reference antenna, use plotms to plot phase vs frequency for a single scan on a bright source (usually the fringe finder or bandpass calibrator). For this observation, use the second scan on J0102+5824, which is scan ?? from the listobs output. Inspect each of the 4 central antennas and use the one which has the most well-behaved phases. We start with scan ?? because it is near the middle of our observation, so J0102+5824 should be at relatively high elevations at all stations during that scan.
#In CASA
plotms(vis='tl016g.ms', xaxis='frequency', yaxis='phase', field='J0102+5824', scan='??', correlation='rr,ll', antenna='FD,KP,LA,PT', iteraxis='baseline', coloraxis='corr')
Use the GUI to step though the baselines. Our observation does not include the cross-hand polarizations (lr and rl), and you usually won't worry about inspecting them anyway since they will usually be too low to see anything useful.
We should also look at our back-up fringe finder, J0217+7349, to see how it compares to J0102+5824. In the GUI, change the "Field" parameter to J0217+7349, and change the scan number to ?? (the second scan on J0217+7349).
In the GUI, go to the "Axes" tab and change the "Data" axis to "Amp" (this is equivalent to setting yaxis='amplitude') to inspect the shapes of the bandpasses and look for RFI. If an antenna has strong RFI, you should not use it as the reference antenna. Also look at the bandpass shapes for J0217+7349.
Should see that J0102 is the better calibration target
For spectral line observations, it is often a good idea to also inspect the autocorreltion data for signs of strong RFI at each site. You can plot the autocorrelations by setting the antenna parameter to *&&&.
#In CASA
plotms(vis='tl016g.ms', xaxis='frequency', yaxis='amp', field='J0102+5824', scan='??', correlation='rr,ll', antenna='*&&&', iteraxis='antenna', coloraxis='corr')
WILL WANT AN IMAGE HERE OF THE PT AUTOCORRELATIONS SHOWING THAT NASTY CRUD
For this observation, we will use Fort Davis (FD) as the reference antenna. Its phases are fairly smooth, it is free of RFI, and the bandpasses look pretty clean.
Identifying a Good Time Range for the Instrumental Delay Calibration
Next, we will look for a good time interval (about 1 minute long) for the instrumental delay calibration. This will be done using the fringe finder (J0102+5824), so we should inspect each scan on it: ??1, ??2, ??3, and ??4. We will want to plot amplitude vs time for each of these scans. To save time, we will just plot the baselines to the reference antenna, FD. Start by plotting scan 1.
#In CASA
plotms(vis='tl016b.ms', xaxis='frequency', yaxis='amp', field='J0102+5824', antenna='FD', scan='1', correlation='rr,ll', iteraxis='baseline', coloraxis='corr')
Use the GUI to step through each baseline for the scan. Then, use the "scan" field in the GUI to look at the other scans on J0102+5824 (scans ??, ??, ??, and ??). Look for outliers in the amplitude data. We want to use a time range with no RFI and with all antennas on source.
Scan ?? is close to the middle of the observation, and looks mostly clean. Use the plotms GUI to set the "scan" to "??" (in the "Data" tab) and the "X Axis" to "Time" (in the "Axes" tab) and step through the baselines again. Notice that some baselines have some higher amplitude points early on in this scan. We will want to avoid these early times. The time range TIME1 to TIME2 looks nice for each baseline to the reference antenna, so we will use that for our instrumental delay calibration later.
Looking for the Maser Lines
This is a spectral line dataset, so we should take a minute to check that we the maser lines on the science target.
Again, use plotms to plot the data.
#In CASA
plotms(vis='tl016b.ms', xaxis='frequency', yaxis='amp', field='W3OH', antenna='FD', scan='??', correlation='rr,ll', iteraxis='baseline', coloraxis='corr')
Stepping through the baselines, you should see that there is indeed a large amplitude spike near the expected maser frequency of 6668.519 MHz on all baselines.
Flagging Data
Should not need to quack these data
Flagging Problem Antennas
We know from our inspection of the fringe finders that PT has some major RFI near the maser line. Therefore, we should not include PT in our calibration. In fact, we should not use PT at all since it may be contributing bad data and/or noise to our observations.
Flag the PT antenna for the entire observation by using the flagdata command
#In CASA
flagdata(vis='tl016g.ms', mode='manual', antenna='PT')
Automated Flagging
The the flagdata task includes a useful automated flagging mode known as "TFCrop". For more details on using the automated TFcrop mode, see the VLA CASA Flagging topical tutorial. Because this is a spectral line data set, we want to limit the automatic flagging to only the time domain.
To begin, you just want to get a feel for what TFCrop will do to the data, which means you should set action='calculate' .
#In CASA
flagdata(vis='tl016g.ms', mode='tfcrop', datacolumn='data', timecutoff=4.0, flagdimension='time', action='calculate', display='both')
This will open a new window with a GUI. The top row of the GUI will display the data as it currently is. The bottom row shows what the data will look like after the flags are applied (blue regions are flagged). Step through the antennas by clicking "Next Baseline". Step through the spectral windows by clicking "Next SPW". Step through at least a few scans by clicking "Next Scan". With timecutoff=4.0, it should really only flag those data points that are 4-times the standard deviation in the time dimension. This should be sufficient to remove the worst of the time-dependent RFI in our observation.
To generate the flags, we will need to run flagdata again with action='apply' .
#In CASA
flagdata(vis='tl016g.ms', mode='tfcrop', datacolumn='data', timecutoff=4.0, flagdimension='time', action='apply', display='both')
You should double check that the automatic flagging will behave as you expect. Once you are satisfied that the flagging will do what you want, click "Stop Display". The program will then look through the entire observation and apply the flags.
NOTE: You can always exit the TFCrop mode without applying any flags by clicking "Quit", even when you set action='apply' .
Flagging "By Hand"
As telescope arrays become larger and more complex, automated flagging will become the primary means of excluding bad data (just look at what ALMA does). However, the current VLBA is still small enough that inspecting and flagging data oneself is not too painful.
You can look for obvious bad data using plotms. Be aware that while you can flag data in plotms, a major drawback to using plotms to flag data is that it makes undoing flags very difficult. The preferred method of flagging by hand is to identify the bad data with plotms, and then generate the flags with flagdata. It is recommended to be somewhat careful with flagging early on in the calibration process and only flag those data points that are obviously bad.
To start with, let's plot the bandpasses for each baseline to FD for all 5 scans on 4C39.25.
#In CASA
plotms(vis='tl016g.ms', xaxis='frequency', yaxis='amp', field='4C39.25', antenna='FD', scan='', correlation='rr,ll', iteraxis='baseline', coloraxis='scan')
Calibrating the Data
Amplitude Corrections from Autocorrelations
Determine the amplitude corrections from the autocorrelations with accor.
#In CASA
accor(vis='tl016g.ms', caltable='tl016g.accor', solint='30s')
NOTE: This step is not required for EVN data, because the EVN correlator performs it during correlation.
Inspect the tl016g.accor solution table with plotms.
#In CASA
plotms(vis='tl016g.accor', xaxis='time', yaxis='amp', iteraxis='antenna', coloraxis='spw')
Ideally, all of the solutions should be very close to 1.0. Stepping through the antennas, it is clear that pretty much every station has significant outliers.
The AIPS VLBA utility script VLBACCOR smooths the autocorrelation corrections by default (with a smoothing time of 30 minutes). It is possible to do this smoothing in CASA 6.3 and later using the smoothcal task. Based on what we saw in our plotting, these data would likely benefit from some smoothing.
#In CASA
smoothcal(vis='tl016g.ms', tablein='tl016bgaccor', caltable='tl016g_smooth.accor', smoothtype='median', smoothtime=1800.0)
Remember to check the smoothed solutions with plotms to make sure it was an improvement.
#In CASA
plotms(vis='tl016g_smooth.accor', xaxis='time', yaxis='amp', iteraxis='antenna', coloraxis='spw')
A Priori Calibration
Unlike the VLA, the VLBA cannot rely on bootstrapping absolute flux density calibration from a well-modeled calibrator. Instead, the VLBA relies on a combination of the system temperature and the known gain curve of the antennas (how the gain of the antenna changes with elevation). Both the system temperatures and gain curves are included in the FITS-IDI file, and both are imported into the Measurement Set with importfitsidi. To get the system temperature and gain curve information into a form that CASA can use for calibration, we use the gencal task to generate calibration tables.
System temperature:
gencal(vis='tl016g.ms', caltable='tl016g.tsys', caltype='tsys', uniform=False)
You will see several messages appear informing you that the Tsys values for many times were negative and will be flagged.
Check the system temperature table with plotms.
#In CASA
plotms(vis='tl016g.tsys', xaxis='time', yaxis='tsys', iteraxis='antenna', coloraxis='spw')
Make sure to also plot the solutions with xaxis='freq' .
DO THE MASERS SHOW UP IN TSYS VS FREQ?
Gain curve:
#In CASA
gencal(vis='tl016g.ms', caltable='tl016g.gcal', caltype='gc')
This observation was done at 6.7 GHz and with a bandwidth of only 32 MHz (??). The gain curve is very uninteresting, so we will not bother looking at the solutions now (although you can if you really want to). If your observation is at 12 GHz or higher, you will probably want to inspect the gain curve table with plotms.
NOTE: For EVN observations, the system temperatures and gain curves are provided in ANTAB files rather than in the FITS-IDI file. If you are working with EVN data, you must append the system temperature data to the FITS-IDI file before running importfitsidi. EVN gain curves must also be extracted from the ANTAB file. For more information, see the EVN Data Reduction Guide.
Instrumental Delay Calibration
Solve for the instrumental delays by using fringefit on a bright source (the "fringe finder"). In our case, the fringe finder is ????. Set the timerange to the one-minute time span you identified while inspecting the data earlier. Because we are only solving for the instrumental delays (changes in phase as a function of frequency), we will set zerorates=True to set all of the delay-rates (changes in phase as a function of time) to be zero.
#In CASA
fringefit(vis='tl016g.ms', caltable='tl016g.sbd', field='???', timerange='??:??:??~??:??:??', solint='inf', zerorates=True, refant='FD', minsnr=10, gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys'], interp=['nearest', 'nearest', 'nearest,nearest'], parang=True)
Watch the logger for any reports of low SNR and failures to converge on a solution. Ideally, the SNR for each station should be very high (for our data, all reported SNR values are larger than 3900, which is very good).
A Note on Table Names: In the Basic Phase-referencing tutorial, we followed the convention of naming tables as described in the EVN CASA tutorials. The instrumental delay table was called "sbd", which stood for "single band delay". It is not necessary to follow this convention. You can name a calibration table anything you want. For these data, perhaps we should use caltable='tl106g.instdly' for the instrumental delay. Just be sure to keep track of your tables and to apply them appropriately in the calibration tasks.
Once fringefit has finished, you should see the following in the logger:
Finished solving. Calibration solve statistics per spw: (expected/attempted/succeeded): Spw 0: 1/1/1
This is exactly what we expect to see. We are only attempting to get one solution per spw because we are only solving over a single solution interval (one minute of data from a single scan).
Apply the instrumental delay corrections using applycal.
#In CASA
applycal(vis='tl016g.ms', gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest'], parang=True)
NOTE: In applycal, the order in which you specify the list for gaintable sets the order for both interp and spwmap. If you change the order of gaintable, be sure to also change the order of interp and spwmap! Also, if you smoothed any of the solutions, be sure to use the appropriate filename for the smoothed table.
It may take a little while (~2 minutes or more) for applycal to finish. If you are watching the logger carefully, you will see the message "Adding CORRECTED_DATA column(s)". Because this is the first time that we have run applycal, CASA needs to create a new column containing the data with the corrections applied.
Check that applying the solutions resulted in improvements:
#In CASA
plotms(vis='tl016g.ms', field='????', xaxis='frequency', yaxis='phase', ydatacolumn='corrected', timerange='??:??:??~??:??:??', correlation='rr,ll', antenna='*&*', iteraxis='antenna', coloraxis='antenna2')
All of the phase values should now be clustered around 0 degrees. There may be some outliers near the band edges, especially on MK. This is nothing to worry about, yet.
Global Fringe Fitting
Now, we will solve for the time and frequency-dependent effects in phase using fringefit, also known as "global fringe fitting" (or sometimes "multi-band delay", if you combine all of the bands). We will need to pick a solution interval that is appropriate for our data. It should be at least 10 seconds, and no longer than the scan length on the phase reference calibrator. For this observation, we will use 30 seconds which will give us one solution per scan on the phase reference calibrator, and 7 solutions per scan on the fringe finder.
For refant, enter a list of antennas to try as the reference antenna. The preferred refant should be listed first, followed by the second choice, then third choice, and so on. It is not recommended to include Mauna Kea (MK) or Saint Croix (SC) in the list, unless the phase reference calibrator is very bright on the longest baselines. For our observation, BR had some serious issues in two of our four spectral windows, so we will omit it from our list of possible reference antennas, too. Set field to the fringe finder and phase reference calibrator.
#In CASA
fringefit(vis='tl016g.ms', caltable='tl016g.mbd', field='????, J????', solint='30s', minsnr=5, zerorates=False, refant='FD,PT,LA,KP,OV,NL,HN', gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest'], parang=True)
This step may take quite a while (probably at least 20 minutes), so this is a good opportunity to stand up, stretch, and go get a tasty beverage.
When the fringefit task is done, check the logger for the solution statistics ("expected/attempted/succeeded"). You want all three numbers to be the same, or at least very similar.
For our project, you should see something like this in the logger: For our project, you should see something like this in the logger:
Finished solving. Calibration solve statistics per spw: (expected/attempted/succeeded): Spw 0: 263/256/256
Because fringefit ran successfully and we are satisfied that the number of solutions is appropriate, we should now take a look at the solutions with plotms.
#In CASA
plotms(vis='tl016g.mbd', xaxis='time', yaxis='phase', iteraxis='antenna', coloraxis='spw')
Use the GUI to switch the yaxis between 'phase', 'delay', and 'delayrate'. The phase and delay solutions should both vary smoothly with time. If they do not, you may need to smooth the table before applying it. The delay-rates should be clustered around zero with some scatter. You should see that the solutions on all antennas look very good.
Now that we are confident that the global/multi-band solutions are good, we will apply them with applycal.
#In CASA
applycal(vis='tl016g.ms', gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd', 'tl016g.mbd'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest', 'linear'], parang=True)
It will probably take a while (about 3 minutes or longer) for applycal to run this time.
Take a look at the calibrated data with plotms to make sure the corrections are improving the phases.
#In CASA
plotms(vis='tl016g.ms', xaxis='frequency', yaxis='phase', ydatacolumn='data', field='????', antenna='*&*', scan=??', correlation='rr,ll', iteraxis='antenna', coloraxis='corr')
Use the GUI to step through the antennas. The uncalibrated data is pretty much a mess. Now, go to the "Axes" tab and change the Data Column from "data" to "corrected". The corrected data should look much better, with the phases centered around zero degrees with some scatter.
Bandpass Calibration
Now we will correct for the shape of the bandpass using the bandpass task. This step requires a very bright source (perferably >1 Jy on all baselines), so we will use our fringe finder 4C39.25. However, unlike when we solved for the instrumental delay (above), we will use all of the scans on the source.
#In CASA
bandpass(vis='tl016g.ms', caltable='tl016g.bpass', field='????', solint='inf', refant='FD', solnorm=True, bandtype='B', gaintable=['tl016b_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd', 'tl016g.mbd'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest', 'linear'], parang=True)
Inspect the solutions with plotms.
#In CASA
plotms(vis='tl016g.bpass', xaxis='frequency', yaxis='amp', iteraxis='antenna', coloraxis='corr')
For spectral line obsevrations, the bandpass solutions need to be very smooth. If you notice any outliers, you delete that table, go back to the data and look carefully for any outliers that may have led to poor bandpass solutions and flag them. Then run the bandpass calibration again and check that the solutions look better. While looking over the solutions, also make sure to check the phases (uses the GUI to set yaxis='phase'). You should see that the phase solutions are all near zero.
Apply the solutions with applycal.
#In CASA
applycal(vis='tl016g.ms', gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd', 'tl016g.mbd', 'tl016g.bpass'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest', 'linear', 'linear,linear'], parang=True)
It will probably take a while for applycal to run again, since you are also applying the solutions from the global fringe fitting.
After applying the bandpass solutions, take a look at the calibrated data with plotms.
#In CASA
plotms(vis='tl016g.ms', xaxis='frequency', yaxis='amp', ydatacolumn='corrected', field='????', antenna='*&*', correlation='rr,ll', iteraxis='baseline', coloraxis='corr')
In the Basic Phase-referincing tutorial, the next step was to do a final amplitude calirbation. However, that step is only necessary for data with bandwidths of 256 MHz or larger. That definitely is not the case for this project, so we will skip that step. (It should not hurt to do it, it just is not necessary.)
Also, in the Basic Phase-referencing tutorial, we split the calibrated data at this point. For this tutorial, we will continue with self-calibration using the orignal data MS. The reason for this will become obvious soon. (Students are encouraged to think about why that might be and see if they are correct later!)
Self-Calibration of Phase Reference Calibrator
Despite our best efforts with the initial calibration, VLBA observations will almost always require additional steps to improve the calibration. Because the VLBA antennas are separated by such large distances, the phases may not be perfectly calibrated by fringe fitting alone. This is partly because each antenna may be looking through vastly different troposphere/ionosphere, and partly because the phase reference calibrators are rarely point sources on long baselines (and fringe-fitting usually assumes the calibrators are point sources on all baselines).
We will improve the calibration using "self-calibration". The process of self-calibration involves building models of the sources, refining the calibration based on those models, and then repeating until you converge to a good solution. The models are created during the imaging process. Each iteration of self-calibration should produce images with lower noise.
More details on self-calibration can be found in the Basic Phase-referencing tutorial. This tutorial will only cover the steps.
Imaging the Calibrator
To create images of our sources, we will use tclean.
#In CASA
tclean(vis='tl016g.ms', field='J????', imagename='J????_sc1', imsize=[640], cell=['0.2mas'], stokes='I', deconvolver='clark', weighting='natural', niter=1000, interactive=True, savemodel='modelcolumn')
Setting interactive=True will allow us to define the areas where there is real flux from the source ("clean windows"). The value for niter needs to be large enough that we can go through enough clean iterations to recover most of the source flux. You can set niter to any arbitrarily large number (a few hundred should be plenty), but 1000 is an easy value to remember. Setting savemodel='modelcolumn' will save the model we create for the image to the model column of the measurement set so we can use it for self-calibration.
Running tclean with interactive=True will open a new Viewer Display Panel where you can control the cleaning of the image. You should see that J1154+6022 is fairly point-like. Use the zoom tool to select the region around the source. Use the "polygon drawing" tool to draw a clean window around the source, excluding the stuff to the left and right of the bright source. Click on the "Continue deconvolution" button (the green circular arrow on the right). Check the residuals when the display updates. Stop cleaning when the area around the source looks noise-like. For example, if you notice positive and negative regions with similar magnitudes (e.g., +0.003 and -0.003) inside your clean window, that is probably a good time to stop cleaning. To stop the cleaning process, click the "stop deconvolving now" button (the red circle with the white "X"). It may take tclean a while to finish and close.
NOTE: If you need to delete a clean window, select "Erase" at the far left of the green region in the GUI and draw a box around the clean window(s) you want to delete, then double click. Remember to select "Add" again to define new clean windows.
Tracking Improvement
To ensure that self-calibration is actually helping to improve the quality of the images, you should track some image statistics. The most reliable way to do this is to define 2 boxes in the image; one box that contains the source ("on-source"), and one box that does not contain any part of the source ("off-source"). You can use imview to inspect the image and determine where to place the on-source and off-source boxes.
#In CASA
imview
In the GUI, go to Data->Open->J1154_sc1.image, then click "raster image" to load the image we just created. Use the zoom tool to select the area around the source. Use the "rectangle drawing" tool to create a box around the source and record the pixel locations of the lower left and upper right corners. The on-source box should be fairly small and fit tightly around the source. For example, lower left = [305,303], upper right = [335,340] looks like it should work. Now, zoom out and draw another box somewhere off the source. You can make the off-source box pretty big. For example, lower left = [71,409], upper right = [588,604]. To make things a little easier, we will assign the values of our on- and off-source boxes to some variables in CASA:
#In CASA
j1154_onbox='305,303,335,340'
j1154_offbox='71,409,588,604'
Once you have your on- and off-source boxes defined, use the task imstat to gather some statistics about those regions.
#In CASA
imstat(imagename='J1154_sc1.image', box=j1154_onbox)
imstat(imagename='J1154_sc1.image', box=j1154_offbox)
Running imstat will return several values in both the logger and the CASA terminal. The values of interest for the self-calibration process are the on-source peak value ("max"), the on-source total flux density ("flux"), and the off-source RMS value ("rms"). It is a very good idea to record each of these values for each image (including the first image without any self-calibration applied) as you proceed with self-calibration. Ideally, the on-source peak and total flux values will increase slightly and the off-source RMS will decrease significantly as you improve the calibration. However, be very suspicious of large increases (>10%) in flux density values, especially when the off-source RMS does not improve dramatically.
For the first image of J1154+6022, the values should be fairly close to:
On-source: 'flux': 0.195 'max': 0.190 Off-source: 'rms': 0.00068
Record these values in your notes so you have a way to check if self-calibration is leading to improvements in the images.
Phase Self-Calibration
Now that we have an initial model of the phase reference calibrator (from the image), we can begin the self-calibration process. In CASA, this is done with the gaincal task.
Next, we will refine the phases. For this step, we will set the solution interval to 20 seconds (solint='20s' ) to give us two solutions per scan on the phase reference calibrator.
#In CASA
gaincal(vis='tl016g.ms', field='J1154+6022', caltable='tl016g.pcal', solint='20s', refant='FD', minblperant=3, gaintype='G', calmode='p', gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd', 'tl016g.mbd', 'tl016g.bpass'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest', 'linear', 'linear,linear'], parang=True)
Note that we are pye-applying all of the initial calibration when we do the phase self-calibration.
Checking the logger, you should see:
Finished solving. Calibration solve statistics per spw: (expected/attempted/succeeded): Spw 0: 266/238/238
Inspect the phase solution table with plotms.
#In CASA
plotms(vis='tl016g.pcal', xaxis='time', yaxis='phase', iteraxis='antenna', coloraxis='spw')
As with the delays, all of the solutions should be clustered around zero.
Notice that when refining both the delays and the phases we set minblperant=3 because closure phases require at least 3 baselines.
Apply the phase self-calibration solutions to the phase reference calibrator.
#In CASA
applycal(vis='tl016g.ms', field='J????', gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd', 'tl016g.mbd', 'tl016g.bpass', 'tl016g.pcal'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest', 'linear', 'linear,linear', 'linear'], parang=True)
Make a new image after the improved phase calibration using tclean just as we did previously. Remember to use a new imagename for this new image of the phase self-calibrated data.
#In CASA
tclean(vis='tl016g.ms', field='J????', imagename='J????_sc2', imsize=[640], cell=['0.2mas'], stokes='I', deconvolver='clark', weighting='natural', niter=1000, interactive=True, savemodel='modelcolumn')
Once you have created the new image, check for improvement with imstat. As long as you did not change the imsize or cell parameters, you can use the same on- and off-source boxes that we defined earlier.
#In CASA
imstat(imagename='J1154_sc2.image', box=j1154_onbox)
imstat(imagename='J1154_sc2.image', box=j1154_offbox)
For the second image of J1154+6022, the values should be fairly close to:
On-source: 'flux': 0.195 'max': 0.191 Off-source: 'rms': 0.00068
Record these values in your notes.
Amplitude Self-Calibration
#In CASA
gaincal(vis='tl016.ms', field='J????', caltable='tl016g.apcal', solint='inf', refant='FD', minblperant=4, gaintype='G', calmode='ap', solnorm=True, gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd', 'tl016g.mbd', 'tl016g.bpass','tl016g.pcal'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest', 'linear', 'linear,linear', 'linear'], parang=True)
Notice that we have set minblperant=4 because closure amplitudes require 4 baselines, and solint='inf' to use the full scan for each solution.
In the logger, you should see:
Finished solving. Calibration solve statistics per spw: (expected/attempted/succeeded): Spw 0: 114/114/114
All of the numbers match, which is what we like to see.
We should take a look at the solutions on each telescope using plotms.
#In CASA
plotms(vis='tl016g.apcal', xaxis='time', yaxis='amp', iteraxis='antenna', coloraxis='spw')
Apply the amplitude solutions to the phase reference calibrator.
#In CASA
applycal(vis='tl016g.ms', field='J????', gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd', 'tl016g.mbd', 'tl016g.bpass','tl016g.pcal','tl016g.apcal'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest', 'linear', 'linear,linear', 'linear', 'linear'], parang=True)
Make one last image of the phase reference calibrator, just to check for improvements. Since we will not be doing any more self-calibration on this source, we will set savemodel='none' to save some time and disk space. Remember to use a new imagename again.
#In CASA
tclean(vis='tl016g.ms', field='J????', imagename='J????_sc3', imsize=[640], cell=['0.2mas'], stokes='I', deconvolver='clark', weighting='natural', niter=1000, interactive=True, savemodel='none')
Check the image statistics to make sure the self-calibration is making an improvement
#In CASA
imstat(imagename='J????_sc3.image', box=j1154_onbox)
imstat(imagename='J????_sc3.image', box=j1154_offbox)
For the final image of J1154+6022, the values should be fairly close to:
On-source: 'flux': 0.215 'max': 0.209 Off-source: 'rms': 0.00016
This is a significant improvement over the phase self-calibrated image!
We are now done with the self-calibration of the data.
The Zoom Band Data
Next, we will load in the data where we have zoomed in on the maser lines.
# In CASA
importfitsidi(fitsidifile='TL016G_zoom.idifits', vis='tl016g_zoom.ms', constobsid=True, scanreindexgap_s=6)
These data have much higher spectral resolution than the full data set. Now that we have a MS, we will apply all of the solutions we drived from the full data set to the zoom band data.
#In CASA
applycal(vis='tl016g_zoom.ms', field='', gaintable=['tl016g_smooth.accor', 'tl016g.gcal', 'tl016g.tsys', 'tl016g.sbd', 'tl016g.mbd', 'tl016g.bpass','tl016g.pcal','tl016g.apcal'], interp=['nearest', 'nearest', 'nearest,nearest', 'nearest', 'linear', 'linear,linear', 'linear', 'linear'], parang=True)
Plot some of the data to make sure the calibration applied properly.
Now that we are satisfied the calibration is good, we can start making images of the maser lines. There are a few differences between making these images and when we made continuum images of the phase reference calibrator. Namely, we will use something something something cube mode.
There are actually many masers in our field of view. To make one image with all of them, you will need a large amount of memory (RAM). This should only be attempted on an NRAO cluster node, or a computer with at least 32 GB of RAM (preferably 64 GB).
Instead of making one large image, we will make a few small images of maser regions.
