VLA Self-calibration Tutorial-CASA5.7.0: Difference between revisions
| Line 156: | Line 156: | ||
EM: I think it would be really useful to explain in some detail when self-cal can help and when it can't as we have done below. But also why did we know that self-cal would work in this case. | EM: I think it would be really useful to explain in some detail when self-cal can help and when it can't as we have done below. But also why did we know that self-cal would work in this case. | ||
There are instances in which self-cal can improve your image quality and | There are instances in which self-cal can improve your image quality and others when it will not. | ||
A couple typical cases in which self-cal can help improve the image of the target source: | A couple typical cases in which self-cal can help improve the image of the target source: | ||
* Extensive artifacts from the source of interest due to calibration errors | * Extensive artifacts from the source of interest due to calibration errors | ||
| Line 168: | Line 169: | ||
EM: could we expound upon how one could tell the difference between the cases in which self-cal WILL help and bullets 1,2, and 4 when self-cal will not help? What do those actually look like? People who are just starting may not know how to tell the difference between these different types of errors. To be honest, I am not sure I could tell you what the last 4 bullets looks like. | EM: could we expound upon how one could tell the difference between the cases in which self-cal WILL help and bullets 1,2, and 4 when self-cal will not help? What do those actually look like? People who are just starting may not know how to tell the difference between these different types of errors. To be honest, I am not sure I could tell you what the last 4 bullets looks like. | ||
In this guide, | In the case of this guide, we knew that self-cal would improve the image quality because after an initial clean of the data, there were extensive artifacts from the source of interest due to calibration errors (). These errors manifested as. | ||
== The Initial Model == | == The Initial Model == | ||
Revision as of 09:14, 11 February 2020
This page is currently under construction.
Introduction
After the initial calibration of a dataset using the observed calibrator sources, there are likely to be residual phase and/or amplitude errors in the data of the target source. Self-calibration (self-cal) is the process of using a model of the target source to reduce the phase and amplitude errors in the visibilities of the same target source. This CASA guide describes the process of self-calibration and how to choose parameters to achieve the best result.
Fundamentally, self-calibration is almost identical to standard calibration. Both standard calibration and self-cal work by comparing the visibility data with model visibilities to solve for calibration solutions. With standard calibration, we are usually provided a model of our calibrator source by the observatory (e.g., VLA Flux-density calibrators) or we adopt a simple model (e.g., a 1 Jy point source at the phase center is a common assumption for VLA phase calibrators). With self-calibration we need to create a model of our target source, e.g., by imaging the target visibilities as we did in the above section. Then for both standard calibration and self-cal we solve for calibration solutions after making choices about the solution interval, signal-to-noise, etc. When applying the standard calibration solutions we use interpolation to correct the target, but for self-cal we apply the calibration solutions directly to the target field from which they were derived.
In this guide, we will create a model using the target data itself (using tclean), use this model to solve for and apply calibration solutions (using gaincal and applycal), then iteratively improve this model with further rounds of self-cal. This is the most common procedure, but there are other variants that are outside the scope of this guide. For example, your initial model for the target may come from fitting a model to the visibilities instead of imaging, or may be based on a priori' knowledge of the target field. Some applications of self-cal may use other calibration tasks, e.g., bandpass, instead of or in addition to gaincal.
Each "round" of self-calibration will follow the same general procedure:
- Conservatively clean the target and save the model.
- Use gaincal with a set of parameters to calculate a calibration table.
- Inspect the calibration solutions using plotms.
- Determine if the solutions should be applied, re-calculated with different parameters, or if it is time to stop.
- Use applycal to apply the table of solutions to the data.
- Use split to write out the calibrated data with the applied solutions, the starting point for the next round of self-cal.
- Start the next round of self-cal.
Note: There is an alternative and equally valid variant of the self-calibration procedure that is outlined in Appendix ... but not used in this guide.
Overall, this guide will cover:
- How to run the general self-calibration process
- How to determine optimal self-cal parameters (e.g. solution interval)
- When to stop the self-cal process
- Advanced topics related to self-cal (e.g., peeling)
The dataset in this guide is a VLA observation of a massive galaxy cluster, MOO J1506+5137, at z=1.09 and is part of the Massive and Distant Clusters of WISE Survey (MaDCoWS: Gonzalez et al. 2019). MOO J1506+5137 stands out in the MaDCoWS sample due to its high radio activity. From the 1300 highest significance MaDCoWS clusters in the FIRST footprint, a sample of 51 clusters with extended radio sources defined as having at least one FIRST source with a deconvolved size exceeding 6.5" within 1' of the cluster center was identified. This sample was observed with the VLA (PI: Gonzalez, 16B-289, 17B-197; PI: Moravec, 18A-039) as a part of a larger study (Moravec et al. 2020). Through these follow-up observations, it was discovered that MOO J1506+5137 is the cluster in the sample with the most radio sources within ~1' of the cluster center. It is unique in having five radio sources, with three of the sources having extended emission. The VLA dataset showcased in this tutorial combined with other datasets suggest that the radio activity among the massive galaxy population appears to be linked to the dynamical state of the cluster (Moravec et al. in prep).
Data for this Tutorial
Obtaining the Data
The original observation has been calibrated using the VLA CASA Pipeline. It can be downloaded from the archive by searching for the following SDM name: 17B-197.sb34290063.eb34589992.58039.86119096065.
We will instruct the archive to apply the calibration solutions derived by the original pipeline execution (EM: is this done automatically or do the readers need to fill in a field?). Therefore, the measurement set (MS) we will be downloading will contain both the raw (uncalibrated) visibilities and the calibrated visibilities, which will appear in the 'DATA' and 'CORRECTED_DATA' columns of the MS, respectively. The raw data alone is 11.4 GB, and this will grow to 20.5 GB after applying the calibration.
As an alternative to requesting these data from the archive you may instead download the calibrated MS directly here:
-- link to data download --
A smaller data set is also available (see Section Splitting the Target Visibilities)
Observation Details
First, we will start CASA in the directory containing the data and then collect some basic information about the observation. The task listobs can be used to display the individual scans comprising the observation, the frequency setup, source list, and antenna locations.
# in CASA
listobs(vis='17B-197.sb34290063.eb34589992.58039.86119096065.ms')
A portion of the listobs output appears below:
================================================================================
MeasurementSet Name: 17B-197.sb34290063.eb34589992.58039.86119096065.ms MS Version 2
================================================================================
Observer: Prof. Anthony H. Gonzalez Project: uid://evla/pdb/34052589
Observation: EVLA
Data records: 5290272 Total elapsed time = 2853 seconds
Observed from 13-Oct-2017/20:40:09.0 to 13-Oct-2017/21:27:42.0 (UTC)
ObservationID = 0 ArrayID = 0
Date Timerange (UTC) Scan FldId FieldName nRows SpwIds Average Interval(s) ScanIntent
13-Oct-2017/20:40:09.0 - 20:45:03.0 1 0 J1549+5038 550368 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [SYSTEM_CONFIGURATION#UNSPECIFIED]
20:45:06.0 - 20:50:03.0 2 0 J1549+5038 555984 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_PHASE#UNSPECIFIED]
20:50:06.0 - 20:59:27.0 3 1 MOO_1506+5136 1050192 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [OBSERVE_TARGET#UNSPECIFIED]
20:59:30.0 - 21:00:51.0 4 0 J1549+5038 151632 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_PHASE#UNSPECIFIED]
21:00:54.0 - 21:10:15.0 5 1 MOO_1506+5136 1050192 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [OBSERVE_TARGET#UNSPECIFIED]
21:10:18.0 - 21:11:39.0 6 0 J1549+5038 151632 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_PHASE#UNSPECIFIED]
21:11:42.0 - 21:21:03.0 7 1 MOO_1506+5136 1050192 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [OBSERVE_TARGET#UNSPECIFIED]
21:21:06.0 - 21:22:27.0 8 0 J1549+5038 151632 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [CALIBRATE_AMPLI#UNSPECIFIED,CALIBRATE_PHASE#UNSPECIFIED]
21:22:30.0 - 21:27:03.0 9 2 3C286 511056 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [CALIBRATE_BANDPASS#UNSPECIFIED,CALIBRATE_FLUX#UNSPECIFIED]
21:27:06.0 - 21:27:42.0 10 2 3C286 67392 [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3] [CALIBRATE_BANDPASS#UNSPECIFIED,CALIBRATE_FLUX#UNSPECIFIED]
(nRows = Total number of rows per scan)
Fields: 3
ID Code Name RA Decl Epoch SrcId nRows
0 NONE J1549+5038 15:49:17.468534 +50.38.05.78820 J2000 0 1561248
1 NONE MOO_1506+5136 15:06:20.353700 +51.36.53.63460 J2000 1 3150576
2 NONE 3C286 13:31:08.287984 +30.30.32.95886 J2000 2 578448
Spectral Windows: (16 unique spectral windows and 1 unique polarization setups)
SpwID Name #Chans Frame Ch0(MHz) ChanWid(kHz) TotBW(kHz) CtrFreq(MHz) BBC Num Corrs
0 EVLA_C#A0C0#0 64 TOPO 4488.000 2000.000 128000.0 4551.0000 12 RR RL LR LL
1 EVLA_C#A0C0#1 64 TOPO 4616.000 2000.000 128000.0 4679.0000 12 RR RL LR LL
2 EVLA_C#A0C0#2 64 TOPO 4744.000 2000.000 128000.0 4807.0000 12 RR RL LR LL
3 EVLA_C#A0C0#3 64 TOPO 4872.000 2000.000 128000.0 4935.0000 12 RR RL LR LL
4 EVLA_C#A0C0#4 64 TOPO 5000.000 2000.000 128000.0 5063.0000 12 RR RL LR LL
5 EVLA_C#A0C0#5 64 TOPO 5128.000 2000.000 128000.0 5191.0000 12 RR RL LR LL
6 EVLA_C#A0C0#6 64 TOPO 5256.000 2000.000 128000.0 5319.0000 12 RR RL LR LL
7 EVLA_C#A0C0#7 64 TOPO 5384.000 2000.000 128000.0 5447.0000 12 RR RL LR LL
8 EVLA_C#B0D0#8 64 TOPO 5488.000 2000.000 128000.0 5551.0000 15 RR RL LR LL
9 EVLA_C#B0D0#9 64 TOPO 5616.000 2000.000 128000.0 5679.0000 15 RR RL LR LL
10 EVLA_C#B0D0#10 64 TOPO 5744.000 2000.000 128000.0 5807.0000 15 RR RL LR LL
11 EVLA_C#B0D0#11 64 TOPO 5872.000 2000.000 128000.0 5935.0000 15 RR RL LR LL
12 EVLA_C#B0D0#12 64 TOPO 6000.000 2000.000 128000.0 6063.0000 15 RR RL LR LL
13 EVLA_C#B0D0#13 64 TOPO 6128.000 2000.000 128000.0 6191.0000 15 RR RL LR LL
14 EVLA_C#B0D0#14 64 TOPO 6256.000 2000.000 128000.0 6319.0000 15 RR RL LR LL
15 EVLA_C#B0D0#15 64 TOPO 6384.000 2000.000 128000.0 6447.0000 15 RR RL LR LL
Initial Data Inspection
Since we have obtained the calibrated visibilites for the calibrator fields, we can now take this opportunity to investigate at high signal-to-noise the magnitude and timescale of the phase fluctuations we will be trying to correct for with selfcal.
Looking at the output of listobs we see that there is a long scan on the amplitude calibrator, 3C 286 (field ID 2). We will plot the calibrated phase vs. time on a single baseline:
# in CASA
plotms(vis='17B-197.sb34290063.eb34589992.58039.86119096065.ms', xaxis='time', yaxis='phase', ydatacolumn='corrected', field='2', antenna='ea05', correlation='RR,LL', avgchannel='64', iteraxis='baseline', coloraxis='spw')
- vis='17B-197.sb34290063.eb34589992.58039.86119096065.ms' : To plot visibilities from the pipeline calibrated MS.
- xaxis='time', yaxis='phase' : To set time as the x-axis and phase as the y-axis of the plot.
- ydatacolumn='corrected' : To plot the calibrated data (from the CORRECTED_DATA column).
- field='2' : To select visibilities from field ID 2, i.e., the amplitude calibrator 3C 286.
- antenna='ea05', iteraxis='baseline' : To view a single baseline at a time.
- correlation='RR,LL' : To plot both parallel-hand correlation products.
- avgchannel='64' : To average all channels in each SPW to increase signal-to-noise. Since the bandpass solutions have been applied to these data the channels will average coherently.
- coloraxis='spw' : To plot each SPW as a different color, which will make it easier to distinguish them.
Use the 'Next Iteration' button of the plotms GUI to cycle through additional baselines. You should see plots that look similar to the example image of the ea04-ea05 baseline. The plotted data have a mean of zero phase because the pipeline calibration solutions have already been applied. The phase is seen to vary with time over a large range (in some cases more than +/- 100 degrees) and the variations appear to be smooth over time scales of a few integrations. All of the spectral windows and both correlations follow the same trend with time. Additionally, the magnitude of the phase variations is larger for the higher frequency spectral windows, a pattern that is consistent with a change in atmospheric density.
Optional extra steps: Create and inspect similar plots using scan 2 of the phase calibrator field (J1549+5038). Repeat for baselines to other antennas.
Splitting the Target Visibilities
It is essential that we split the calibrated visibilities for the target we want to self-calibrate, meaning that the visibilities of the target source get copied from the CORRECTED_DATA column of the pipeline calibrated MS to the DATA column of a new measurement set. CASA calibration tasks always operate by comparing the visibilities in the DATA column to the source model, where the source model is given by either the MODEL_DATA column, a model image or component list, or the default model of a 1 Jy point source at the phase center. In the same way that the calibration pipeline used the raw visibilities in the DATA column to solve for calibration tables and then created the CORRECTED_DATA column by dividing the DATA column by these tables, self-calibration will work by comparing the pipeline calibrated visibilities (in the DATA column of the split MS) to a model, solving for self-calibration tables, and then creating a new CORRECTED_DATA column by applying the self-calibration tables.
# in CASA
split(vis='17B-197.sb34290063.eb34589992.58039.86119096065.ms',datacolumn='corrected',field='1', correlation='RR,LL', outputvis='obj.ms')
- vis='17B-197.sb34290063.eb34589992.58039.86119096065.ms' : The input visibilities for split. Here, these are the visibilities produced by the pipeline.
- datacolumn='corrected' : To copy the calibrated visibilities from the input MS.
- field='1' : The field ID of the target we want to self-calibrate.
- correlation='RR,LL' : To select only the parallel hand correlations. This will make the output data set smaller by about a factor of two.
- outputvis='obj.ms' : The name of the new measurement set that split will create.
-- The output MS can be can be directly downloaded here: obj.ms (3.4 GB) --
When to Use Self-Cal
EM: I think it would be really useful to explain in some detail when self-cal can help and when it can't as we have done below. But also why did we know that self-cal would work in this case.
There are instances in which self-cal can improve your image quality and others when it will not.
A couple typical cases in which self-cal can help improve the image of the target source:
- Extensive artifacts from the source of interest due to calibration errors
- Extensive artifacts from a background source due to direction-independent calibration errors
Some cases in which self-cal will *not* improve the image of the target source:
- When the image artifacts are due to errors in creating the image (e.g., ignoring wide-field effects)
- When the image artifacts are due to errors in deconvolution (e.g., ignoring wide-band effects)
- When the source(s) are too weak to achieve sufficient signal-to-noise in the calibration solutions
- When there is a bright outlying source with direction-dependent calibration errors
EM: could we expound upon how one could tell the difference between the cases in which self-cal WILL help and bullets 1,2, and 4 when self-cal will not help? What do those actually look like? People who are just starting may not know how to tell the difference between these different types of errors. To be honest, I am not sure I could tell you what the last 4 bullets looks like.
In the case of this guide, we knew that self-cal would improve the image quality because after an initial clean of the data, there were extensive artifacts from the source of interest due to calibration errors (). These errors manifested as.
The Initial Model
Preliminary Imaging
Prior to solving for self-calibration solutions we need to make an initial model of the target field. We will generate this model by deconvolving the target field using the task tclean. But before doing so, there are several imaging considerations that we should address prior to making this model (discussed below):
Image field-of-view: Ideally, we want our self-calibration model to include all of the sources present in the data (e.g., sources near the edge of the primary beam or in the first sidelobe). This is typically achieved by making an image large enough to encompass all of the apparent sources, or by making a smaller image of the target plus one or more outlier fields. We will start with a large dirty image of the entire primary beam (PB) in order to better understand the sources in the galaxy cluster plus any background sources that will need to be cleaned. A rule of thumb for the VLA is that the FWHM of the PB in arcminutes is approximately 45 * (1 GHz / nu). For our C-band observations, the VLA primary beam is ~8' FWHM. In order to image the entire PB and the first sidelobe we need an image field of view that is about four times larger, so we will choose a 32' field-of-view for our initial image. (see ... for a more complete discussion of the calculation of primary beam sizes)
Image cell size: There are a few different ways to estimate the synthesized beam size for these observations taken with the C-band in the B-configuration. For one, we can use this table in the VLA Observational Status Summary, which gives a resolution of 1.0" (see ... for a more complete discussion of the calculation of synthesized beam sizes). It is recommended to choose a cell size that will result in at least 5 image pixels across the FWHM of the synthesized beam, therefore we require a cell size of 0.20"/pixel or smaller.
Image size in pixels: We can convert our desired field-of-view to pixels using the cell size: 32' * (60" / 1') * (1 pixel / 0.20") = 9600 pixels. However, when making large images, CASA will run faster if we choose an image size that is optimized for the FFT algorithm. The recommendation is to choose an image size that can be expressed as 5*2^n*3^m. This is so that after CASA applies an internal padding factor of 1.2, the images being used in the FFTs can be broken down into small 2x2 and 3x3 matrices. So we will choose an optimized image size of 9720 (n=3,m=5) for this image.
Wide-field effects: Large images may require additional consideration due to non-coplanar baselines (the W-term). In CASA, this is usually addressed by turning on the W-project algorithm. See this link to CASAdocs for a more detailed discussion.
We can estimate whether our image requires W-projection by calculating the recommended number of w-planes using this formula taken from page 392 of the NRAO 'white book',
[math]\displaystyle{ N_{wprojplanes} = \left ( \frac{I_{FOV}}{\theta_{syn}} \right ) \times \left ( \frac{I_{FOV}}{1\, \mathrm{radian}} \right ) }[/math]
where I_FOV is the image field-of-view and theta_syn is the synthesized beam size. This evaluates to wprojplanes ~ 18 so we will choose to turn on the w-project algorithm with gridder='widefield' and set wprojplanes=18.
We will now create a preliminary dirty image using these parameters.
# in CASA
tclean(vis='obj.ms',imagename='obj.dirty.9720pix',imsize=9720, cell='0.2arcsec', pblimit=-0.1, gridder='widefield', wprojplanes=18 )
- cell='0.2arcsec': The size of an image pixel (see above).
- imsize=9720: The number of pixels across one side of the (square) image.
- pblimit=-0.1: We set this to a small negative number to turn off the PB mask.
- gridder='widefield': To turn on the w-project algorithm.
- wprojplanes=18: The number of w-planes to use for w-projection.
Figure 2A below shows the resulting dirty image, and Figure 2B shows a zoom-in of the central region of the image. Several outlying sources are detectable; the four brightest are marked with magenta circles.
 |
 |
|---|
We have a few options about how to deal with these outlying sources:
- Proceed with the self-calibration procedure using a large field-of-view that includes all the outlying sources.
- Proceed with the self-calibration procedure using a small field-of-view that includes only the central sources, and add an outlier field on each of the outlying sources.
- Model and UV-subtract the outlying sources.
- Peel the outlying sources.
- Proceed with the self-calibration procedure using a small field-of-view that includes only the central sources, ignoring the outlying sources.
In this guide, we will first choose to ignore the outlying sources in order to present a simplified self-calibration procedure. At the end of this guide, we will re-visit this choice and show examples of more advanced procedures.
Creating the Initial Model
After deciding how to deal with the outlying sources, our next step is to make the initial model that we will use for self-calibration. There are a few additional parameters that apply to this step, discussed below.
Image field-of-view: For this science case we are only interested in sources within ~1.5' of the cluster center. Since we have chosen to ignore the outlying sources at this stage, we will proceed with an image field-of-view of 3'.
Wide-field effects: We repeat the calculation of wprojplanes from the Initial Imaging section, which results in wprojplanes ~ 1. Therefore, we turn off the correction for non-coplanar baselines by setting gridder='standard' .
Wide-band imaging: Our images will combine data from all spectral windows, spanning a frequency range of about 4.5-6.5 GHz (a fractional bandwidth of about 36%). Each source's amplitude may vary substantially over this frequency range, due to either the source's intrinsic spectral variation and/or the frequency dependence of the VLA's primary beam. To mitigate these errors during deconvolution we will use deconvolver='mtmfs' and nterms=2. For further discussion of wide-band imaging, see the CASA documentation for wide-band imaging and the VLA Imaging CASAguide.
Imaging weights: When constructing the initial model, especially when there are large image artifacts, it is recommended to use "robust" imaging weights. In CASA, this is enabled with weighting='briggs' , and then choosing a value for the robust parameter. Large positive values of robust often result in large positive PSF sidelobes while large and negative values of robust often produce large negative PSF sidelobes. Since many types of image artifacts scale with PSF sidelobe levels, a reasonable compromise is often around robust=0.
Image deconvolution: We will need to deconvolve (clean) this image in order to produce a model of the field. We will want to control the cleaning depth and masking interactively, so we set interactive=True. We also must choose the number of clean iterations with the niter parameter-- 1000 iterations is a suggested starting value, and this can be changed interactively after we start cleaning.
Saving the model: After deconvolution, there are a couple options for how to save the model, controlled by the savemodel parameter. The default is none which will not save the model. It is essential that this default is changed or else the selfcal procedure will fail. The option savemodel='virtual' will save a copy of the image to the SOURCE subtable of the measurement set to be used later for on-the-fly model visibility prediction. This option is sometimes recommended for very large data sets. The other option, savemodel='modelcolumn' is the recommended setting and the one that we will use in this guide. This option will do the model visibility prediction after cleaning and save the result to the MODEL_DATA column.
Now we are ready to create our first clean image. This image will provide the starting model that is required by the calibration routines, and it will showcase why we need self-calibration for these data.
# in CASA
tclean(vis='obj.ms',imagename='obj.prelim_clean.3am', gridder='standard', cell='0.2arcsec', imsize=900, pblimit=-0.1, deconvolver='mtmfs', nterms=2, niter=1000, interactive=True, weighting='briggs', robust=0, savemodel='modelcolumn')
- gridder='standard': We select the default gridder (assumes coplanar baselines).
- cell='0.2arcsec': The size of an image pixel.
- imsize=900: To create an image with a 3' field-of-view.
- pblimit=-0.1: To remove the PB mask; see previous sections.
- deconvolver='mtmfs': We will turn on the wide-band deconvolution algorithm as discussed above.
- nterms=2: The number of Taylor terms for wide-band deconvolution.
- niter=1000: Set a relatively large number of iterations as a starting point.
- interactive=True: So we can interactively place the mask.
- weighting='briggs': Turn on 'robust' image weighting
- robust=0: Set the value of the robust weighting.
- savemodel='modelcolumn': To enable writing the MODEL_DATA column after imaging. **important**
During the interactive cleaning, we place circular masks around each of the strong sources in turn:
- The rightmost source (Figure 2A) -- let it run the first set of clean iterations by pressing the green circle arrow in the CASA viewer.
- The leftmost double-lobed source (Figure 2B) -- let it run an additional set of clean iterations by pressing the green circle arrow in the CASA viewer.
- The middle source (Figure 2C) -- continue cleaning inside the masks by pressing the green circle arrow in the CASA viewer.
At this point, the residual emission is at about the same level as the artifacts so we stop cleaning (a total of about 225 iterations; see Figure 2D for an example of the artifacts). It is strongly recommended to only mask and clean emission that is believed to be real so as not to include artifacts in the model. Figure 2E shows the resulting clean image that we will try to improve through the use of self-calibration.
 |
 |
 |
 |
 |
|---|
Verifying the Initial Model
There have been reported instances where CASA fails to save the model visibilities when using interactive clean. It is crucial that the model is saved correctly, otherwise self-calibration will use the 'default' model of a 1 Jy point source at the phase center. The default model may be very different from your target field and we do not want to carry out the self-cal procedure with this incorrect model. Therefore, it is recommended to verify that the model has been saved correctly by inspecting the model visibilities.
# in CASA
plotms(vis='obj.ms', xaxis='UVwave', yaxis='amp', ydatacolumn='model', correlation='RR,LL', avgchannel='64', avgtime='300s')
- vis='obj.ms' : To plot visibilities from the split MS.
- xaxis='UVwave', yaxis='amp' : To set UV-distance in wavelengths as the x-axis and amplitude as the y-axis of the plot.
- ydatacolumn='model' : To plot the model visibilities (from the MODEL_DATA column).
- correlation='RR,LL' : To only plot the parallel-hand visibilities. The model RL=LR=0 since we only made a Stokes I model image.
- avgchannel='64' : To average all channels per SPW.
- avgtime='300s' : To average in time in chunks of 300 seconds.
The resulting plot should resemble Figure 3 on the right. This plot shows that some baselines see up to 15 mJy of flux, but that the source becomes resolved on the longer baselines.
This model that has been plotted is clearly not the default model of a 1 Jy point source and so we have verified that tclean has correctly written the MODEL_DATA column.
First Round of Self-Calibration
First Self-Calibration Table
For this first round of self-cal we will use the model that we just created above, so we are now ready to solve for the first self-cal solutions. We will explore various parameters of the task gaincal in order to learn more about the data and settle on the optimal parameters. The most relavant parameters are discussed below:
Solution interval: This is controlled with the solint parameter and is one of the most fundamental parameters for self-calibration. The value of this parameter can vary between 'int' for a single integration (corresponding to 3 seconds for this data set) up to 'inf' for infinite (meaning an entire scan or the entire observation, depending on the value of the combine parameter) . We typically want to choose the shortest solution interval for which we can achieve adequate signal-to-noise in the calibration solutions.
Data combination: The data can be combined in multiple ways to improve signal-to-noise, but if the target source is bright enough to obtain good calibration solutions in a short timescale without data combination then these options need not be used. One can try to combine multiple SPWs with combine='spw' if the SPWs are at similar frequencies, and can generally expect to increase the solution's signal-to-noise by the square root of the number of SPWs that are combined. Both parallel-hand corrleations, if present, can be combined by setting gaintype='T' instead of gaintype='G' and this will generally increase the signal-to-noise by a factor of root 2. Combining scans is not usually recommended.
Amplitude and phase correction: Because large phase errors will result in incoherent averaging and lead to lower amplitudes, we always want to start with phase-only self-calibration. We achieve this by setting calmode='p' . In later rounds of self-cal, after the phases have been well corrected, we can try calmode='ap' to include an amplitude component to the solutions. We may also want to consider normalizing the amplitudes using the solnorm parameter.
Reference antenna: As with standard calibration, we want to choose a reference antenna for the calibration solutions. One that is near the center of the array, but not flagged due to shadowing, is generally recommended.
Signal-to-noise ratio (SNR): The default minimum SNR in gaincal is 3.0, but this can be adjusted with the minsnr parameter. Solutions below this minimum are flagged in the output calibration table. Sometimes we want to increase this minimum, e.g., to 5.0, to reject noisy solutions. Alternatively, we may want to lower this minimum, e.g., to zero, usually for inspection purposes.
We will now create our initial self-calibration table. This will not be the final table for the first round of self-calibration, but rather, a temporary table that we will inspect to help determine the optimal parameters.
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_initial.tb',solint='int',refant='ea24',calmode='p',gaintype='G',minsnr=0)
- caltable='selfcal_initial.tb':'Name the calibration tables something intuitive to distinguish each one.
- solint='int': We choose a solution interval equal to the integration time (3 seconds) in order to get a sense of the structure and timescale of the variations.
- refant='ea24': The chosen reference antenna.
- calmode='p': To start with phase only calibration.
- gaintype='G': To solve for the polarizations separately.
- minsnr=0: To turn off flagging of low-SNR solutions, so that we can see all the solutions.
You may see several messages printed to the terminal while gaincal is running. There are a couple different types of messages, e.g.,
Found no unflagged data at: (time=2017/10/13/20:50:07.5 field=0 spw=0 chan=0)
This means that all the input data was flagged for this solution interval. This is generally harmless unless there are far fewer solutions in the output table than you were expecting.
Another message you may see printed to the terminal looks like:
13 of 54 solutions flagged due to SNR < 3 in spw=0 at 2017/10/13/20:58:13.5
This is telling us that solutions have been flagged for being below the minimum signal-to-noise ratio set by the minsnr parameter. You will not see these messages in the first execution of gaincal because we have set minsnr=0, but watch out for these in subsequent executions. In the above example, 13 solutions are flagged out of a total of 54 (one per antenna per polarization) in SPW 0 and for a single 3 second solution interval at the reported time. If we had instructed gaincal to combine polarizations with gaintype='T' then there would have been 27 total solutions. You may see several of these messages, in which case try to determine if they correspond to the same antenna, SPW or time as this may indicate the presense of bad data. If these messages appear across antennas, times and SPWs then this likely indicates that the signal-to-noise is too low and that more data needs to be combined.
It is recommended to also check the logger messages written by gaincal to find the total number of solutions, i.e.,
INFO gaincal Calibration solve statistics per spw: (expected/attempted/succeeded): INFO gaincal Spw 0: 561/522/522 INFO gaincal Spw 1: 561/522/522 ... ... ...
This shows that gaincal found solutions for most of the solution intervals.
Plotting the First Self-Calibration Table
To view these solutions, we use plotms.
# in CASA
plotms(vis='selfcal_initial.tb',xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=2, gridcols=2, coloraxis='corr')
- xaxis='time' & yaxis='phase' : View the phase variations over time with respect to antenna 24.
- iteraxis='antenna' : Create separate plots of the corrections for each antenna.
- gridrows=3 & gridcols=3: It can be helpful to view multiple plots at once, as we will be stepping through several plots. In this case, 9 plots per page.
- coloraxis='corr' : To use different colors when plotting different polarizations (R and L will be black and red, respectively).
Iterate through these plots using the 'Next Iteration' button (green triangle) to inspect the solutions for all antennas. Some noteworthy observations include:
- There are large, coherent phase changes of more than 100 degrees.
- The timescale of these changes is fairly short, about 20 seconds.
- The scatter in these signals is high (approximately a few 10s of degrees), indicating low signal-to-noise.
- The phase changes in the two polarizations appear to match each other.
- There is some interesting behavior in the 3rd scan for antennas ea15 and ea27.
Note: The plot of the reference antenna, ea24, is not unusual and should be considered to be consistent with zero phase.
It is apparent from these plots that we can combine polarizations to improve the solution signal-to-noise ratio, since we observed that the solutions for the two polarizations were very similar. The next thing we want to understand is if we can combine SPWs, and if so, which ones. We can plot the previous solutions in a slightly different way to help answer this question.
We will view these solutions again using plotms, but this time we will color the solutions by SPW.
# in CASA
plotms(vis='selfcal_initial.tb',xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3, coloraxis='spw')
- xaxis='time' & yaxis='phase' : View the phase variations over time with respect to antenna 24.
- iteraxis='antenna' : Create separate plots of the corrections for each antenna.
- gridrows=3 & gridcols=3: It can be helpful to view multiple plots at once, as we will be stepping through several plots. In this case, 3 plots per page.
- coloraxis='spw' : To use different colors when plotting different SPWs.
Iterate again through these plots using the 'Next Iteration' button (green triangle) to inspect the solutions for all antennas. When you get to ea15 it should be clear that the solutions are not the same for all SPWs. This is also true but less obvious for ea27 due to the limited number of colors available to plotms.
We can inspect this further in the following plot:
# in CASA
plotms(vis='selfcal_initial.tb',xaxis='time',yaxis='phase',antenna='ea15',iteraxis='spw',gridrows=3, gridcols=3, coloraxis='scan')
- xaxis='time' & yaxis='phase' : View the phase variations over time with respect to antenna 24.
- antenna='ea15' : To select only antenna ea15.
- iteraxis='spw' : Create separate plots of the corrections for each SPW.
- gridrows=3 & gridcols=3: To view multiple plots at once. In this case, 9 plots per page.
- coloraxis='scan' : To use different colors when plotting different scans.
The first set of 9 plots should have a similar pattern. On the next iteration, this pattern should continue for SPWs 9~11, but then change for SPWs 12~15. The signal-to-noise for SPW 13 is also noticeably lower. If we create these plots for ea27 we will see a similar pattern, only this time the pattern is constant over SPWs 0~5 and then it changes to a new pattern that is constant for SPWs 6~15. For both ea15 and ea27, the change only happens in the third of the three scans.
Unfortunately, since all SPWs do not show the same phase solutions, it will not be trivial to combine SPWs to increase signal-to-noise. We will continue without combining SPWs, but in Appendix ... we show a more complicated example of how to handle this situation.
Self-calibration Tables of Various Solution Intervals
Now that we have made the decision about how to handle SPW combination, we will move on to consider time averaging. We observed the previous solutions to display large, coherent phase changes but also to have significant scatter due to low signal-to-noise. We could increase the signal-to-noise by root 2 for each doubling of the solution interval, but it does not make sense to average over timescales larger than the characteristic time over which the phase remains constant. In this section, we will demonstrate this effect by creating and plotting tables over a range of solution intervals. We will also combine both polarizations to improve the solution signal-to-noise ratio, since we observed the two polarizations to measure the same phase changes (with gaintype='T' ).
This table will have a solution interval of 3 seconds (the integration time).
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_pol_solint_3.tb',solint='int',refant='ea24',calmode='p',gaintype='T', minsnr=0)
plotms(vis='selfcal_combine_pol_solint_3.tb',xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3, coloraxis='spw')
This table will have a solution interval of 6 seconds (2 times the integration time).
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_pol_solint_6.tb',solint='6s',refant='ea24',calmode='p',gaintype='T', minsnr=0)
plotms(vis='selfcal_combine_pol_solint_6.tb',xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3, coloraxis='spw')
This table will have a solution interval of 12 seconds (4 times the integration time).
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_pol_solint_12.tb',solint='12s',refant='ea24',calmode='p',gaintype='T', minsnr=0)
plotms(vis='selfcal_combine_pol_solint_12.tb',xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3, coloraxis='spw')
This table will have a solution interval of 24 seconds (8 times the integration time).
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_pol_solint_24.tb',solint='24s',refant='ea24',calmode='p',gaintype='T', minsnr=0)
plotms(vis='selfcal_combine_pol_solint_24.tb',xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3, coloraxis='spw')
This table will have a solution interval of 48 seconds (16 times the integration time).
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_pol_solint_48.tb',solint='48s',refant='ea24',calmode='p',gaintype='T', minsnr=0)
plotms(vis='selfcal_combine_pol_solint_48.tb',xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3, coloraxis='spw')
This table will have a solution interval of 96 seconds (32 times the integration time).
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_pol_solint_96.tb',solint='96s',refant='ea24',calmode='p',gaintype='T', minsnr=0)
plotms(vis='selfcal_combine_pol_solint_96.tb',xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3, coloraxis='spw')
We will now investigate which solution interval (solint) is optimal for this particular dataset by comparing the corrections calculated with several different solution intervals. We can see that several of the shorter timescale solutions encapsulate the structure of the phase variations, and we can see the signal-to-noise of the solutions improving with longer timescales. But in the very long timescale solutions, we can see that the phase is varying faster than the solution interval and so the solutions no longer do a good job of capturing the phase changes that we are trying to correct for.
If we were to apply the low signal-to-noise solutions, then we would, on average, be correcting for the phase changes but we would also introduce random phase errors that would be impossible to remove, and that would reduce the sensitivity of our observation. Alternatively, if we were to apply the longest timescale solutions, then we would be leaving residual phase errors in the corrected data. Usually, the best choice is to identify the longest (highest signal-to-noise) solution interval that is capable of describing the structure of the phase variations. For the examples above, this appears to be between about 12 and 24 seconds.
This is particularly obvious if we overplot the solutions using plotms (is this possible from the command line?). The figure to the right compares the shortest solution interval (3 seconds) with the 12 second solution interval.
Applying the First Self-Calibration Table
# in CASA
applycal(vis='obj.ms',gaintable='sc1_tb_15s.gain')
Other potential intervals that could be useful to use are 5 seconds and 60 seconds. Additionally, if the plotcal reveals no structure in these intervals, you can combine the spws together and apply a combined solution by setting combine='spw' (you will have to set spwmap=[0 x number of spectral windows] in apply cal in order to tell CASA to apply the combined solution to all of the spectral windows -- in this case spwmap = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0])
Our goal is to create a better model of the source using the data itself to create a better image. The process with calculating corrections and applying them has (if done correctly) reduced the errors allowing us to make a better image and create a better model. applycal applies the solutions to the (raw) MS DATA column and writes the calibrated data into the CORRECTED_DATA column, thus for determining the next order/level of corrections during the next round of calibration we need to make the corrected data column the data column.
# in CASA
split(vis='obj.ms',outputvis='sc1_obj.ms')
Second Round of Self-Calibration
We will go through the same process again: image, gaincal, applycal, and image. This time the image we create at first will be an indication if the solutions we applied made a difference.
# in CASA
tclean(vis='obj.ms',imagename='obj.sc1',imsize=750,cell='0.24arcsec',weighting='briggs',deconvolver='mtmfs',niter=500,savemodel='modelcolumn',interactive=True)
Larger mask.
# in CASA
gaincal(vis='obj.ms',caltable='sc2_tb_15s',solint='15s',refant='ea24',calmode='p',gaintype='T')
- solint='15s': Since we found a good interval of 15 seconds previously, we move forward with 15s.
Josh: do we want to suggest/show other solution intervals? Do we want to comment on using various intervals -- when is it useful to do several rounds with one interval versus several rounds with each a different solint?
View the solutions.
# in CASA
plotcal(caltable='sc1_tb_int.gain',xaxis='time',yaxis='phase',iteration='antenna',subplot=321)
We can see our corrections made a difference. In particular notice the amplitude of the phase variations, instead of +_60-80 degress they are fractions of a degree. This indicates that we have corrected for a majority of the variations.
When to stop: (looking at amplitude of corrections, image has improved/dynamical range)
Notes about amplitude corrections (?) (mode='ap' or 'a').
When self-cal helped, but didn't improve the image enough: (suggestions for other things to do)
Appendix
PB calculations
PSF calculations
Variant of the Self-cal Procedure
Note: There is an alternative and equally valid variant of the self-calibration procedure that is outlined below but not used in this guide. This variant avoids the use of split to freeze-in the self-cal solutions and instead builds up a list containing the calibration table(s) from each round of self-cal.
- Conservatively clean the target and save the model.
- Use gaincal with a set of parameters to calculate a calibration table, pre-applying all previous self-cal tables 'on-the-fly' with the gaintable parameter.
- Inspect the new calibration solutions using plotms.
- Determine if the solutions should be applied, re-calculated with different parameters, or if it is time to stop.
- Append the new solutions to the list of calibration tables.
- Use applycal to apply the list of tables to the data.
- Start the next round of self-cal.
Using combine='spw' and spwmap
Appendix about combine='spw'
From listobs, the three target scans are: 3,5,7
we can use combine='spw' for scans 3 and 5
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_spw_scan_3and5.tb',solint='int',refant='ea24',calmode='p',gaintype='T',combine='spw',scan='3,5')
- caltable='selfcal_combine_spw_scan_3and5.tb':'Name the calibration tables something intuitive to distinguish each one.
- solint='int': We choose a solution interval equal to the integration time (3 seconds) in order to get a sense of the structure and timescale of the variations.
- refant='ea24': The chosen reference antenna.
- calmode='p': To start with phase only calibration.
- gaintype='T': To combine polarizations.
- combine='spw' : To also combine SPWs.
- scan='3,5' : To select only scans 3 and 5.
for scan 7, the solutions are similar over SPWs 0~5 and over SPWs 12~15 so these can be combined separately.
We can do this by running gaincal multiple times and setting append=True.
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_spw_scan7.tb',solint='int',refant='ea24',calmode='p',gaintype='T',combine='spw',spw='0~5',scan='7',append=True)
- caltable='selfcal_combine_spw_scan7.tb':'Name the calibration tables something intuitive to distinguish each one.
- solint='int': We choose a solution interval equal to the integration time (3 seconds) in order to get a sense of the structure and timescale of the variations.
- refant='ea24': The chosen reference antenna.
- calmode='p': To start with phase only calibration.
- gaintype='T': To combine polarizations.
- combine='spw' : To also combine SPWs.
- spw='0~5' : To select only SPWs 0~5.
- scan='7' : To select only scan 7.
- append=True: To append solutions to the existing table of the same name.
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_spw_scan7.tb',solint='int',refant='ea24',calmode='p',gaintype='T',combine='spw',spw='6~11',scan='7',append=True)
- caltable='selfcal_combine_spw_scan7.tb':'Name the calibration tables something intuitive to distinguish each one.
- solint='int': We choose a solution interval equal to the integration time (3 seconds) in order to get a sense of the structure and timescale of the variations.
- refant='ea24': The chosen reference antenna.
- calmode='p': To start with phase only calibration.
- gaintype='T': To combine polarizations.
- combine= : To turn off combination of SPWs.
- spw='6~11' : To select only SPWs 6~11.
- scan='7' : To select only scan 7.
- append=True: To append solutions to the existing table of the same name.
# in CASA
gaincal(vis='obj.ms',caltable='selfcal_combine_spw_scan7.tb',solint='int',refant='ea24',calmode='p',gaintype='T',combine='spw',spw='12~15',scan='7',append=True)
- caltable='selfcal_combine_spw_scan7.tb':'Name the calibration tables something intuitive to distinguish each one.
- solint='int': We choose a solution interval equal to the integration time (3 seconds) in order to get a sense of the structure and timescale of the variations.
- refant='ea24': The chosen reference antenna.
- calmode='p': To start with phase only calibration.
- gaintype='T': To combine polarizations.
- combine='spw' : To turn back on combination of SPWs.
- spw='12~15' : To select only SPWs 12~15.
- scan='7' : To select only scan 7.
- append=True: To append solutions to the existing table of the same name.
Note that when combining SPWs, only one set of solutions are written. These solutions are indexed by the lowest numbered SPW that was used. So the solutions created by combining SPWs 0~5 were written out as the solutions of SPW 0, and the solutions created by combining SPWs 12~15 were written out as the solutions of SPW 12. If we were to apply these solutions in the typical fashion, applycal would not see solutions for SPWs 1~5 or 13~15 and the data would not be corrected. Therefore, we must use the spwmap parameter of applycal to map the combined solutions back onto all the SPWs from which they were derived. The following calls to applycal demonstrate this.
First, apply the solutions from scans 3 and 5, which combined all of the SPWs.
# in CASA
applycal(vis='obj.ms',gaintable='selfcal_combine_spw_scan_3and5.tb',interp='nearest',scan='3,5',spw='',spwmap=16*[0])
- gaintable='selfcal_combine_spw_scan_3and5.tb':'Name the calibration tables something intuitive to distinguish each one.
- interp='nearest': Correct the visibilities using the solution from the closest time interval.
- scan='3,5' : To select scans 3 and 5 for correction.
- spw= : To select all SPWs for correction.
- spwmap=16*[0]: To map the solutions from SPW 0 to all 16 SPWs.
Then apply the solutions from scan 7. We will need to map the solutions from SPW 0 to SPWs 0~5, the solutions from SPWs 6~11 to themselves, and the solutions from SPW 12 to SPWs 12~15.
# in CASA
applycal(vis='obj.ms',gaintable='selfcal_combine_spw_scan7.tb',interp='nearest',scan='7',spw='',spwmap=6*[0]+range(6,12)+4*[12])
- gaintable='selfcal_combine_spw_scan7.tb':'Name the calibration tables something intuitive to distinguish each one.
- interp='nearest': Correct the visibilities using the solution from the closest time interval.
- spw= : To select all SPWs for correction.
- scan='7' : To select scans 3 and 5 for correction.
- spwmap=6*[0]+range(6,12)+4*[12]: To map the solutions from SPW 0 to all 16 SPWs.
