EVLA high frequency spectral line tutorial - IRC+10216 - calibration
This guide is under development.
Getting the data
The data for this tutorial can be obtained by anonymous FTP from ftp://ftp.aoc.nrao.edu/staff/gvanmoor/community_day/. Download all 4 files TAR in this directory.
# In UNIX wget 'ftp://ftp.aoc.nrao.edu/staff/gvanmoor/community_day/*'
Inspection and Flagging
List observational parameters
listobs provides almost all relevant observational parameters such as correlator setup (frequencies, bandwidths, channel number and widths, polarization products), sources, scans, scan intents, and antenna locations. The following command writes the listobs output to a file on disk. You can review the output at any time using less from the casapy or UNIX command line.
# In CASA vis='day2_TDEM0003_20s_full' listobs(vis=vis, verbose=True, listfile=vis+'.listobs.txt')
Here is a subset of the listobs output.
Fields: 4 ID Code Name RA Decl Epoch SrcId nVis 2 D J0954+1743 09:54:56.82363 +188.8.131.524 J2000 2 32726 3 NONE IRC+10216 09:47:57.38200 +184.108.40.20600 J2000 3 99540 5 F J1229+0203 12:29:06.69973 +02.03.08.5982 J2000 5 5436 7 E J1331+3030 13:31:08.28798 +220.127.116.1189 J2000 7 2736 (nVis = Total number of time/baseline visibilities per field) Spectral Windows: (2 unique spectral windows and 1 unique polarization setups) SpwID #Chans Frame Ch1(MHz) ChanWid(kHz) TotBW(kHz) Corrs 0 64 TOPO 36387.2295 125 8000 RR RL LR LL 1 64 TOPO 36304.542 125 8000 RR RL LR LL Sources: 10 ID Name SpwId RestFreq(MHz) SysVel(km/s) 0 J1008+0730 0 0.03639232 -0.026 0 J1008+0730 1 0.03639232 -0.026 2 J0954+1743 0 0.03639232 -0.026 2 J0954+1743 1 0.03639232 -0.026 3 IRC+10216 0 0.03639232 -0.026 3 IRC+10216 1 0.03639232 -0.026 5 J1229+0203 0 0.03639232 -0.026 5 J1229+0203 1 0.03639232 -0.026 7 J1331+3030 0 0.03639232 -0.026 7 J1331+3030 1 0.03639232 -0.026 Antennas: 19: ID Name Station Diam. Long. Lat. 0 ea01 W09 25.0 m -18.104.22.168 +22.214.171.124 1 ea02 E02 25.0 m -107.37.04.4 +33.54.01.1 2 ea03 E09 25.0 m -126.96.36.199 +188.8.131.52 3 ea04 W01 25.0 m -107.37.05.9 +33.54.00.5 4 ea05 W08 25.0 m -184.108.40.206 +220.127.116.11 5 ea07 N06 25.0 m -107.37.06.9 +18.104.22.168 6 ea08 N01 25.0 m -107.37.06.0 +33.54.01.8 7 ea09 E06 25.0 m -22.214.171.124 +126.96.36.199 8 ea12 E08 25.0 m -188.8.131.52 +184.108.40.206 9 ea15 W06 25.0 m -220.127.116.11 +18.104.22.168 10 ea19 W04 25.0 m -22.214.171.124 +126.96.36.199 11 ea20 N05 25.0 m -107.37.06.7 +33.54.08.0 12 ea21 E01 25.0 m -107.37.05.7 +188.8.131.52 13 ea22 N04 25.0 m -107.37.06.5 +33.54.06.1 14 ea23 E07 25.0 m -184.108.40.206 +220.127.116.11 15 ea24 W05 25.0 m -18.104.22.168 +22.214.171.124 16 ea25 N02 25.0 m -107.37.06.2 +33.54.03.5 17 ea27 E03 25.0 m -107.37.02.8 +33.54.00.5 18 ea28 N08 25.0 m -107.37.07.5 +126.96.36.199
Note that the Rest Frequency and Systemic Velocity are wrong in the listobs log by a factor 10^6 and 1000, respectively, given the quoted units (MHz) and (km/s). This was due to a temporary error in the EVLA Observing Tool that has subsequently been fixed. Because the sky frequencies are correct, and we set the rest frequency explicitly later in the deconvolution stage, this does not present a problem for the data reductions.
We summarize the observing strategy in this table.
|Gain calibrator||J0954+1743||field id = 2|
|Bandpass calibrator||J1229+0203||field id = 5|
|Flux calibrator||J1331+3030 (3C286)||field id = 7|
|Science target||IRC+10216||field id = 3|
Plot antenna positions
We can use task plotants to generate a plot of the antenna positions. This is useful for choosing a reference antenna.
# In CASA plotants(vis=vis,figfile=vis+'.plotants.png')
We will use ea02 as the reference antenna.
Next, let's look at the elevation as a function of time for all sources. We will use plotms to plot the data. It's not the case for these data, but if the elevation is very low (usually at start or end of track) you may want to flag. Also, how near in elevation your flux calibrator is to your target will impact your ultimate absolute flux calibration accuracy. Unfortunately, the target and flux calibrator are not particularly well-matched for this observation, as you can show by plotting the elevation for each source (each sources has a different colors). Thus we are strongly dependent on the opacity and gaincurve corrections to get the flux scale right for these data. (This is something to keep in mind when planning observations!)
# In CASA plotms(vis=vis, xaxis='time',yaxis='elevation',correlation='RR,LL', avgchannel='64',spw='0:4~60', coloraxis='field')
Next, let's look at all the source amplitudes as a function of time.
# In CASA plotms(vis=vis, xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60', coloraxis='field')
Now zoom in on the region very near zero amplitude for sources J0954+1743 and IRC+10216. To zoom, select the Zoom tool in lower left corner of the plotms GUI, then you can left click to draw a box. Look for the low values (you may want to zoom a few times to really see the suspect points clearly). Now use the Mark Region and Locate buttons (located along the bottom of the GUI) to see which antenna is causing problems. The output is be shown in the logger. Since all the "located" baselines include ea12, this is the responsible antenna.
Now click the clear region button, and then go back to the zoom button to zoom in further to note exactly what the time range is: 03:41:00~04:10:00.
Check the other sideband by changing spw to 1:4~60. You will have to rezoom. If you have trouble, click on the Mark icon and then back to zoom. In spw=1, ea07 is bad from the beginning until after next pointing run: 03:21:40~04:10:00. To see this, compare the amplitudes when antenna is set to 'ea07' and when it is set to one of the other antennas, such as 'ea08'.
If you set antenna to 'ea12' and zoom in on this intial timerange, you can also see that ea12 is bad during the same time range as for spw 0. You can also see this by entering '!ea07' for antenna, which removes ea07 from the plot (in CASA selection, ! deselects).
We can set up a flagging command to get both bad antennas for the appropriate time and spw:
# In CASA flagdata(vis=vis, field=['2,3','2,3'], spw=['','1'], antenna=['ea12','ea07'], timerange=['03:41:00~04:10:00','03:21:40~04:10:00'])
flagdata works by spanning up a matrix. The first entries in each list must be taken as one flagging command, as well as the second entries etc. Lists within lists are fine. In the above example, the first flagging command is issued for fields 2 and 3 for all spws and within the 03:41:00~04:10:00 timerange. A second command is again for the fields 2 and 3 but for spw 1 only and for the second timerange in the list '03:21:40~04:10:00'.
Note that because the chosen timerange is limited to fields 2 and 3, the field parameter is not really needed; however, flagdata will run fastest if you put as many constraints as possible.
Now remove the !ea07 from antenna and replot both spw, zooming in to be sure that all obviously low points are gone. Also zoom in and check 3C286 (J1229+0203 is already obvious because it is so bright!).
Lets look more closely at IRC+10216:
# In CASA plotms(vis=vis,field='3', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0~1:4~60', coloraxis='spw')
You can see a that there are some noisy high points. But now try
# In CASA plotms(vis=vis, field='3', xaxis='uvdist', yaxis='amp', correlation='RR,LL', avgchannel='64', spw='0~1:4~60', coloraxis='spw')
Most of the high points on IRC+10216 are due to large scale emission on short baselines, but there is still some noisy stuff -- for a target like this with extended emission it's best to wait until later to decide what to do about it. We will not be able to get adequate calibration for antennas that are truly bad (even if they don't stand out here) so these will be obvious later.
Set the model for the flux calibrator
Next, we set the model for the flux calibrator. Depending on your observing frequency and angular resolution you can do this several ways. In the past, one typically used a point source (constant flux) model for the flux calibrator, possibly with a uvrange cutoff if necessary. More recently for the VLA/EVLA, model images for the most common flux calibrators have been made available for use in cases where the sources are somewhat resolved. This is most likely to be true at higher frequencies and at higher resolutions (more extended arrays).
Currently, CASA contains models for the most common calibrators at the most common frequencies but not yet the new EVLA frequency bands S, Ku, and Ka. These will be added as soon as they become available. One may check the availability of calibration models in setjy:
# In CASA setjy(vis=vis,listmodimages=T)
The terminal will now show the available models, e.g. 3C286_C.im, 3C48_K.im etc. (setjy will search in the working directory for images that may contain models, as well as in a CASA directory where known calibrator models are stored)
As mentioned above, the Ka band does not yet have a model. We will use the K-band model instead.
# In CASA setjy(vis=vis, field='7', modimage='3C286_K.im')
The logger output for each SPW is:
J1331+3030 (fld ind 7) spw 0 [I=1.7762, Q=0, U=0, V=0] Jy, (Perley-Butler 2010) J1331+3030 (fld ind 7) spw 1 [I=1.7794, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)
The absolute fluxes for the frequencies have now been determined and we can proceed to the bandpass and complex gain calibrations.
Before determining the bandpass solution, we need to inspect phase and amplitude variations with time and frequency on the bandpass calibrator to decide how best to proceed. We limit the number of antennas to make the plot easier to see. We chose ea02 as it seems like a good candidate for the reference antenna.
# In CASA plotms(vis=vis, field='5', xaxis='channel', yaxis='phase', correlation='RR', avgtime='1e8', spw='0:4~60', antenna='ea02', antenna='ea02&ea23')
The phase variation is modest ~10 degrees. Now expand to all baselines that include ea02, then hit the "Plot" button.
# In CASA plotms(vis=vis,field='5', xaxis='channel', yaxis='phase', correlation='RR', avgtime='1e8', spw='0:4~60', antenna='ea02', coloraxis='antenna2')
From this you can see that the phase variation across the bandpass is modest. Next check LL, and spw=1, both correlations. Also check other antennas if you like.
Now look at the phase as a function of time.
# In CASA plotms(vis=vis, field='5', xaxis='time', yaxis='phase', correlation='RR', avgchannel='64', spw='0:4~60', antenna='ea02&ea23')
Expand to all antennas with ea02
# In CASA plotms(vis=vis, field='5', xaxis='time', yaxis='phase', correlation='RR', avgchannel='64', spw='0:4~60', antenna='ea02', coloraxis='antenna2')
You may want to select "Custom" under "unflagged points symbol" and then change Style from "2" to "3" under the "Display" tab, "Unflagged Points Symbol" .
You can see that the phase variations are smooth, but do vary significantly over the 5 minutes of observation -- in most cases by a few 10s of degrees. Zoom in to see this better if you want.
The conclusion from this investigation is that you need to correct the phase variations with time before solving for the bandpass to prevent decorrelation of the vector averaged bandpass solution. Since the phase variation as a function of channel is modest, you can average over several channels to increase the signal to noise of the phase vs. time solution. If the phase variation as a function of channel is larger you may need to use only a few channels to prevent introducing delay-based closure errors as can happen from averaging over non-bandpass corrected channels with large phase variations.
Since the bandpass calibrator is quite strong we do the phase-only solution on the integration time of 20 seconds (solint='int').
# In CASA bpphase = vis + '.bpphase.gcal' gaincal(vis=vis, caltable=bpphase, field='5', spw='0~1:20~40', solint='int', refant='ea02', gaintype='G', calmode='p', gaincurve=True)
Plot the solutions.
# In CASA plotcal(caltable=bpphase, xaxis='time',yaxis='phase', iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
These solutions will appear in the CASA plotter gui. If you closed it after plotting the antennas above, it should reopen. If it is still open from before, the new plots should just appear. After you are done looking at the first set of plots, push the "Next" button on the GUI to see the next set of antennas.
Next we can apply this phase solution on the fly while determining the bandpass solutions on the timescale of the bandpass calibrator scan (solint='inf' ).
# In CASA bpass = vis + '.bandpass.bcal' bandpass( vis=vis, caltable=bpass, field='5', solint='inf', refant='ea02', solnorm=True, gaintable=bpphase, gaincurve=True )
Plot the solutions, amplitude and phase.
# In CASA plotcal(caltable=bpass, xaxis='chan', yaxis='amp', iteration='antenna', subplot=331)
# In CASA plotcal(caltable=bpass, xaxis='chan', yaxis='phase', iteration='antenna', subplot=331)
Note that phases on ea12 look noiser than on other antennas. This jitter could indicate bad pointing; note that ea12 had just come back in the array.
Now that we have a bandpass solution to apply we can solve for the antenna-based phase and amplitude gain calibration. Since the phase changes on a much shorter timescale than the amplitude, we will solve for them separately. In particular, if the phase changes significantly over a scan time, the amplitude would be decorrelated, if the un-corrected phase were averaged over this timescale. Note that we re-solve for the gain solutions of the bandpass calibrator, so we can derive new solutions that are corrected for the bandpass shape. Since the bandpass calibrator will not be used again, this is not strictly necessary, but is useful to check its calibrated flux density for example. We use a minimum signal-to-noise of 2 here as it seems to be a good compromise for using good data without rejecting too many solutions (minsnr=2).
# In CASA gaincal(vis=vis, caltable=intphase, field='2,5,7', spw='0~1:4~60', refant='ea02', calmode='p', solint='int', minsnr=2.0, gaintable=bpass, gaincurve=T)
Here solint='int' coupled with calmode='p' will derive a single phase solution for each 20 second integration. Note that the bandpass table is applied on-the-fly before solving for the phase solutions, however the bandpass is NOT applied to the data permanently until applycal is run later on.
Note that quite a few solutions are rejected due to SNR<2 (printed to terminal). For the most part it is only one or two solutions out of >30 so this isn't too worrying. Take note if you see large numbers of rejected solutions per integration. This is likely an indication that solint is too short for the S/N of the data.
Now look at the phase solution, and note the obvious scatter within a scan time.
# In CASA plotcal(caltable=intphase, xaxis='time',yaxis='phase', iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
Although solint='int' (i.e. the integration time of 10 seconds) is the best choice to apply before for solving for the amplitude solutions, it is not a good idea to use this to apply to the target. This is because the phase-scatter within a scan can dominate the interpolation between calibrator scans. Instead, we also solve for the phase on the scan time, solint='inf' (but combine=' ', since we want one solution per scan) for application to the target later on. Unlike the bandpass task, for gaincal, the default of the combine parameter is combine=' '.
# In CASA scanphase = vis + '.scanphase.gcal' gaincal(vis=vis, caltable=scanphase, field='2,5,7',spw='0~1:4~60', refant='ea02',calmode='p',solint='inf',minsnr=2.0, gaincurve=T)
# In CASA plotcal(caltable=scanphase,xaxis='time',yaxis='phase', iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
Note that there are no failed solutions here because of the added S/N afforded by the longer solint. Alternatively, instead of making a separate phase solution for application to the target, one can also run smoothcal to smooth the solutions derived on the integration time.
Next we apply the bandpass and solint='int' phase-only calibration solutions on-the-fly to derive amplitude solutions. Here the use of solint='inf', but combine=' ' will result in one solution per scan interval.
# In CASA amp_gcal = vis + '.amp.gcal' gaincal(vis=vis, caltable=amp_gcal, field='2,5,7',spw='0~1:4~60', refant='ea02',calmode='ap',solint='inf',minsnr=2.0, gaintable=[bpass, intphase], gaincurve=T)