The workhorse for synthesis calibration is the cb tool.
This chapter explains how to calibrate interferometer data within the CASA task system. Calibration is the process of determining the net complex correction factors that must be applied to each visibility in order to make them as close as possible to what an idealized interferometer would measure, such that when the data is imaged an accurate picture of the sky is obtained. This is not an arbitrary process, and there is a philosophy behind the CASA calibration methodology (see § 4.2.1 for more on this). For the most part, calibration in CASA using the tasks is not too different than calibration in other packages such as AIPS or Miriad, so the user should not be alarmed by cosmetic differences such as task and parameter names!
Alert: The calibration table format changed in CASA 3.4. CASA 4.2 is the last version that will support the caltabconvert function that provides conversions from the pre-3.4 caltable format to the modern format; it will be removed for CASA 4.3. In general, it is best to recalculate calibration using CASA 3.4 or later.
Alert: In CASA 4.2 the gaincurve and opacity parameters have been removed from all calibration tasks (as advertised in 4.1). These calibration types are supported via the gencal task.
Alert: As part of continuing development of a more flexible and improved interface for specifying calibration for apply, a new parameter has been introduced in applycal and the solving tasks: docallib. This parameter toggles between use of the traditional calibration apply parameters (gaintable, gainfield, interp, spwmap, and calwt), and a new callib parameter which currently provides access to the experimental Cal Library mechanism, wherein calibration instructions are stored in a file. The default remains docallib=False in CASA 4.5, and this reveals the traditional apply parameters which continue to work as always, and the remainder of this chapter is still written using docallib=F. Users interested in the Cal Library mechanism’s flexibility are encouraged to try it and report any problems; see Appendix G for information on how to use it, including how to convert traditional applycal to Cal Library format. Note also that plotms and mstransform now support use of the Cal Library to enable on-the-fly calibration when plotting and generating new MSs.
The standard set of calibration solving tasks (to produce calibration tables) are:
There are helper tasks to create, manipulate, and explore calibration tables:
There are some development versions of calibration and utility tasks that are recently added to the suite:
These are not yet full-featured, and may have only rudimentary controls and options.
The following sections outline the use of these tasks in standard calibration processes.
Information on other useful tasks and parameter setting can be found in:
A work-flow diagram for CASA calibration of interferometry data is shown in Figure 4.1. This should help you chart your course through the complex set of calibration steps. In the following sections, we will detail the steps themselves and explain how to run the necessary tasks and tools.
Figure 4.1: Flow chart of synthesis calibration operations. Not shown are use of table manipulation and plotting tasks accum, plotcal, and smoothcal (see Figure 4.2).
This can be broken down into a number of discrete phases:
The flow chart and the above list are in a suggested order. However, the actual order in which you will carry out these operations is somewhat fluid, and will be determined by the specific data-reduction use cases you are following. For example, you may need to do an initial Gain Calibration on your bandpass calibrator before moving to the Bandpass Calibration stage. Or perhaps the polarization leakage calibration will be known from prior service observations, and can be applied as a constituent of Prior Calibration.
Calibration is not an arbitrary process, and there is a methodology that has been developed to carry out synthesis calibration and an algebra to describe the various corruptions that data might be subject to: the Hamaker-Bregman-Sault Measurement Equation (ME), described in Appendix E. The user need not worry about the details of this mathematics as the CASA software does that for you. Anyway, it’s just matrix algebra, and your familiar scalar methods of calibration (such as in AIPS) are encompassed in this more general approach.
There are a number of “physical” components to calibration in CASA:
At its most basic level, Calibration in CASA is the process of taking “uncalibrated” data, setting up the operation of calibration tasks using parameters, solving for new calibration tables, and then applying the calibration tables to form “calibrated” data. Iteration can occur as necessary, with the insertion of other non-calibration steps (e.g. imaging to generate improved source models for “self-calibration”).
Figure 4.2: Chart of the table flow during calibration. The parameter names for input or output of the tasks are shown on the connectors. Note that from the output solver through the accumulator only a single calibration type (e.g. ’B’, ’G’) can be smoothed, interpolated or accumulated at a time. accum is optional (and not recommended as of v4.0). The final set of cumulative calibration tables of all types (accumulated or as a list of caltables) are then input to applycal as shown in Figure 4.1.
The calibration tables are the currency that is exchanged between the calibration tasks. The “solver” tasks (gaincal, bandpass, blcal, polcal) take in the MS (which may have a calibration model attached) and previous calibration tables, and will output an “incremental” calibration table (it is incremental to the previous calibration, if any). This table can then be smoothed using smoothcal if desired.
You can optionally accumulate the incremental calibration onto previous calibration tables with accum, which will then output a cumulative calibration table. This task will also interpolate onto a different time scale. See § 4.5.6 for more on accumulation and interpolation.
Figure 4.2 graphs the flow of these tables through the sequence
solve => smooth => accumulate
Note that this sequence applied to separate types of tables (e.g. ’B’, ’G’) although tables of other types can be previous calibration input to the solver.
The final set of cumulative calibration tables is what is applied to the data using applycal. You will have to keep track of which tables are the intermediate incremental tables, and which are cumulative, and which were previous to certain steps so that they can also be previous to later steps until accumulation. This can be a confusing business, and it will help if you adopt a consistent table naming scheme (see Figure 4.2) for an example naming scheme).
CASA supports the calibration of traditional VLA data that is imported from the Archive through the importvla task. See § 2.2.3 for more information.
ALERT: Data taken both before and after the Modcomp turn-off in late June 2007 will be handled automatically by importvla. You do not need to set special parameters to do so, and it will obey the scaling specified by applytsys.
You can also import VLA data in UVFITS format with the importuvfits task (§ 2.2.7.1). However, in this case, you must be careful during calibration in that some prior or previous calibrations (see below) may or may not have been done in AIPS and applied (or not) before export.
For example, the default settings of AIPS FILLM will apply VLA gaincurve and approximate (weather-based) atmospheric optical depth corrections when it generates the extension table CL 1. If the data is exported immediately using FITTP, then this table is included in the UVFITS file. However, CASA is not able to read or use the AIPS SN or CL tables, so that prior calibration information is lost and must be applied during calibration here (i.e. using gaincurve=True and setting the opacity parameter).
On the other hand, if you apply calibration in AIPS by using the SPLIT or SPLAT tasks to apply the CL tables before exporting with FITTP, then this calibration will be in the data itself. In this case, you do not want to re-apply these calibrations when processing in CASA.
Jansky VLA data can be loaded into CASA either via importevla or by using the task importasdm. Both tasks will convert ASDM raw data files into Measurement Sets. importasdm will convert the data itself and the majority of the metadata. importevla will run importasdm followed by Jansky VLA-specific corrections, like the application of the on-line flags (e.g. times when the subreflector was not in place or the an antenna was not on source), an option to clip values that are exactly zero (as of 2010, such values still may appear in the VLA raw data), and flagging for shadowing.
There are a number of “a priori” calibration quantities that may need to be initialized or estimated before further calibration solving is carried out. These include
These are pre-determined effects and should be applied (if known) as priors when solving for other calibration terms, and included in the final application of all calibration. If unknown, then they will be solved for or subsumed in other calibration such as bandpass or gains.
We now deal with these in turn.
See Appendix F for a more complete description of weight accounting in CASA.
CASA 4.3 introduced initial experimental support for spectral weights. At this time, this is mainly relevant to ALMA processing for which spectral Tsys corrections, which faithfully reflect spectral sensitivity, are available. In most other cases, sensitivity is, to a very good approximation, channel-independent after bandpass calibration (and often also before), except perhaps at the very edges of spectral windows (and for which analytic expressions of the sensitivity loss are generally unavailable). Averaging of data with channel-dependent flagging which varies on sufficiently short timescales will also generate channel-dependent net weights (see split or mstransform for more details).
By default, CASA’s weight accounting scheme maintains unchannelized weight information that is appropriately updated when calibration is applied. In the case of spectral calibrations (Tsys and bandpass), an appropriate spectral average is used for the weight update. This spectral average is formally correct for weight update by bandpass. For Tsys, traditional treatments used a single measurement per spectral window; ALMA has implemented spectral Tsys to better track sensitivity as a function of channel, and so should benefit from spectral weight accounting as described here, especially where atmospheric emmission lines occur. If spectral weight accounting is desired, users must re-initialize the spectral weights using the initweights task:
initweight(vis='mydata.ms', wtmode='nyq', dowtsp=True)
In this task, the wtmode parameter controls the weight initialization convention. Usually, when initializing the weight information for a raw dataset, one should choose wtmode=’nyq’ so that the channel bandwidth and integration time information are used to initialize the weight information. The dowtsp parameter controls whether (T) or not (F) the spectral weights (WEIGHT_SPECTRUM column) are initialized. The default is dowtsp=False, wherein only the non-spectral weights (WEIGHT column) will be initialized. If the spectral weights have been initialized, then downstream processing that supports spectral weights will use/update them. In CASA 4.3 and later, this includes applycal, clean, and split/mstransform; use of spectral weights in calibration solving (e.g., gaincal and other solve tasks) is scheduled for the CASA 4.5 release.
Note that importasdm and importevla currently initialize the non-spectral weights using channel bandwidth and integration time information (equivalent to initweights(vis=’mydata.ms’,wtmode=’nyq’,dowtsp=F)). In general, it only makes sense to run initweights on a raw dataset which has not yet been calibrated, and it should only be necessary if the filled weights are inappropriate, or if spectral weight accounting is desired in subsequent processing. It is usually not necessary to re-initialize the weight information when redoing calibration from scratch (the raw weight information is preserved in the SIGMA/SIGMA_SPECTRUM columns). (Re-)initializing the weight information for data that has already been calibrated (with calwt=T, presumably) is formally incorrect and is not recommended.
When combining datasets from different epochs, it is generally preferable to have used the same version of CASA (most recent is best), and with the same weight information conventions and calwt settings. Doing so will minimize the likelihood of arbitrary weight imbalances that might lead to net loss of sensitivity, and maximize the likelihood that real differences in per-epoch sensitivity (e.g., due to different weather conditions and instrumental setups) will be properly accounted for. Modern instruments support more variety in bandwidth and integration time settings, and so use of these parameters in weight initialization is preferred (c.f. use of unit weight initialization, which has often been the traditional practice).
Alert: Full and proper weight accounting for the EVLA formally depends on the veracity of the switched power calibration scheme (§ 4.3.7). As of mid-2015, use of the EVLA switched power is not yet recommended for general use, and otherwise uniform weights are carried through the calibration process. As such, spectral weight accounting is not yet meaningful. Facilities for post-calibration estimation of spectral weights are planned for a future release.
Some telescopes, including the old VLA, ALMA, and the VLBA, record the visibilities in the form of raw correlation coefficient with weights proportional to the number of bits correlated. The correlation coefficient is the fraction of the total signal that is correlated, and thus multiplication by the system temperature (Tsys) and the antenna gain (in Jy/K) will produce visibilities with units of correlated flux density. ALMA records Tsys(K) information in the MS which can be extracted as a caltable using gencal with calmode=’tsys’, and applied to data to yield units of K. Calibration to flux density in Jy is achieved via reference to sources of known power.
Alert: Note that the old VLA system did this initial calibration on-line. The modern VLA does not record normalized visibilities. Instead, the correlations are delivered in raw engineering units that are proportional to power. The actual total power received is continuously monitored during the observation, with a calibration signal of known temperature (K) switched in at a rate of 10 Hz. This is the so-called “switched-power” calibration system on the VLA. This enables a continuous record of the Tsys(K), as well as net electronic gain variation of each antenna’s receiving system. The correlator requantizer gain is also monitored. These data are recorded in MS subtables and appropriate calibration factors can be derived from them by gencal with caltype=’swpow’, and stored in a caltable for application. This calibration is not a “Tsys” calibration of the traditional sort; the switched-power gain is used to correct the visibility amplitude, and the Tsys is used to set the weights. This system is still being commissioned (as of early 2014). Observations using 8-bit sampling are usually reasonably calibrated; 3-bit-sampled switched-power data are subject to compression effects that are not yet completely understood, and the switched power calibration is not recommended (instead, correction only by the requantizer gain is recommended, using caltype=’rq’).
See § 4.3.6 for more information on use of gencal.
Large antennas (such as the 25-meter antennas used in the VLA and VLBA) have a forward gain and efficiency that changes with elevation. Gain curve calibration involves compensating for the effects of elevation on the amplitude of the received signals at each antenna. Antennas are not absolutely rigid, and so their effective collecting area and net surface accuracy vary with elevation as gravity deforms the surface. This calibration is especially important at higher frequencies where the deformations represent a greater fraction of the observing wavelength. By design, this effect is usually minimized (i.e., gain maximized) for elevations between 45 and 60 degrees, with the gain decreasing at higher and lower elevations. Gain curves are most often described as 2nd- or 3rd-order polynomials in zenith angle.
Gain curve calibration has been implemented in CASA for the modern VLA and old VLA (only), with gain curve polynomial coefficients available directly from the CASA data repository. To make gain curve and antenna efficiency corrections for VLA data, use gencal with caltable=’gceff’. See § 4.3.6 for more information on use of gencal.
ALERT: If you are not using VLA data, do not use gaincurve corrections. A general mechanism for incorporating gaincurve information for other arrays will be made available in future releases. The gain-curve information available for the VLA is time-dependent (on timescales of months to years, at least for the higher frequencies), and CASA will automatically select the date-appropriate gain curve information. Note, however, that the time-dependence was poorly sampled prior to 2001, and so gain curve corrections prior to this time should be considered with caution.
The troposphere is not completely transparent. At high radio frequencies (>15 GHz), water vapor and molecular oxygen begin to have a substantial effect on radio observations. According to the physics of radiative transmission, the effect is threefold. First, radio waves from astronomical sources are absorbed (and therefore attenuated) before reaching the antenna. Second, since a good absorber is also a good emitter, significant noise-like power will be added to the overall system noise. Finally, the optical path length through the troposphere introduces a time-dependent phase error. In all cases, the effects become worse at lower elevations due to the increased air mass through which the antenna is looking. In CASA, the opacity correction described here compensates only for the first of these effects, tropospheric attenuation, using a plane-parallel approximation for the troposphere to estimate the elevation dependence.
To make opacity corrections in CASA, an estimate of the zenith opacity is required (see observatory-specific chapters for how to measure zenith opacity). This is then supplied to the caltype=’opac’ parameter in gencal which creates a calibration table with all the information. E.g. for data with two spectral windows, the inputs are like:
gencal(vis='dataset.ms',
caltable='opacity.cal',
caltype='opac',
spw='0,1',
parameter=[0.0399,0.037])
If you do not have an externally supplied value for opacity, for example from a VLA tip procedure, then you should either use an average value for the telescope, or leave it at zero and let your gain calibration compensate as best it can (e.g. that your calibrator is at the same elevation as your target at approximately the same time. As noted above, there are no facilities yet to estimate this from the data (e.g. by plotting Tsys vs. elevation).
Below, we give instructions for determining opacity for Jansky VLA data from weather statistics and VLA observations where tip-curve data is available. It is beyond the scope of this cookbook to provide information for other telescopes.
For the VLA site, weather statistics and/or seasonal models that average over many years of weather statistics prove to be reasonable good ways to estimate the opacity at the time of the observations. The task plotweather calculates the opacity as a mix of both actual weather data and seasonal model. It has the following inputs:
# plotweather :: Plot elements of the weather table; estimate opacity.
vis = '' # MS name
seasonal_weight = 0.5 # weight of the seasonal model
doPlot = True # set this to True to create a plot
The task plots the weather statistics if doPlot=T, like shown in Figure 4.3. The bottom panel displays the calculated opacities for the run as well as a seasonal model. The parameter seasonal_weight can be adjusted to calculate the opacities as a function of the weather data alone seasonal_weight=0, only the seasonal model seasonal_weight=1, or a mix of the two (values between 0 and 1). Calculated opacities are shown in the logger output, one for each spectral window. plotweather can also assign a python variable to a list of calculated opacities (one entry for each spw) when run as:
myTau = plotweather(vis='myvladata.ms')
In this example, myTau will be returned with a list of per-spw opacities, e.g. myTau=[0.02,0.03] and can later be used as input for gencal in caltype=’opac’ in the parameter setting, e.g.,
# opac for spws 0,1 in myTau
gencal(vis='myvladata.ms',caltype='opac',spw='0,1',parameter=myTau)
Note that it is important to explicitly specify the spws that are covered by the opacity values stored in myTau. For most modern VLA data there will be more than two spws, probably.
See § 4.3.6 for more information on use of gencal.
For VLA data, zenith opacity can be measured at the frequency and during the time observations are made using a VLA tipping scan in the observe file. Historical tipping data are available at:
http://www.vla.nrao.edu/astro/calib/tipper
Choose a year, and click Go to get a list of all tipping scans that have been made for that year.
If a tipping scan was made for your observation, then select the appropriate file. Go to the bottom of the page and click on the button that says Press here to continue.. The results of the tipping scan will be displayed. Go to the section called ’Overall Fit Summary’ to find the fit quality and the fitted zenith opacity in percent. If the zenith opacity is reported as 6%, then the actual zenith optical depth value is opacity=0.060 for gaincal and other calibration tasks.
If there were no tipping scans made for your observation, then look for others made in the same band around the same time and weather conditions. If nothing is available here, then at K and Q bands you might consider using an average value (e.g. 6% in reasonable weather). See the VLA memo
http://www.vla.nrao.edu/memos/test/232/232.pdf
for more on the atmospheric optical depth correction at the VLA, including plots of the seasonal variations.
When solving for visibility-plane calibration, CASA calibration applications compare the observed DATA column with the Fourier transform of calibrator model when it is provided (if no model is specified, a point source at the phase center is assumed).
The setjy task is used to set the proper flux density and attaches a model image (if specified) of the calibrator to the MS. For sources that are recognized flux calibrators (listed in Tables 4.1 and 4.2, see also § C.1), setjy can calculate the flux densities as a function of frequency (and time, for Solar System objects). Otherwise, the flux densities should be manually specified ( standard=’manual’). This mode has the following options:
standard = 'manual' # Flux density standard
fluxdensity = [1, 0, 0, 0] # Specified flux density [I,Q,U,V]; (-1 will lookup values)
spix = [] # Spectral index (including higher terms) of I fluxdensity
reffreq = '1GHz' # Reference frequency for spix
polindex = [] # Coefficients of an expansion of frequency-dependent linear polarization fraction expression
polangle = [] # Coefficients of an expansion of frequency-dependent polarization angle expression
rotmeas = 0.0 # Rotation measure (in
rad/m^2)
In the simplest version, the flux for Stokes I can be set via the fluxdensity subparameter as the first entry in a list. In the above example [1,0,0,0] setjy sets the flux to 1 Jy. Additional Stokes specifications can be set in remaining list members. A spectral index can be provided via the spix and reffreq parameters. Finally it is also possible to provide coefficients for the polarization fraction and angle as well as a rotation measure to define the model (polindex, polangle, rotmeas parameters).
For the VLA, the default source models are customarily point sources defined by the ’Baars’, ’Perley 90’, ’Perley-Taylor 99’, ’Perley-Butler 2010’, time-variable ’Perley-Butler 2013’,’Scaife-Heald 2012’, or ’Stevens-Reynolds 2016’ flux density scales (§ C.1.1; ’Perley-Butler 2013’ is the current standard by default), or point sources of unit flux density if the flux density is unknown. In fact, the model can be any image in Jy/pixel units (models typically generated by the clean task).
Optionally, the MODEL column can be filled with the Fourier transform of (option usescratch=T is setjy, ft, and clean). But for most Measurement Sets, the performance and data storage requirements are less demanding without the MODEL_DATA column.
The inputs for setjy are:
# setjy :: Fills the model column with the visibilities of a calibrator
vis = '' # Name of input visibility file
field = '' # Field name(s)
spw = 'all' # Spectral window identifier (list)
selectdata = True # Other data selection parameters
timerange = '' # Time range to operate on (for usescratch=T)
scan = '' # Scan number range (for usescaratch=T)
intent = '' # Observation intent
observation = '' # Observation ID range (for usescratch=T)
scalebychan = True # scale the flux density on a per channel basis or else on
# a per spw basis
standard = 'Perley-Butler 2013' # Flux density standard
model = '' # File location for field model
listmodels = False # List the available modimages for VLA calibrators or Tb
# models for Solar System objects
usescratch = False # Will create if necessary and use the MODEL_DATA
Table 4.1: Recognized Flux Density Calibrators. Note that the VLA uses J2000 calibrator names. CASA accepts all strings that contain the names below. E.g. ’PKS 1934-638’ will be recognized Standards are: (1) Perley-Butler 2010, (2) Perley-Butler 2013, (3) Perley-Taylor 99, (4) Perley-Taylor 95, (5) Perley 90, (6) Baars (Baars, J. W. M., et al. 1977, A&A, 61, 99); (7) Scaife-Heald 2012, (8) Stevens-Reynolds 2016; see § C.1.1 for details.
3C Name B1950 Name J2000 Name Alt. J2000 Name Standards 3C48 0134+329 0137+331 J0137+3309 1,3,4,5,6,7 3C123 0433+295 0437+296 J0437+2940 2 3C138 0518+165 0521+166 J0521+1638 1,3,4,5,6 3C147 0538+498 0542+498 J0542+4951 1,3,4,5,6,7 3C196 0809+483 0813+482 J0813+4813 1,2,7 3C286 1328+307 1331+305 J1331+3030 1,2,3,4,5,6,7 3C295 1409+524 1411+522 J1411+5212 1,2,3,4,5,6,7 – 1934-638 – J1939-6342 1,3,4,5,6,8 3C380 1828+487 1829+487 J1829+4845 7
1 Venus: model for ∼300 MHz to 350 GHz, no atmospheric lines (CO,H2O,HDO, etc.)
Planets Venus1, Mars2, Jupiter3, Uranus4, Neptune5 Moons Jupiter: Io, Europa, Ganymede, Callisto Saturn: Titan7 Asteroids Ceres, Pallas8, Vesta8, Juno8
2 Mars: tabulated as a function of time and frequency (30 - 1000GHz) based on Rudy et al. (1988), no atmospheric lines (CO, H20, H2O2, HDO, etc.)
3 Jupiter: model for 30-1020GHz, does not include synchrotron emission
4 Uranus: model for 60-1800GHz, contains no rings or synchrotron.
5 Neptune: model for 2-2000GHz, the broad CO absorption line is included, but contains no rings or synchrotron.
7 Titan: model for 53.3-1024.1GHz, include many spectral lines
8 not recommended (The temperature is not yet adjusted for varying distance from the Sun. The model data can be scaled after running setjy, but it is an involved process.)
Details are described in ALMA Memo 594 available on https://science.nrao.edu/facilities/alma/aboutALMA/Technology/ALMA_Memo_Series/alma594/abs594.
By default the setjy task will cycle through all fields spectral windows and channels, (one solution per spw with scalebychan = False) , setting the flux density either to 1 Jy (unpolarized), or if the source is recognized as one of the calibrators in the above table, to the flux density (assumed unpolarized) appropriate to the observing frequency. For example, to run setjy on a Measurement Set called data.ms:
setjy(vis='data.ms') # This will set all fields and spectral windows
Models of available calibrator sources can be listed by setting listmodels=True. setjy will then come up with all images that are in the paths where calibrator models for known telescopes are stored. It will also show all images in the working directory - any image there could potentially be a calibrator model. If the calibrator model is found by listmodels it can be used in the modimage parameter without a path.
The fluxdensity parameter can be used to specify the flux of the calibrator in all Stokes parameters. It it thus a list of values [I,Q,U,V], e.g. [’12Jy’,’13mJy’,’0Jy’,’0Jy’]. In addition, a spectral index can be specified via spix, a reference frequency reffreq (using the definition: S = fluxdensity× freq/reffreqspix), as well as a polarization index ( polindex), angle (polangle) and a rotation measure (rotmeas).
Most calibrator sources are based on radio emission from quasars and jets. The spectral indices of these sources are such that at (sub)mm wavelengths the majority of these sources become too weak and variable to be reliable flux estimators. Alternatives are thermal objects such as planets, moons, and asteroids. Those sources, however, are all Solar System objects, which implies that they move and may be (strongly) resolved. The standard=’Butler-JPL-Horizons 2010’ and the recommended standard=’Butler-JPL-Horizons 2012’ (for more information on the implemented models, see ALMA Memo 594 soon available on https://science.nrao.edu/facilities/alma/aboutALMA/Technology/ALMA_Memo_Series/alma594/abs594.) option of setjy includes flux density calibration using Solar System objects. For ’Butler-JPL-Horizons 2012’ CASA currently supports the objects listed in Table 4.2 to be applied to ALMA data. These names are recognized when they are used in the ’field’ parameter in setjy. In that case, setjy will obtain the geocentric distance and angular diameter at the time of the observation from a (JPL–Horizons) ephemeris and calculate model visibilities. Currently the objects are modeled as uniform temperature disks, but effects like primary beam attenuation and limb darkening will be accounted for soon. Note that this model may oversimplify the real structure, in particular asteroids.
An example, using data from the M99 tutorial in http://casaguides.nrao.edu/index.php?title=CARMA_spectral_line_mosaic_M99:
setjy(vis=’c0104I’, field=’MARS’, spw=’0 2’, standard=’Butler-JPL-Horizons 2012’)
Tip: Running casalog.filter(’INFO1’) before running setjy with a Solar System object may send the logger a reference to the temperature measurement. Use casalog.filter(’INFO’) to restore the normal logging level.
The source model will be attached to the MS and applied to all calibration steps when usescratch=False. usescratch=True fills the MODEL_DATA column with the Fourier transform of the model. As of CASA 3.4. we found that under some circumstances, creation of the MODEL column may prevent memory issues and if tasks fail, we recommend to set usescratch=True. Note that currently setjy will not transform a full-Stokes model image such that all polarizations are applied correctly. You need to use ft for this.
To limit this operation to certain fields and spectral windows, use the field and/or spw parameters, which take the usual data selection strings (§ 2.3). For example, to set the flux density of the first field (all spectral windows)
setjy(vis='data.ms',field='0')
or to set the flux density of the second field in spectral window 17
setjy(vis='data.ms',field='1',spw='17')
The full-polarization flux density (I,Q,U,V) may also be explicitly provided:
setjy(vis='data.ms',
field='1',spw='16', # Run setjy on field id 1, spw id 17
fluxdensity=[3.5,0.2,0.13,0.0]) # and set I,Q,U,V explicitly
ALERT: The apparent brightness of objects in the Solar System will vary with time because of the Earth’s varying distance to them, if nothing else. If the field index of a flux calibrator spans several days, you should run setjy more than once, limiting each run to a suitable timerange by using the timerange, scan, and/or observation selection parameters. Note that it is the field index that matters, not the name. Typically concat assigns moving objects a new field index for each observation, so usually it is not necessary to select a time range in setjy. However, it is worth checking with listobs, especially for planets.
For observations of solar system objects using the ’Butler-JPL-Horizons 2010’ and ’Butler-JPL-Horizons 2012’ models (§ 4.3.5) setjy will know and apply the flux distribution across the extended structure of the calibrators.
For other sources, namely VLA calibrator sources, a flux density calibrator can be resolved at the observing frequency and the point source model generated by setjy will not be appropriate. If available, a model image of the resolved source at the observing frequency may be used to generate the appropriate visibilities using the modimage parameter (or in older versions explicitly with the ft task). To use this, provide modimage with the path to the model image. Remember, if you just give the file name, it will assume that it is in the current working directory. Note also that setjy using a model image will only operate on that single source, thus you would run it multiple times (with different field settings) for different sources.
Otherwise, you may need to use the uvrange selection (§ 4.4.1.2) in the calibration solving tasks to exclude the baselines where the resolution effect is significant. There is not hard and fast rule for this, though you should consider this if your calibrator is shows a drop of more than 10% on the longest baselines (use plotxy, § 3.3.2, to look at this). You may need to do antenna selection also, if it is heavily resolved and there are few good baselines to the outer antennas. Note that uvrange may also be needed to exclude the short baselines on some calibrators that have extended flux not accounted for in the model. Note: the calibrator guides for the specific telescopes usually indicate appropriate min and max for uvrange. For example, see the VLA Calibration Manual at:
https://science.nrao.edu/facilities/vla/observing/callist
for details on the use of standard calibrators for the VLA.
Model images for some flux density calibrators are provided with CASA (note that the exactl location may change depending on your installation directory):
e.g., these are found in the data/nrao/VLA/CalModels sub-directory of the CASA installation. For example, just point to the repository copy, e.g.
modimage = '/usr/lib/casapy/data/nrao/VLA/CalModels/3C48_C.im'
or if you like, you can copy the ones you wish to use to your working directory.
The models available are:
3C138_L.im 3C147_L.im 3C286_L.im 3C48_L.im
3C138_C.im 3C147_C.im 3C286_C.im 3C48_C.im
3C138_X.im 3C147_X.im 3C286_X.im 3C48_X.im
3C138_U.im 3C147_U.im 3C286_U.im 3C48_U.im
3C138_K.im 3C147_K.im 3C286_K.im 3C48_K.im
3C138_Q.im 3C147_Q.im 3C286_Q.im 3C48_Q.im
(more calibrator models for the VLA are available at
https://science.nrao.edu/facilities/vla/data-processing/models
These are all un-reconvolved images of AIPS CC lists.
It is important that the model image not be one
convolved with a finite beam; it must have units of Jy/pixel (not
Jy/beam).
Note that setjy will rescale the flux in the models for known sources (e.g. those in Table 4.1) to match those it would have calculated. It will thus extrapolated the flux out of the frequency band of the model image to whatever spectral windows in the MS are specified (but will use the structure of the source in the model image).
ALERT: The reference position in the modimage is currently used by setjy when it does the Fourier transform, thus differences from the positions for the calibrator in the MS will show up as phase gradients in the uv-plane. If your model image position is significantly different but you don’t want this to affect your calibration, then you can doctor either the image header using imhead (§ 6.2) or in the MS (using the ms tool) as appropriate. In an upcoming release we will put in a toggle to use or ignore the position of the modimage. Note that this will not affect the flux scaling (only put in erroneous model phases); in any event small position differences, such as those arising by changing epoch from B1950 to J2000 using regridimage (§ 6.14), will be inconsequential to the calibration.
This illustrates the use of uvrange for a slightly resolved calibrator:
# Import the data
importvla(archivefiles='AS776_A031015.xp2', vis='ngc7538_XBAND.ms',
freqtol=10000000.0, bandname='X')
# Flag the ACs
flagautocorr('ngc7538_XBAND.ms')
# METHOD 1: Use point source model for 3C48, plus uvrange in solve
# Use point source model for 3C48
setjy(vis='ngc7538_XBAND.ms',field='0');
# Limit 3C48 (fieldid=0) solutions to uvrange = 0-40 klambda
gaincal(vis='ngc7538_XBAND.ms', caltable='cal.G', field='0',
solint=60.0, refant='10', selectdata=True, uvrange='0~40klambda',
append=False)
# Append phase-calibrator's solutions (no uvrange) to the same table
gaincal(vis='ngc7538_XBAND.ms', caltable='cal.G', field='2',
solint=60.0, refant='10', selectdata=True, uvrange='',
append=True)
# Fluxscale
fluxscale(vis='ngc7538_XBAND.ms', caltable='cal.G', reference=['0137+331'],
transfer=['2230+697'], fluxtable='cal.Gflx', append=False)
while the following illustrates the use of a model:
# METHOD 2: use a resolved model copied from the data repository
# for 3C48, and no uvrange
# (NB: detailed freq-dep flux scaling TBD)
# Copy the model image 3C48_X.im to the working directory first!
setjy(vis='ngc7538_XBAND.ms', field='0', modimage='3C48_X.im')
# Solutions on both calibrators with no uvrange
gaincal(vis='ngc7538_XBAND.ms', caltable='cal.G2', field='0,2',
solint=60.0, refant='10',
append=False)
# Fluxscale
fluxscale(vis='ngc7538_XBAND.ms', caltable='cal.G2', reference=['0137+331'],
transfer=['2230+697'], fluxtable='cal.G2flx', append=False)
# Both methods give 2230 flux densities ~0.7 Jy, in good agreement with
# AIPS
The gencal task provides a means of specifying antenna-based calibration values manually. The values are put in designated tables and can be applied to the data on-the-fly in solving tasks and using applycal.
The gencal task has the inputs:
# gencal :: Specify Calibration Values of Various Types
vis = '' # Name of input visibility file
caltable = '' # The new/existing calibration table
caltype = 'tecim' # The calibration type: 'amp','ph','sbd','mbd','antpos','an
# tposvla','tsys','evlagain','opac','gc','gceff','eff','te
# cim'
infile = '' # Input ancilliary file
spw = 'all' # Calibration spw(s) selection
antenna = '1' # Calibration antenna(s) selection
pol = '' # Calibration polarizations(s) selection
parameter = [] # The calibration values
Current antenna-based gencal options (caltype) are:
The calibration parameter specifications cannot be time-variable in the present implementation (though some of them will introduce implicit time-dependence upon evaluation in the apply). Calibration values can be assigned to each spw, antenna and pol selection, where applicable. The list of calibration values specified in parameter must conform to the range of spectral windows, antennas, and polarizations specified in spw, antenna and pol, with the values specified in order of the specified polarizations (fastest), antennas, and spectral windows (slowest). If any of spw, antenna, or pol are left unspecified (empty strings), the values specified in parameter will be assumed applicable to all values of the unspecified data axes. The output caltable will otherwise assume nominal calibration values for unspecified spectral windows, antennas, and polarizations. Note that antenna position corrections formally do not have spectral-window or polarization dependence; such specifications should not be used with ’antpos’.
The same caltable can be specified for multiple runs of gencal, in which case the specified parameters will be incorporated cumulatively. E.g., amplitude parameters (caltype=’amp’) multiply and phase-like parameters (’ph’, ’sbd’,’mbd’,’antpos’) add. Parameters for ’amp’ and ’ph’ corrections can be incorporated into the same caltable (in separate runs), but each of the other types require their own unique caltable. A mechanism for specifying manual corrections via a text file will be provided in the future.
Two kinds of delay corrections are supported. For caltype=’sbd’, the specified delays (in nanoseconds) will be applied locally to each spectral window, referring the derived phase corrections to each spectral window’s reference frequency (where the phase correction will be zero). The phases in each spectral window will nominally be flattened, but any phase offsets between spectral windows will remain. (These can be corrected using caltype=’phase’, or via ordinary spectral-window-dependent phase calibration.) For caltype=’mbd’, the evaluated phase corrections are referred to zero frequency. This causes a correction that is coherent over many spectral windows. If the data are already coherent over many spectral windows and share a common multi-band delay (e.g., VLA data, per baseband), caltype=’mbd’ corrections will maintain this coherence and flatten the frequency-dependent phase. Using caltype=’sbd’ in this instance will introduce phase offsets among spectral windows that reflect the multi-band delay.
For antenna position corrections (caltype=’antpos’), the antenna position offsets are specified in the ITRF frame. If the antenna field is left empty, gencal will try to look up the appropriate antenna position offsets at the time of the observation from the VLA baseline webpage http://www.vla.nrao.edu/astro/archive/baselines/.
For VLA position corrections in the VLA-centric frame, use caltype=’antposvla’, and gencal will rotate them to ITRF before storing them in the output caltable.
The sign and scale convention for gencal corrections (indeed for all CASA caltables) is such that the specified parameters (and as stored in caltables) are the factors that corrupt ideal data to yield the observed data. Thus, when applied to correct the data, their effective inverse will automatically be taken. I.e., amplitude factors will be divided into the data on correction. Phase-like parameters adopt the convention that the complex factor for the second antenna in the baseline is conjugated, and then both antenna factors are divided into the data on correction. (These conventions differ from AIPS in that multiplying correction factors are stored in AIPS calibration tables; however, the phase convention ends up being the same since AIPS conjugates the complex factor for the first antenna in the baseline.)
The following series of examples illustrate the use of gencal.
For the dataset ’data.ms’, the following sequence of gencal runs introduces, into a single caltable (’test.G’), (1) an antenna-based amplitude scale correction of 3.0 for all polarizations, antennas, and spectral windows, (2) phase corrections for all spectral windows and polarizations of 45 and 120 degrees to antennas EA03 and EA04, respectively, (3) phase corrections for all spectral windows of 63 and -34 in R (only) for antennas EA05 and EA06, respectively, and (4) phase corrections for all spectral windows of 14, -23, -130, and 145 degrees for antenna/polarizations EA09/R, EA09/L, EA10/R, and EA10/L, respectively:
gencal(vis='data.ms',caltable='test.G',caltype='amp', \
spw='',antenna='',pol='', \
parameter=[3])
gencal(vis='data.ms',caltable='test.G',caltype='ph', \
spw='',antenna='EA03,EA04',pol='', \
parameter=[45,120])
gencal(vis='data.ms',caltable='test.G',caltype='ph', \
spw='',antenna='EA05,EA06',pol='R', \
parameter=[63,-34])
gencal(vis='data.ms',caltable='test.G',caltype='ph', \
spw='',antenna='EA09,EA10',pol='R,L', \
parameter=[14,-23,-130,145])
In the following example, delay corrections in both polarizations will be adjusted for antenna EA09 by 14 nsec in spw 2 and -130 nsec in spw 3, and for antenna EA10 by -23 nsec in spw 2 and 145 nsec in spw 3:
gencal(vis='test.ms',caltable='test.sbd',caltype='sbd', \
spw='2,3',antenna='EA09,EA10',pol='', \
parameter=[14,-23,-130,145])
In the following example, antenna position corrections in meters (in
ITRF) for antenna EA09 (dBx=0.01, dBy=0.02, dBz=0.03) and for antenna
EA10 (dBx=-0.03, dBy=-0.01, dBz=-0.02) are introduced. Note that
three parameters are required for each antenna. The antenna offsets
can be obtained for the ’Jansky VLA/ old VLA Baseline Corrections’ web page:
://www.vla.nrao.edu/astro/archive/baselines. The table
given on this webpage has a format like:
; 2010 BASELINE CORRECTIONS IN METERS
;ANT
;MOVED OBSDATE Put_In_ MC(IAT) ANT PAD Bx By Bz
;
JAN27 FEB12 FEB21 01:57 11 E04 0.0000 0.0000 0.0000
JAN27 FEB12 FEB21 01:57 26 W03 -0.0170 0.0204 0.0041
MAR24 MAR25 MAR26 18:28 17 W07 -0.0061 -0.0069 -0.0055
APR21 MAY02 MAY04 23:25 12 E08 -0.0072 0.0045 -0.0017
If your observations fall in between the ’Antenna Moved’ and ’Put_In_’ dates of a given antenna, you may choose to apply the offsets in that table; the ’Put_In_’ time stamp marks the date where the more accurate solution was introduced in the data stream directly and no correction is required anymore. In gencal the offsets will be inserted as:
gencal(vis='test.ms',caltable='test.antpos',caltype='antpos', \
antenna='EA09,EA10', \
parameter=[0.01,0.02,0.03, -0.03,-0.01,-0.02])
In the following example, antenna position corrections (in the traditional VLA-centric frame) will be introduced in meters for antenna EA09 (dBx=0.01, dBy=0.02, dBz=0.03) and for antenna EA10 (dBx=-0.03, dBy=-0.01, dBz=-0.02) These offsets will be rotated to the ITRF frame before storing them in the caltable.
gencal(vis='test.ms',caltable='test.antposvla',caltype='antposvla', \
antenna='EA09,EA10', \
parameter=[0.01,0.02,0.03, -0.03,-0.01,-0.02])
gencal is also the task to generate gaincurve, antenna efficiency, and opacity tables. The first two items can be determined together with caltype=’gceff’ and the latter with caltype=’opac’. These tables are treated just like any other calibration table and will be carried through the calibration steps. This method replaces the older method where ’gaincurve’ and ’opacity’ keywords were present in calibration tasks such as gaincal, bandpass, or applycal.
Noise diodes in the Jansky VLA antennas can be used to pre-calibrate the data. The diodes follow an ON-OFF cycle and the power for both states is measured and recorded. This is called the ’VLA switched power’ calibration. To apply the switched power data, one needs to create a calibration table with gencal using caltype=’evlagain’, like
gencal(vis='test.ms',caltable='VLAswitchedpower.cal',caltype='evlagain')
For ALMA the calibration of system temperature is done via hot loads and the data recorded similar to the VLA in the Measurement Set (ALMA will provide Measurement Sets where these data are available. To derive the calibration table from it, use caltype=’tsys’:
gencal(vis='test.ms',caltable='ALMAtsys.cal',caltype='tsys')
This calibration tables created for ALMA or VLA are then carried along all further calibration steps in the gaintable parameter.
# wvrgcal :: Generate a gain table based on Water Vapour Radiometer data
vis = '' # Name of input visibility file
caltable = '' # Name of output gain
calibration table
toffset = 0 # Time offset (sec) between
interferometric and WVR data
segsource = True # Do a new coefficient
calculation for each source
tie = [] # Prioritise tieing the phase
of these sources as well as possible (requires
# segsource=True)
sourceflag = [] # Regard the WVR data for
these source(s) as bad and do not produce corrections for it
# (requires segsource=True)
disperse = False # Apply correction for dispersion
wvrflag = [''] # Regard the WVR data for these
antenna(s) as bad and replace its data with
# interpolated values from
neighbouring antennas
statfield = '' # Compute the statistics (Phase
RMS, Disc) on this field only
statsource = '' # Compute the statistics (Phase
RMS, Disc) on this source only
smooth = '' # Smooth calibration solution
on the given timescale
scale = 1.0 # Scale the entire phase
correction by this factor
spw = [] # List of the spectral window
IDs for which solutions should be saved into the
# caltable
wvrspw = [] # List of the spectral window
IDs from which the WVR data should be taken
reversespw = '' # Reverse the sign of the
correction for the listed SPWs (only needed for early ALMA
# data before Cycle 0)
cont = False # Estimate the continuum (e.g.,
due to clouds) (experimental)
maxdistm = 500.0 # maximum distance (m) of an
antenna used for interpolation for a flagged antenna
minnumants = 2 # minimum number of near
antennas (up to 3) required for interpolation
mingoodfrac = 0.8 # If the fraction of unflagged
data for an antenna is below this value (0. to 1.), the
# antenna is flagged.
usefieldtab = False # derive the antenna AZ/EL
values from the FIELD rather than the POINTING table
refant = [''] # use the WVR data from this
antenna for calculating the dT/dL parameters (can give
# ranked list)
The task wvrgcal generates a gain table based on Water Vapor Radiometer (WVR) data and is used for ALMA data reduction. It is an interface to the executable ‘‘wvrgcal’’, which is part of the CASA 4.6 distribution and can also be called from outside CASA. The wvrgcal software is based on the libair and libbnmin libraries which were developed by Bojan Nikolic at the University of Cambridge as part of EU FP6 ALMA Enhancement program.
CASA 4.6 contains version 2.1 of wvrgcal. The algorithmic core of wvrgcal is described in three ALMA memos (number 587, 588, and 593) which describe the algorithms implemented in the software. They can be found at at the ALMA Memo Series1. Since wvrgcal version 2.0 (2013), maintenance and extensions of this software are done by ESO (contact: D. Petry).
Briefly, wvrgcal follows a Bayesian approach to calculate the coefficients that convert the outputs of the ALMA 183 GHz water-vapor radiometers (mounted on each antenna) into estimates of path fluctuations which can then be used to correct the observed interferometric visibilities.
The CASA task interface to wvrgcal follows closely the interface of the shell executable at the same time staying within the CASA task parameter conventions.
In ALMA data, the WVR measurements belonging to a given observation are contained in the ASDM for that observation. After conversion to an MS using importasdm, the WVR information can be found in separate spectral windows. As of April 2016, it is still one single spectral window for all WVRs, however, the ID of the spectral window may vary between datasets. wvrgcal identifies the SPW autonomously but it can also be specified via the parameter "wvrspw" (see below). The specified spectral window(s) must be present in the MS for wvrgcal to work.
The various features of wvrgcal are then controlled by a number of task parameters (see the list above). They have default values which will work for ALMA data. An example for a typical wvrgcal call can be found in the ALMA CASA guide for the NGC 3256 analysis:
wvrgcal(vis='uid___A002_X1d54a1_X5.ms',
caltable='cal-wvr-uid___A002_X1d54a1_X5.W', toffset=-1,
segsource=True, tie=["Titan,1037-295,NGC3256"], statsource="1037-295",
spw=[17,19,21,23])
Here, vis is the name of input visibility file (which as mentioned above also contains the WVR data in spectral window 0) and caltable is the name of the output gain calibration table.
toffset is the known time offset in seconds between the WVR measurements and the visibility integrations they are valid for. For ALMA, this offset is presently -1 s (which is also the default value).
The parameter segsource (segregate source) controls whether separate coefficients are calculated for each source. The default value True is the recommended one for ALMA. When segsource is True, the subparameter tie is available. It permits to form groups of sources for which common coefficients are calculated as well as possible. The tie parameter ensures best possible phase transfer between a group of sources. In general it is recommended to tie together all of the sources in a single Science Goal (in ALMA speak) and their phase calibrator(s). The recommended maximum angular distance up to which two sources can be tied is 15∘.
The parameter statsource controls for which sources statistics are calculated and displayed in the logger. This has no influence on the generated calibration table.
Via the parameter spw, one can control for which of the input spectral windows wvrgcal will calculate phase corrections and store them in the output calibration table. By default, solutions for all spectral windows are written except for the ones containing WVR data.
wvrgcal respects the flags in the Main and ANTENNA table of the MS. The parameter mingoodfrac lets the user set a requirement on the minimum fraction of good measurements for accepting the WVR data from an antenna. If antennas are flagged, their WVR solution is interpolated from the three nearest neighboring antennas. This process can be controlled with the new parameters maxdistm and minnumants. The former sets the maximum distance an antenna used for interpolation may have from the flagged one. And minnumants sets how many near antennas there have to be for interpolation to take place.
For more details on the WVR Phase correction, see also the the ALMA Memo “Quality Control of WVR Phase Correction Based on Differences Between WVR Channels” by B. Nikolic, R. C. Bolton & J. S. Richer2, see also ALMA memo #5933.
wvrgcal writes out a variety of information to the logger, including various statistical measures of the performance. This allows the user to judge whether WVR correction is appropriate for the ms, to check whether any antennas have problematic WVR values, and to examine the predicted performance of the WVR correction when applied.
For each set of correction coefficients which are calculated (the number of coefficient sets are controlled by the parameters nsol, segsource and tie), the wvrgcal output to the logger first of all shows the time sample, the individual temperatures of each of the four WVR channels and the elevation of the source in question at that time. For each of these coefficient sets, it then gives the evidence of the bayesian parameter estimation, the calculated precipitable water vapor (PWV) and its error in mm, and the correction coefficients found for each WVR channel (dTdL).
The output then shows the statistical information about the observation. First of all it gives the start and end times for the parts of the observation used to calculate these statistics (controlled by segsource). It then shows a break down for each of the antennas in the data set. This gives the antenna name and number; whether or not it has a WVR (column WVR); whether or not it has been flagged (column FLAG); the RMS of the path length variation with time towards that antenna (column RMS); and the discrepancy between the RMS path length calculated separately for different WVR channels (column Disc.). These values allow the user to see if an individual WVR appears to have been suffering from problems during the observation, and to flag that antenna using wvrflag if necessary.
This discrepancy value, Disc., can in addition be used as a simple diagnostic tool to evaluate whether or not the WVR correction caltable created by wvrgcal should be applied. In the event of the WVR observations being contaminated by strong cloud emission in the atmosphere, the attempt by wvrgcal to fit the water vapor line may not be successful, and applying the produced calibration table can in extreme cases reduce the quality of the data. However, these weather conditions should identified by a high value in the discrepancy column produced when running wvrgcal.
Discrepancy values of greater than a 1000 microns usually indicate strong cloud contamination of the WVR data, and the output calibration table should probably not be applied. If the values are between 100 and 1000 microns, then the user should manually examine the phases before and after applying the caltable to decide if WVR correction is appropriate. Work is underway at ALMA to provide additional routines to at least partially remove the cloud component from the WVR data before calculating phase corrections. CASA 4.7 will contain a first tested version of such a tool.
After the antenna-by-antenna statistics, the output then displays some estimates of the performance of the wvrgcal correction. These are the thermal contribution from the water vapor to the path fluctuations per antenna (in microns), the largest path fluctuation found on a baseline (in microns), and the expected error on the path length calculated for each baseline due to the error in the coefficients (in microns).
The information about antenna pointing direction is by default taken from the POINTING table. Should this table not be present for some reason, the user can instead switch to determining the antenna positions from the phase directions in the FIELD table (under the assumption that all antennas were pointing ideally). The switch is performed by setting the parameter usefieldtab to True.
By default wvrgcal puts solutions for all spectral windows of the MS into the output calibration table. Since usually only the spectral windows are of interest in which the science target and the calibrators were observed, it is not necessary to store solutions for other spectral windows.
The spectral windows for which solutions are stored can be selected
with the parameter spw, e.g. spw = [17,19,21,23] will
make wvrgcal write only solutions for spectral windows 17, 19, 21,
and 23.
W.r.t to the input WVR spectral windows, wvrgcal will by default regard all windows with 4 channels as WVR data. In typical ALMA data there is only one such spectral window in each ASDM. This may change in the future. In any case, the input WVR spectral window(s) can be selected with the optional parameter wvrspw. The syntax is the same as for the parameter spw above.
CASA 4.3 introduces initial support for on-axis ionospheric corrections, using time- and direction-dependent total electron content (TEC) information obtained from the internet. The correction includes the dispersive delay (∝ ν−1) delay and Faraday rotation (∝ ν−2) terms. These corrections are most relevant at observing frequencies less than ∼5 GHz. When relevant, the ionosphere correction table should be generated at the beginning of a reduction along with other calibration priors (antenna position errors, gain curve, opacity, etc.), and carried through all subsequent calibration steps. Formally, the idea is that the ionospheric effects (as a function of time and on-axis direction) will be nominally accounted for by this calibration table, and thus not spuriously leak into gain and bandpass solves, etc. In practice, the quality of the ionospheric correction is limited by the relatively sparse sampling (in time and direction) of the available TEC information. Especially active ionospheric conditions may not be corrected very well. Also, direction-dependent ionosphere corrections are not yet supported. (Various improvements are under study for future releases.)
To generate the ionosphere correction table, first import a helper function from the casapy recipes repository:
from recipes import tec_maps
Then, generate a TEC surface image:
tec_maps.create(vis='mydata.ms',doplot=T,imname='iono')
This function goes to the web to obtain TEC information for the observing date and location, and generates a time-dependent CASA image containing this information. The string specified for imname is used as a prefix for two output images, with suffixes .IGS_TEC.im (the actual TEC image) and .IGS_RMS_TEC.im (a TEC error image). If imname is unspecified, the MS name (from vis) will be used as the prefix.
The quality of the retrieved TEC information improves with time after the observing date, becoming optimal 1-2 weeks later. Both images can be viewed as a movie in the CASA viewer. If doplot=T, the function will also produce a plot of the TEC as a function of time in a vertical direction over the observatory.
Finally, to generate the ionosphere correction caltable, pass the .IGS_TEC.im image into gencal, using caltype=’tecim’:
gencal(vis='mydata.ms',caltable='tec.cal',caltype='tecim',infile='iono.IGS_TEC.im')
This iterates through the dataset and samples the zenith angle-dependent projected line-of-sight TEC for all times in the observation, storing the result in a standard CASA caltable. Plotting this caltable will show how the TEC varies between observing directions for different fields and times, in particular how it changes as zenith angle changes, and including the nominal difference between science targets and calibrators.
This caltable should then be used as a prior in all subsequent calibration solves, and included in the final applycal.
A few warnings:
Special thanks are due to Jason Kooi (UIowa) for his contributions to ionospheric corrections in CASA.
Other a priori calibrations will be added to the calibrater (cb) tool in the near future. These will include instrumental line-length corrections, etc. Where appropriate, solving capabilities for these effects will also be added.
The gaincal, bandpass, polcal, and blcal tasks actually solve for the unknown calibration parameters from the visibility data obtained on calibrator sources, placing the results in a calibration table. They take as input an MS, and a number of parameters that specify any prior calibration or previous calibration tables to pre-apply before computing the solution. These are placed in the proper sequence of the Measurement Equation automatically.
We first discuss the parameters that are in common between many of the calibration tasks. Then we describe each solver in turn.
There are a number of parameters that are in common between the calibration “solver” tasks. These also appear in some of the other calibration manipulation and application tasks.
The input Measurement Set and output table are controlled by the following parameters:
vis = '' # Name of input visibility file
caltable = '' # Name of output calibration table
The MS name is input in vis. If it is highlighted red in the inputs (§ 1.4.5.4) then it does not exist, and the task will not execute. Check the name and path in this case.
The output table name is placed in caltable. Be sure to give a unique name to the output table, or be careful. If the table exists, then what happens next will depend on the task and the values of other parameters (e.g. § 4.4.1.6). The task may not execute giving a warning that the table already exists, or will go ahead and overwrite the solutions in that table, or append them. Be careful.
Selection is controlled by the parameters:
field = '' # field names or index of calibrators: ''==>all
spw = '' # spectral window:channels: ''==>all
intent = '' # Select observing intent
selectdata = False # Other data selection parameters
Field and spectral window selection are so often used, that we have made these standard parameters field and spw respectively. intent is the scan intent that was specified when the observations were set up. They typically describe what was intended with a specific scan, i.e. a flux or phase calibration, a bandpass, a pointing, an observation of your target, or something else or a combination. The format for the scan intents of your observations are listed in the logger when you run listobs. Minimum matching with wildcards will work, like ’*BANDPASS*’. This is especially useful when multiple intents are attached to scans. Finally, observation is an identifier to distinguish between different observing runs, mainly used for ALMA.
The selectdata parameter expands as usual, uncovering other selection sub-parameters:
selectdata = True # data selection parameters
field = '' # field names or field index
# numbers (blank for all)
spw = '' # spectral windows:channels (blank for all)
timerange = '' # time range (blank for all)
uvrange = '' # uv range (blank for all)
antenna = '' # antenna/baselines (blank for all)
scan = '' # scan numbers (blank for all)
correlation = '' # correlations (blank for all)
array = '' # (sub)array numbers (blank for all)
observation = '' # Select by observation ID(s)
msselect = '' # MS selection (blank for all)
Note that if selectdata=False these parameters are not used when the task is executed, even if set underneath.
The most common selectdata parameter to use is uvrange, which can be used to exclude longer baselines if the calibrator is resolved, or short baselines of the calibrator contains extended flux not accounted for in the model (e.g. § 4.3.5.1).
See § 2.3 for more on the selection parameters.
These parameters control the on-the-fly application of various calibration or effect-based corrections prior to the solving process.
The parang parameter turns on the application of the antenna-based parallactic angle correction (’P’) in the measurement equation. This is necessary for polarization calibration and imaging, or for cases where the parallactic angles are different for geographically spaced antennas and it is desired that the ordinary gain calibration not absorb the inter-antenna parallactic angle phase. When dealing with only the parallel-hand data (e.g. RR, LL, XX, YY), and an unpolarized calibrator model for a co-located array (e.g. the VLA or ALMA), you can set parang=False and save some computational effort. Otherwise, set parang=True to apply this correction.
See § 4.3 for more on Prior Calibration, including how to invoke gaincurve and opacity correction using gencal.
Calibration tables that have already been determined can also be applied before solving for the new table:
docallib = False # Use traditional cal apply parameters
gaintable = [] # Gain calibration table(s) to apply on the fly
gainfield = [] # Select a subset of calibrators from gaintable(s)
interp = [] # Interpolation mode (in time) to use for each gaintable
spwmap = [] # Spectral windows combinations to form for gaintable(s)
This is controlled by the gaintable parameter, which takes a string or list of strings giving one or more calibration tables to pre-apply. For example,
gaintable = ['ngc5921.bcal','ngc5921.gcal']
specifies two tables, in this case bandpass and gain calibration tables respectively.
The other parameters key off gaintable, taking single values or lists, with an entry for each table in gaintable. The order is given by that in gaintable.
The gainfield parameter specifies which fields from the respective gaintable to select for apply. This is a list, with each entry a string or list of strings. The default ’’ for an entry means to use all in that table. For example,
gaintable = ['ngc5921.bcal','ngc5921.gcal']
gainfield = [ '1331+305', ['1331+305','1445+099'] ]
or using indices
gainfield = [ '0', ['0','1'] ]
to specify the field ’1331+305’ from the table ’ngc5921.bcal’ and fields ’1331+305’ and ’1445+099’ from the second table ’ngc5921.gcal’. We could also have wildcarded the selection, e.g.
gainfield = [ '0', '*' ]
taking all fields from the second table. And of course we could have used the default
gainfield = [ '0', '' ]
or even
gainfield = [ '0' ]
which is to take all for the second table in gaintable. In addition, gainfield can be specified by
gainfield = [ 'nearest' ]
which selects the calibrator that is the spatially closest (in sky coordinates) to each of the selected MS fields specified in the field parameter. This is particularly useful for running applycal with a number of different sources to be calibrated in a single run.
The interp parameter chooses the interpolation scheme to be used when pre-applying the solution in the tables. Interpolation in both time and frequency (for channel-dependent calibrations) are supported. The choices are currently ’nearest’ and ’linear’, and ’nearest’, ’linear’, cubic, and spline for frequency-dependent interpolation. Frequency-dependent interpolation is only relevant for channel-dependent calibration tables (like bandpasses) that are undersampled in frequency relative to the data.
For each gaintable, specify the interpolation style in quotes, with the frequency-dependent interpolation style specified after a comma, if relevant.
If the uncalibrated phase is changing rapidly, a ’nearest’ interpolation is not desirable. Usually, interp=’linear’ is the best choice. For example,
gaintable=['gain','bandpass']
interp = [ 'nearest', 'linear,spline' ]
uses nearest “interpolation” on the first table, and linear (in time) and spline (in freq) on the second.
The spwmap parameter sets the spectral window combinations to form for the gaintable(s). This is a list, or a list of lists, of integers giving the spw IDs to map. There is one list for each table in gaintable, with an entry for each ID in the MS. For example,
spwmap=[0,0,1,1] # apply from spw=0 to 0,1 and 1 to 2,3
for an MS with spw=0,1,2,3. For multiple gaintable, use lists of lists, e.g.
spwmap=[ [0,0,1,1], [0,1,0,1] ] # 2nd table spw=0 to 0,2 and 1 to 1,3
The parameters controlling common aspects of the solution are:
solint = 'inf' # Solution interval: egs. 'inf', '60s' (see help)
combine = 'scan' # Data axes which to combine for solve (obs, scan,
# spw, and/or field)
preavg = -1.0 # Pre-averaging interval (sec) (rarely needed)
refant = '' # Reference antenna name(s)
minblperant = 4 # Minimum baselines _per antenna_ required for solve
minsnr = 3.0 # Reject solutions below this SNR
The time and frequency (if relevant) solution interval is given by solint. Optionally a frequency interval for each solution can be added after a comma, e.g. solint=’60s,300Hz’. Time units are in seconds unless specified differently. Frequency units can be either channels or Hz and only make sense for bandpass of frequency dependent polarization calibration. The special values ’inf’ and -1 specify an “infinite” solution interval encompassing the entire dataset, while ’int’ or zero specify a solution every integration. You can use time quanta in the string, e.g. solint=’1min’ and solint=’60s’ both specify solution intervals of one minute. Note that ’m’ is a unit of distance (meters); ’min’ must be used to specify minutes. The solint parameter interacts with combine to determine whether the solutions cross scan or field boundaries.
The parameter controlling the scope of the solution is combine. For the default combine=’’ solutions will break at obsId, scan, field, and spw boundaries. Specification of any of these in combine will extend the solutions over the boundaries (up to the solint). For example, combine=’spw’ will combine spectral windows together for solving, while combine=’scan’ will cross scans, and combine=’obs,scan’ will use data across different observation IDs and scans (usually, obsIds consist of many scans, so it is not meaningful to combine obsIds without also combining scans). Thus, to do scan-based solutions (single solution for each scan), set
solint = 'inf'
combine = ''
while
solint = 'inf'
combine = 'scan'
will make a single solution for the entire dataset (for a given field and spw).
solint = 'inf,30ch'
will calculate a bandpass solution for each scan, averaging over 30 channels.
You can specify multiple choices for combination:
combine = 'scan,spw'
for example.
Alert: Care should be exercised when using combine=’spw’ in cases where multiple groups of concurrent spectral windows are observed as a function of time. Currently, only one aggregate spectral window can be generated in a single calibration solve execution, and the meta-information for this spectral window is calculated from all selected MS spectral windows. To avoid incorrect calibration meta-information, each spectral window group should be calibrated independently (also without using append=True. Additional flexibility in this area will be supported in a future version.
The reference antenna is specified by the refant parameter. A list of antennas can be provided to this parameter and if the first antenna is not present in the data, the next antenna in the list will be used, etc. It is useful to “lock” the solutions with time, effectively rotating (after solving) the phase of the gain solutions for all antennas such that the reference antenna’s phase is constant at zero. If the selected antenna drops out, another antenna will be selected for ongoing consistency in time (at its “current” value) until the refant returns, usually at a new value (not zero), which will be kept fixed thenceforth. You can also run without a reference antenna, but in this case the solutions will formally float with time; in practice, the first antenna will be approximately constant near zero phase. It is usually prudent to select an antenna in the center of the array that is known to be particularly stable, as any gain jumps or wanders in the refant will be transferred to the other antenna solutions. Also, it is best to choose a reference antenna that never drops out.
Setting a preavg time (only needed in polcal) will let you average data over periods shorter than the solution interval first before solving on longer timescales.
The minimum signal-to-noise ratio allowed for an acceptable solution is specified in the minsnr parameter. Default is minsnr=3. The minblperant parameter sets the minimum number of baselines to other antennas that must be preset for each antenna to be included in a solution. This enables control of the constraints that a solution will require for each antenna.
The following parameters control some things that happen after solutions are obtained:
solnorm = False # Normalize solution amplitudes post-solve.
append = False # Append solutions to (existing) table.
# False will overwrite.
The solnorm parameter toggles on the option to normalize the solution after the solutions are obtained. The exact effect of this depends upon the type of solution. Not all tasks use this parameter.
One should be aware when using solnorm that if this is done in the last stage of a chain of calibration, then the part of the calibration that is “normalized” away will be lost. It is best to use this in early stages (for example in a first bandpass calibration) so that later stages (such as final gain calibration) can absorb the lost normalization scaling. It is not strictly necessary to use solnorm=True at all, but is sometimes helpful if you want to have a normalized bandpass for example.
The append parameter, if set to True, will append the solutions from this run to existing solutions in caltable. Of course, this only matters if the table already exists. If append=False and caltable exists, it will overwrite.
The append parameter should be used with care, especially when also using combine in non-trivial ways. E.g., calibration solves will currently (CASA 4.5) refuse to append incongruent aggregate spectral windows (e.g., observations with more than one group of concurrent spectral windows). This limitation arises from difficulty determining the appropriate spectral window fan-out on apply, and will be relaxed in a future version.
For channelized data, it is usually desirable to solve for the gain variations in frequency as well as in time. Variation in frequency arises as a result of non-uniform filter passbands or other frequency-dependent effects in signal transmission. It is usually the case that these frequency-dependent effects vary on timescales much longer than the time-dependent effects handled by the gain types ’G’ and ’T’. Thus, it makes sense to solve for them as a separate term: ’B’, using the bandpass task.
The inputs to bandpass are:
# bandpass :: Calculates a bandpass calibration solution
vis = '' # Name of input visibility file
caltable = '' # Name of output gain calibration table
field = '' # Select field using field id(s) or field name(s)
spw = '' # Select spectral window/channels
intent = '' # Select observing intent
selectdata = True # Other data selection parameters
timerange = '' # Select data based on time range
uvrange = '' # Select data within uvrange (default units meters)
antenna = '' # Select data based on antenna/baseline
scan = '' # Scan number range
observation = '' # Select by observation ID(s)
msselect = '' # Optional complex data selection (ignore for now)
solint = 'inf' # Solution interval in time[,freq]
combine = 'scan' # Data axes which to combine
# for solve (obs, scan, spw, and/or field)
refant = '' # Reference antenna name(s)
minblperant = 4 # Minimum baselines _per
# antenna_ required for solve
minsnr = 3.0 # Reject solutions below this
# SNR (only applies for bandtype = B)
solnorm = False # Normalize average solution amplitudes to 1.0
bandtype = 'B' # Type of bandpass solution (B or BPOLY)
fillgaps = 0 # Fill flagged solution channels by interpolation
smodel = [] # Point source Stokes parameters for source model.
append = False # Append solutions to the (existing) table
docallib = False # Use callib or traditional cal apply parameters
gaintable = [] # Gain calibration table(s) to apply on the fly
gainfield = [] # Select a subset of calibrators from gaintable(s)
interp = [] # Interpolation mode (in
# time) to use for each gaintable
spwmap = [] # Spectral windows
# combinations to form for gaintable(s)
parang = False # Apply parallactic angle correction
Many of these parameters are in common with the other calibration tasks and are described above in § 4.4.1.
The bandtype parameter selects the type of solution used for the bandpass. The choices are ’B’ and ’BPOLY’. The former solves for a complex gain in each channel in the selected part of the MS. See § 4.4.2.2 for more on ’B’. The latter uses a polynomial as a function of channel to fit the bandpass, and expands further to reveal a number of sub-parameters See § 4.4.2.3 for more on ’BPOLY’.
It is usually best to solve for the bandpass in channel data before solving for the gain as a function of time. However, if the gains of the bandpass calibrator observations are fluctuating over the timerange of those observations, then it can be helpful to first solve for the gains of that source with gaincal , and input these to bandpass via gaintable. See more below on this strategy.
We now describe the issue of bandpass normalization, followed by a description of the options bandtype=’B’ and bandtype=’BPOLY’.
The solnorm parameter (§ 4.4.1.6) deserves more explanation in the context of the bandpass. Most users are used to seeing a normalized bandpass, where the mean amplitude is unity and fiducial phase is zero. The toggle solnorm=True allows this. However, the parts of the bandpass solution normalized away will be still left in the data, and thus you should not use solnorm=True if the bandpass calibration is the end of your calibration sequence (e.g. you have already done all the gain calibration you want to). Note that setting solnorm=True will NOT rescale any previous calibration tables that the user may have supplied in gaintable.
You can safely use solnorm=True if you do the bandpass first (perhaps after a throw-away initial gain calibration) as we suggest above in § 4.2, as later gain calibration stages will deal with this remaining calibration term. This does have the benefit of isolating the overall (channel independent) gains to the following gaincal stage. It is also recommended for the case where you have multiple scans on possibly different bandpass calibrators. It may also be preferred when applying the bandpass before doing gaincal and then fluxscale (§ 4.4.4), as significant variation of bandpass among antennas could otherwise enter the gain solution and make (probably subtle) adjustments to the flux scale.
We finally note that solnorm=False at the bandpass step in the calibration chain will still in the end produce the correct results. It only means that there will be a part of what we usually think of the gain calibration inside the bandpass solution, particularly if bandpass is run as the first step.
Calibration type ’B’ differs from ’G’ only in that it is determined for each channel in each spectral window. It is possible to solve for it as a function of time, but it is most efficient to keep the ’B’ solving timescale as long as possible, and use ’G’ or ’T’ for frequency-independent rapid time-scale variations.
The ’B’ solutions are limited by the signal-to-noise ratio available per channel, which may be quite small. It is therefore important that the data be coherent over the time-range of the ’B’ solutions. As a result, ’B’ solutions are almost always preceded by an initial ’G’ or ’T’ solve using gaincal (§ 4.4.3). In turn, if the ’B’ solution improves the frequency domain coherence significantly, a ’G’ or ’T’ solution following it will be better than the original.
For example, to solve for a ’B’ bandpass using a single short scan on the calibrator, then
default('bandpass')
vis = 'n5921.ms'
caltable = 'n5921.bcal'
gaintable = '' # No gain tables yet
gainfield = ''
interp = ''
field = '0' # Calibrator 1331+305 = 3C286 (FIELD_ID 0)
spw = '' # all channels
selectdata = False # No other selection
bandtype = 'B' # standard time-binned B (rather than BPOLY)
solint = 'inf' # set solution interval arbitrarily long
refant = '15' # ref antenna 15 (=VLA:N2) (ID 14)
bandpass()
On the other hand, we might have a number of scans on the bandpass calibrator spread over time, but we want a single bandpass solution. In this case, we could solve for and then pre-apply an initial gain calibration, and let the bandpass solution cross scans:
gaintable = 'n5921.init.gcal' # Our previously determined G table
gainfield = '0'
interp = 'linear' # Do linear interpolation
solint = 'inf' # One interval over dataset
combine = 'scan' # Solution crosses scans
Note that we obtained a bandpass solution for all channels in the MS. If explicit channel selection is desired, for example some channels are useless and can be avoided entirely (e.g. edge channels or those dominated by Gibbs ringing), then spw can be set to select only these channels, e.g.
spw = '0:4~59' # channels 4-59 of spw 0
This is not so critical for ’B’ solutions as for ’BPOLY’, as each channel is solved for independently, and poor solutions at edges can be ignored.
If you have multiple time solutions, then these will be applied using whatever time interpolation scheme is specified in later tasks.
The combine parameter (§ 4.4.1.5) can be used to combine data across spectral windows, scans, and fields.
For some observations, it may be the case that the SNR per channel is insufficient to obtain a usable per-channel ’B’ solution. In this case it is desirable to solve instead for a best-fit functional form for each antenna using the bandtype=’BPOLY’ solver. The ’BPOLY’ solver naturally enough fits (Chebychev) polynomials to the amplitude and phase of the calibrator visibilities as a function of frequency. Unlike ordinary ’B’, a single common ’BPOLY’ solution will be determined for all spectral windows specified (or implicit) in the selection. As such, it is usually most meaningful to select individual spectral windows for ’BPOLY’ solves, unless groups of adjacent spectral windows are known a priori to share a single continuous bandpass response over their combined frequency range (e.g., PdBI data).
The ’BPOLY’ solver requires a number of unique sub-parameters:
bandtype = 'BPOLY' # Type of bandpass solution (B or BPOLY)
degamp = 3 # Polynomial degree for BPOLY amplitude solution
degphase = 3 # Polynomial degree for BPOLY phase solution
visnorm = False # Normalize data prior to BPOLY solution
maskcenter = 0 # Number of channels in BPOLY to avoid in center of band
maskedge = 0 # Percent of channels in BPOLY to avoid at each band edge
The degamp and degphase parameters indicate the polynomial degree desired for the amplitude and phase solutions. The maskcenter parameter is used to indicate the number of channels in the center of the band to avoid passing to the solution (e.g., to avoid Gibbs ringing in central channels for PdBI data). The maskedge drops beginning and end channels. The visnorm parameter turns on normalization before the solution is obtained (rather than after for solnorm).
The combine parameter (§ 4.4.1.5) can be used to combine data across spectral windows, scans, and fields.
Note that bandpass will allow you to use multiple fields, and can determine a single solution for all specified fields using combine=’field’. If you want to use more than one field in the solution it is prudent to use an initial gaincal using proper flux densities for all sources (not just 1Jy) and use this table as an input to bandpass because in general the phase towards two (widely separated) sources will not be sufficiently similar to combine them, and you want the same amplitude scale. If you do not include amplitude in the initial gaincal, you probably want to set visnorm=True also to take out the amplitude normalization change. Note also in the case of multiple fields, that the ’BPOLY’ solution will be labeled with the field ID of the first field used in the ’BPOLY’ solution, so if for example you point plotcal at the name or ID of one of the other fields used in the solution, plotcal does not plot.
For example, to solve for a ’BPOLY’ (5th order in amplitude, 7th order in phase), using data from field 2, with G corrections pre-applied:
bandpass(vis='data.ms', # input data set
caltable='cal.BPOLY', #
spw='0:2~56', # Use channels 3-57 (avoid end channels)
field='0', # Select bandpass calibrator (field 0)
bandtype='BPOLY', # Select bandpass polynomials
degamp=5, # 5th order amp
degphase=7, # 7th order phase
gaintable='cal.G', # Pre-apply gain solutions derived previously
refant='14') #
The bandpass calibrator can have a spectral slope that will change the spectral properties of the solutions. If the slope is significant, the best way is to model the slope and store that model in the bandpass calibrator MS. To do so, go through the normal steps of bandpass and the gaincal runs on the bandpass and flux calibrators, followed by setjy of the flux calibrator. The next step would be to use fluxscale on the bandpass calibrator to derive the slope of it. fluxscale can store this information in a python dictionary which is subsequently fed into a second setjy run, this time using the bandpass calibrator as the source and the derived slope (the python dictionary) as input. This step will create a source model with the correct overall spectral slope for the bandpass. Finally, rerun bandpass and all other calibration steps again, making use of the newly created internal bandpass model.
The fundamental calibration to be done on your interferometer data is to calibrate the antenna-based gains as a function of time. Some of these calibrations are known beforehand (“a priori”) and others must be determined from observations of calibrators, or from observations of the target itself (“self-calibration”).
It is best to have determined a (constant or slowly-varying) “bandpass” from the frequency channels by solving for the bandpass (see above). Thus, the bandpass calibration table would be input to gaincal via the gaintable parameter (see below).
The gaincal task has the following inputs:
# gaincal :: Determine temporal gains from calibrator observations
vis = '' # Name of input visibility file
caltable = '' # Name of output gain calibration table
field = '' # Select field using field id(s) or field name(s)
spw = '' # Select spectral window/channels
intent = '' # Select observing intent
selectdata = True # Other data selection parameters
timerange = '' # Select data based on time range
uvrange = '' # Select data within uvrange (default units meters)
antenna = '' # Select data based on antenna/baseline
scan = '' # Scan number range
observation = '' # Select by observation ID(s)
msselect = '' # Optional complex data selection (ignore for now)
solint = 'inf' # Solution interval: egs. 'inf', '60s' (see help)
combine = '' # Data axes which to combine
# for solve (obs, scan, spw, and/or field)
preavg = -1.0 # Pre-averaging interval (sec) (rarely needed)
refant = '' # Reference antenna name(s)
minblperant = 4 # Minimum baselines _per antenna_ required for solve
minsnr = 3.0 # Reject solutions below this SNR
solnorm = False # Normalize average solution
# amplitudes to 1.0 (G, T only)
gaintype = 'G' # Type of gain solution (G,T,GSPLINE,K,KCROSS)
smodel = [] # Point source Stokes parameters for source model.
calmode = 'ap' # Type of solution: ('ap', 'p', 'a')
append = False # Append solutions to the (existing) table
docallib = False # Use callib or traditional cal apply parameters
gaintable = [] # Gain calibration table(s) to apply on the fly
gainfield = [] # Select a subset of calibrators from gaintable(s)
interp = [] # Temporal interpolation for
# each gaintable (=linear)
spwmap = [] # Spectral windows
# combinations to form for gaintable(s)
parang = False # Apply parallactic angle correction on the fly
Data selection is done through the standard field, spw, intent, and selectdata expandable sub-parameters (see § 2.3). The bulk of the other parameters are the standard solver parameters. See § 4.4.1 above for a description of these.
The gaintype parameter selects the type of gain solution to compute. The choices are ’T’, ’G’, and ’GSPLINE’. The ’G’ and ’T’ options solve for independent complex gains in each solution interval (classic AIPS style), with ’T’ enforcing a single polarization-independent gain for each co-polar correlation (e.g. RR and LL, or XX and YY) and ’G’ having independent gains for these. See § 4.4.3.1 for a more detailed description of ’G’ solutions, and § 4.4.3.2 for more on ’T’. The ’GSPLINE’ fits cubic splines to the gain as a function of time. See § 4.4.3.3 for more on this option.
Systematic time-dependent complex gain errors are almost always the dominant calibration effect, and a solution for them is almost always necessary before proceeding with any other calibration. Traditionally, this calibration type has been a catch-all for a variety of similar effects, including: the relative amplitude and phase gain for each antenna, phase and amplitude drifts in the electronics of each antenna, amplitude response as a function of elevation (gain curve), and tropospheric amplitude and phase effects. In CASA, it is possible to handle many of these effects separately, as available information and circumstances warrant, but it is still possible to solve for the net effect using calibration type G.
Generally speaking, type G can represent any per-spectral window multiplicative polarization- and time-dependent complex gain effect downstream of the polarizers. (Polarization- and time- independent effects upstream of the polarizers may also be treated implicitly with G.) Multi-channel data (per spectral window) will be averaged in frequency before solving (use calibration type B to solve for frequency-dependent effects within each spectral window).
To solve for G on, say, fields 1 & 2, on a 90s timescale, and do so relative to gaincurve corrections:
gaincal('data.ms',
caltable='cal.G', # Write solutions to disk file 'cal.G'
field='0,1', # Restrict field selection
solint=90.0, # Solve for phase and amp on a 90s timescale
gaintable=['cal.gc'] # a gain curve table from gencal
refant='3') #
plotcal('cal.G','amp') # Inspect solutions
These G solution will be referenced to antenna 4. Choose a well-behaved antenna that is located near the center of the array and is ever-present for the reference antenna. For non-polarization datasets, reference antennas need not be specified although you can if you want. If no reference antenna is specified, an effective phase reference that is an average over the data will be calculated and used. For data that requires polarization calibration, you must choose a reference antenna that has a constant phase difference between the right and left polarizations (e.g. no phase jumps or drifts). If no reference antenna (or a poor one) is specified, the phase reference may have jumps in the R–L phase, and the resulting polarization angle response will vary during the observation, thus corrupting the polarization imaging.
To apply this solution, along with the gain curve correction, to the calibrators (fields 0,1) and the target source (field 2):
applycal('data.ms',
field='0,1,2', # Restrict field selection (cals + src)
gaintable=['cal.gc','cal.G']) # Apply gc and G solutions to correct data
The calibrated data is written to the CORRECTED_DATA column, with calwt=True by default. This parameter can also be a list of Boolean values for which each entry then controls the calculation of weights based on each individual input calibration table. calwt=False will recompute the weights form the SIGMA column, thus resetting the weights to their original value.
Alert: Current (as of February 2014) Jansky VLA data has no calibrated weights (unless they are computed from switched power calibration). To avoid trouble, calwt=False should be set for those data sets. Older, pre-upgrade VLA data should still be calibrated with calwt=True.
At high frequencies, it is often the case that the most rapid time-dependent gain errors are introduced by the troposphere, and are polarization-independent. It is therefore unnecessary to solve for separate time-dependent solutions for both polarizations, as is the case for ’G’. Calibration type ’T’ is available to calibrate such tropospheric effects, differing from ’G’ only in that a single common solution for both polarizations is determined. In cases where only one polarization is observed, type ’T’ is adequate to describe the time-dependent complex multiplicative gain calibration.
In the following example, we assume we have a ’G’ solution obtained on a longish timescale (longer than a few minutes, say), and we want a residual ’T’ solution to track the polarization-independent variations on a very short timescale:
gaincal('data.ms', # Visibility dataset
caltable='cal.T', # Specify output table name
gaintype='T', # Solve for T
field='0,1', # Restrict data selection to calibrators
solint=3.0, # Obtain solutions on a 3s timescale
gaintable='cal120.G') # Pre-apply prior G solution
For dual-polarization observations, it will always be necessary to obtain a ’G’ solution to account for differences and drifts between the polarizations (which traverse different electronics), but solutions for rapidly varying polarization-independent effects such as those introduced by the troposphere will be optimized by using ’T’. Note that ’T’ can be used in this way for self-calibration purposes, too.
At high radio frequencies, where tropospheric phase fluctuates rapidly, it is often the case that there is insufficient signal-to-noise ratio to obtain robust ’G’ or ’T’ solutions on timescales short enough to track the variation. In this case it is desirable to solve for a best-fit functional form for each antenna using the ’GSPLINE’ solver. This fits a time-series of cubic B-splines to the phase and/or amplitude of the calibrator visibilities.
The combine parameter (§ 4.4.1.5) can be used to combine data across spectral windows, scans, and fields. Note that if you want to use combine=’field’, then all fields used to obtain a ’GSPLINE’ amplitude solution must have models with accurate relative flux densities. Use of incorrect relative flux densities will introduce spurious variations in the ’GSPLINE’ amplitude solution.
The ’GSPLINE’ solver requires a number of unique additional parameters, compared to ordinary ’G’ and ’T’ solving. The sub-parameters are:
gaintype = 'GSPLINE' # Type of solution (G, T, or GSPLINE)
splinetime = 3600.0 # Spline (smooth) timescale (sec), default=1 hours
npointaver = 3 # Points to average for phase wrap (okay)
phasewrap = 180 # Wrap phase when greater than this (okay)
The duration of each spline segment is controlled by splinetime. The actual splinetime will be adjusted such that an integral number of equal-length spline segments will fit within the overall range of data.
Phase splines require that cycle ambiguities be resolved prior to the fit; this operation is controlled by npointaver and phasewrap. The npointaver parameter controls how many contiguous points in the time-series are used to predict the cycle ambiguity of the next point in the time-series, and phasewrap sets the threshold phase jump (in degrees) that would indicate a cycle slip. Large values of npointaver improve the SNR of the cycle estimate, but tend to frustrate ambiguity detection if the phase rates are large. The phasewrap parameter may be adjusted to influence when cycles are detected. Generally speaking, large values (>180∘) are useful when SNR is high and phase rates are low. Smaller values for phasewrap can force cycle slip detection when low SNR conspires to obscure the jump, but the algorithm becomes significantly less robust. More robust algorithms for phase-tracking are under development (including fringe-fitting).
For example, to solve for ’GSPLINE’ phase and amplitudes, with splines of duration 600 seconds,
gaincal('data.ms',
caltable='cal.spline.ap',
gaintype='GSPLINE' # Solve for GSPLINE
calmode='ap' # Solve for amp & phase
field='0,1', # Restrict data selection to calibrators
splinetime=600.) # Set spline timescale to 10min
ALERT’: The ’GSPLINE’ solutions cannot yet be used in fluxscale. You should do at least some ’G’ amplitude solutions to establish the flux scale, then do ’GSPLINE’ in phase before or after to fix up the short timescale variations. Note that the “phase tracking” algorithm in ’GSPLINE’ needs some improvement.
gaintype=’K’ solves for simple antenna-based delays via Fourier transforms of the spectra on baselines to the reference antenna. This is not a global fringe fit but will be useful for deriving delays from data of reasonable snr. If combine includes ’spw’, multi-band delays solved jointly from all selected spectral windows will be determined, and will be identified with the first spectral window id in the output caltable. When applying a multi-band delay table, spwmap is required to distribute the solutions to all spectral windows.
After solving for delays, a subsequent bandpass is recommended to describe higher-order channel-dependent variation in the phase (and amplitude).
gaintype=’KCROSS’ solves for a global cross-hand delay. Use parang=T and apply prior gain and bandpass solutions. Alert: Multi-band delays are not yet supported for KCROSS solutions.
The ’G’ or ’T’ solutions obtained from calibrators for which the flux density was unknown and assumed to be 1 Jansky are correct in a time- and antenna- relative sense, but are mis-scaled by a factor equal to the inverse of the square root of the true flux density. This scaling can be corrected by enforcing the constraint that mean gain amplitudes determined from calibrators of unknown flux density should be the same as determined from those with known flux densities. The fluxscale task exists for this purpose.
The inputs for fluxscale are:
# fluxscale :: Bootstrap the flux density scale from standard calibrators
vis = '' # Name of input visibility file (MS)
caltable = '' # Name of input calibration table
fluxtable = '' # Name of output, flux-scaled calibration table
reference = [''] # Reference field name(s) (transfer flux scale FROM)
transfer = [''] # Transfer field name(s) (transfer flux scale TO), '' ->
# all
listfile = '' # Name of listfile that contains the fit information.
# Default is (no file).
append = False # Append solutions?
refspwmap = [-1] # Scale across spectral window boundaries. See help
# fluxscale
gainthreshold = -1.0 # Threshold (% deviation from the median) on gain
# amplitudes to be used in the flux scale calculation
antenna = '' # antennas to include/exclude
incremental = False # incremental caltable
fitorder = 1 # order of spectral fitting
display = False # display some statistics of flux scaling
Before running fluxscale, one must have first run setjy for the reference sources and run a gaincal that includes reference and transfer fields. After running fluxscale the output fluxtable caltable will have been scaled such that the correct scaling will be applied to the transfer sources.
For example, given a ’G’ table, e.g. ’cal.G’, containing solutions for a flux density calibrator (in this case ’3C286’) and for one or more gain calibrator sources with unknown flux densities (in this example ’0234+285’ and ’0323+022’):
fluxscale(vis='data.ms',
caltable='cal.G', # Select input table
fluxtable= 'cal.Gflx', # Write scaled solutions to cal.Gflx
reference='3C286', # 3C286 = flux calibrator
transfer='0234+258, 0323+022') # Select calibrators to scale
The output table, ’cal.Gflx’, contains either the scaling factors alone (incremental=T) to be used alongside with the input gain table ’cal.G’, or a scaled version of the gain table (incremental=F), that replaces it for the execution of applycal.
Note that the assertion that the gain solutions are independent of the calibrator includes the assumption that the gain amplitudes are strictly not systematically time-dependent in any way. While synthesis antennas are designed as much as possible to achieve this goal, in practice, a number of effects conspire to frustrate it. When relevant, it is advisable to pre-apply gaincurve and opacity corrections when solving for the ’G’ solutions that will be flux-scaled (see § 4.3 and § 4.4.1.3). When the ’G’ solutions are essentially constant for each calibrator separately, the fluxscale operation is likely to be robust.
fluxscale will report the fluxes of each spw for each source. In addition, it will attempt a fit across the spws of each source and report a spectral index and curvature (S∝(ν/ν0)α+log(ν/ν0))). This information can be subsequently used to build up a model for the spectral slope of a calibrator with the setjy task if required.
The fluxscale task can be executed on either ’G’ or ’T’ solutions, but it should only be used on one of these types if solutions exist for both and one was solved relative to the other (use fluxscale only on the first of the two).
ALERT: The ’GSPLINE’ option is not yet supported in fluxscale (see § 4.4.3.3).
If the reference and transfer fields were observed in different spectral windows, the refspwmap parameter may be used to achieve the scaling calculation across spectral window boundaries.
The refspwmap parameter functions similarly to the standard spwmap parameter (§ 4.4.1.4), and takes a list of indices indicating the spectral window mapping for the reference fields, such that refspwmap[i]=j means that reference field amplitudes from spectral window j will be used for spectral window i.
Note: You should be careful when you have a dataset with spectral windows with different bandwidths, and you have observed the calibrators differently in the different spw. The flux-scaling will probably be different in windows with different bandwidths.
For example,
fluxscale(vis='data.ms',
caltable='cal.G', # Select input table
fluxtable= 'cal.Gflx', # Write scaled solutions to cal.Gflx
reference='3C286', # 3C286 = flux calibrator
transfer='0234+258,0323+022' # Select calibrators to scale
refspwmap=[0,0,0]) # Use spwid 0 scaling for spwids 1 & 2
will use spw=0 to scale the others, while in
fluxscale(vis='data.ms',
caltable='cal.G', # Select input table
fluxtable='cal.Gflx', # Write scaled solutions to cal.Gflx
reference='3C286', # 3C286 = flux calibrator,
transfer='0234+285, 0323+022', # select calibrators to scale,
refspwmap=[0,0,1,1]) # select spwids for scaling,
the reference amplitudes from spectral window 0 will be used for spectral windows 0 and 1 and reference amplitudes from spectral window 2 will be used for spectral windows 2 and 3.
If the flux density calibrator is resolved, the assumption that it is a point source will cause solutions on outlying antennas to be biased in amplitude. In turn, the fluxscale step will be biased on these antennas as well. In general, it is best to use model for the calibrator, but if such a model is not available, it is important to limit the solution on the flux density calibrator to only the subset of antennas that have baselines short enough that the point-source assumption is valid. This can be done by using antenna and uvrange selection when solving for the flux density calibrator. For example, if antennas 1 through 8 are the antennas among which the baselines are short enough that the point-source assumption is valid, and we want to be sure to limit the solutions to the use of baselines shorter than 15000 wavelengths, then we can assemble properly scaled solutions for the other calibrator as follows (note: specifying both an antenna and a uvrange constraint prevents inclusion of antennas with only a small number of baselines within the specified uvrange from being included in the solution; such antennas will have poorly constrained solutions):
As an example, we first solve for gain solutions for the flux density calibrator (3C286 observed in field 0) using a subset of antennas
gaincal(vis='data.ms',
caltable='cal.G', # write solutions to cal.G
field='0' # Select the flux density calibrator
selectdata=True, # Expand other selectors
antenna='0~7', # antennas 0-7,
uvrange='0~15klambda', # limit uvrange to 0-15klambda
solint=90) # on 90s timescales, write solutions
# to table called cal.G
Now solve for other calibrator (0234+285 in field 1) using all antennas (implicitly) and append these solutions to the same table
gaincal(vis='data.ms',
caltable='cal.G', # write solutions to cal.G
field='1',
solint=90,
append=T) # Set up to write to the same table
Finally, run fluxscale to adjust scaling
fluxscale(vis='data.ms',
caltable='cal.G', # Input table with unscaled cal solutions
fluxtable='cal.Gflx', # Write scaled solutions to cal.Gflx
reference='3C286', # Use 3c286 as ref with limited uvrange
transfer='0234+285') # Transfer scaling to 0234+285
The fluxscale calculation will be performed using only the antennas common to both fields, but the result will be applied to all antennas on the transfer field. Note that one can nominally get by only with the uvrange selection, but you may find that you get strange effects from some antennas only having visibilities to a subset of the baselines and thus causing problems in the solving.
Full support for instrumental polarization calibration for the circular feed basis (e.g., VLA) is provided in CASA. Support for the linear feed basis (e.g., ALMA) is now practical (as of v4.0) and is also described below. The linear feed basis treatment will continue to be expanded and streamlined for the v4.3 release.
The inputs to polcal are:
# polcal :: Determine instrumental polarization calibrations
vis = '' # Name of input visibility file
caltable = '' # Name of output gain calibration table
field = '' # Select field using field id(s) or field name(s)
spw = '' # Select spectral window/channels
intent = '' # Select observing intent
selectdata = True # Other data selection parameters
timerange = '' # Select data based on time range
uvrange = '' # Select data within uvrange (default units meters)
antenna = '' # Select data based on antenna/baseline
scan = '' # Scan number range
observation = '' # Select by observation ID(s)
msselect = '' # Optional complex data selection (ignore for now)
solint = 'inf' # Solution interval
combine = 'obs,scan' # Data axes which to combine
# for solve (obs, scan, spw, and/or field)
preavg = 300.0 # Pre-averaging interval (sec)
refant = '' # Reference antenna name(s)
minblperant = 4 # Minimum baselines _per antenna_ required for solve
minsnr = 3.0 # Reject solutions below this SNR
poltype = 'D+QU' # Type of instrumental
# polarization solution (see help)
smodel = [] # Point source Stokes parameters for source model.
append = False # Append solutions to the (existing) table
docallib = False # Use callib or traditional cal apply parameters
gaintable = [] # Gain calibration table(s) to apply
gainfield = [] # Select a subset of calibrators from gaintable(s)
interp = [] # Interpolation mode (in
# time) to use for each gaintable
spwmap = [] # Spectral windows
# combinations to form for gaintable(s)
The polcal task uses many of the standard calibration parameters as described above in § 4.4.1.
The key parameter controlling polcal is poltype. The choices are:
There are channelized solution modes for the above options. For example, substitute ’Df’ for ’D’ in the ’D*’ modes described above to get a channelized D-term solution; substitute ’Xf’ for ’X’ to get channelized position angle correction.
ALERT: polcal will obtain a separate D-term solution for each field supplied to it. This limitation will be relaxed in the future, enabling more sensitive solutions.
ALERT: This section concentrates on polarization calibration for the circular feed basis. It will be generalized to include the linear feed basis for the v4.3 release. See § 4.4.5.4 for the currently supported processing steps for the linear feed basis.
Fundamentally, with good ordinary gain (and bandpass, if relevant) calibration already in hand, good polarization calibration must deliver both the instrumental polarization and position angle calibration. An unpolarized source can deliver only the first of these, but does not require parallactic angle coverage. A polarized source can only deliver the position angle calibration also if its polarization is known a priori. Sources that are polarized, but with unknown polarization, must always be observed with sufficient parallactic angle coverage, where "sufficient" is determined by SNR and the details of the solving mode.
These principles are stated assuming the instrumental polarization solution is solved using the "linear approximation" where cross-terms in more than a single product of the instrumental or source polarizations are ignored in the Measurement Equation (see § E). A more general non-linearized solution, with sufficient SNR, may enable some relaxation of the requirements indicated here, and modes supporting such an approach are currently under development.
For instrumental polarization calibration, there are 3 types of calibrator choice:
| Cal Polarization | Parallactic Angles | model | polmode | Result |
| unpolarized | any | set Q=U=0 | ’D’ or ’Df’ | D-terms only |
| known non-zero | 2+ scans | set Q,U | ’D+X’ or ’Df+X’ | D-terms and PA |
| unknown | 2+ scans | ignored | ’D+QU’ or ’Df+QU’ | D-terms and source |
Note that the parallactic angle ranges spanned by the scans in the modes that require this should be large enough to give good separation between the components of the solution. In practice, 60∘ is a good target.
Each of these solutions should be followed with a ’X’ solution on a source with known polarization position angle (and correct Q+iU in the model). ALERT: polmode=’D+X’ will soon be enhanced to deliver this automatically.
The polcal task will solve for the ’D’ or ’X’ terms using the model visibilities that are in the model attached to the MS. Calibration of the parallel hands must have already been carried out using gaincal and/or bandpass in order to align the phases over time and frequency. This calibration must be supplied through the gaintable parameters, but any cal-tables to be used in polcal must agree (e.g. have been derived from) the data in the DATA column and the FT of the model. Thus, for example, one would not use the cal-table produced by fluxscale as the rescaled amplitudes would no longer agree with the contents of the model.
Be careful when using resolved calibrators for polarization calibration. A particular problem is if the structure in Q and U is offset from that in I. Use of a point model, or a resolved model for I but point models for Q and U, can lead to errors in the ’X’ calibration. Use of a uvrange will help here. The use of a full-Stokes model with the correct polarization is the only way to ensure a correct calibration if these offsets are large.
When your data has more than one channel per spectral window, it is important to note that the calibrator polarization estimate currently assumes the source polarization signal is coherent across each spectral window. In this case, it is important to be sure there is no large cross-hand delay still present in your data. Unless the online system has accounted for cross-hand delays (typically intended, but not always achieved), the gain and bandpass calibration will only correct for parallel-hand delay residuals since the two polarizations are referenced independently. Good gain and bandpass calibration will typically leave a single cross-hand delay (and phase) residual from the reference antenna. Plots of cross-hand phases as a function of frequency for a strongly polarized source (i.e., that dominates the instrumental polarization) will show the cross-hand delay as a phase slope with frequency. This slope will be the same magnitude on all baselines, but with different sign in the two cross-hand correlations. This cross-hand delay can be estimated using the gaintype=’KCROSS’ mode of gaincal (in this case, using the strongly polarized source 3C286):
default('gaincal')
vis = 'polcal_20080224.cband.all.ms'
caltable = 'polcal.xdelcal'
field = '3C286'
spw = ''
solint = 'inf'
combine = 'scan'
refant = 'VA15'
smodel = [1.0,0.11,0.0,0.0]
gaintype = 'KCROSS'
gaintable = ['polcal.gcal','polcal.bcal']
gaincal()
Note that smodel is used to specify that 3C286 is polarized; it is not important to specify this polarization stokes parameters correctly, as only the delay will be solved for (not any absolute position angle or amplitude scaling). The resulting solution should be carried forward and applied along with the gain (.gcal) and bandpass (.bcal) solutions in subsequent polarization calibration steps.
In the following example, we do a standard ’D+QU’ solution on the bright source BLLac (2202+422) which has been tracked through a range in parallactic angle:
default('polcal')
vis = 'polcal_20080224.cband.all.ms'
caltable = 'polcal.pcal'
field = '2202+422'
spw = ''
solint = 'inf'
combine = 'scan'
preavg = 300.0
refant = 'VA15'
minsnr = 3
poltype = 'D+QU'
gaintable = ['polcal.gcal','polcal.bcal','polcal.xdelcal]
gainfield = ['']
polcal()
This assumes setjy and gaincal have already been run. Note that the original gain-calibration table is used in gaintable so that what is in the model is in agreement with what is in the gaintable, rather than using the table resulting from fluxscale.
Now, we need to set the R-L phase using a scan on 3C48 (0137+331):
default('polcal')
vis = 'polcal_20080224.cband.all.ms'
caltable = 'polcal.polx'
field = '0137+331'
refant = 'VA15'
minsnr = 3
poltype = 'X'
smodel = [1.0,-0.0348,-0.0217,0.0] # the fractional Stokes for 3C48
gaintable = ['polcal.gcal','polcal.bcal','polcal.xdelcal','polcal.pcal']
polcal()
Note that the fractional polarization of 3C48 has been properly specified in smodel here.
If, on the other hand, we had a scan on an unpolarized bright source, for example 3C84 (0319+415), we could use this to calibrate the leakages:
default('polcal')
vis = 'polcal_20080224.cband.all.ms'
caltable = 'polcal.pcal'
field = '0319+415'
refant = 'VA15'
poltype = 'D'
gaintable = ['polcal.gcal','polcal.bcal','polcal.xdelcal]
polcal()
We would then do the ’X’ calibration as before (but using this D-table in gaintable).
CASA v4.0.0 introduces supports for instrumental polarization calibration for the linear feed basis at a level that is now practical for the general user. Some details remain to be implemented with full flexibility, and much of what follows will be streamlined for the v4.1 release.
Calibrating the instrumental polarization for the linear feed basis is somewhat more complicated than the circular feed basis because the polarization effects (source and instrument) appear in all four correlations at first or zeroth order (whereas for circular feeds, the polarization information only enters the parallel hand correlations at second order). As a result, e.g., the time-dependent gain calibration will be distorted by any non-zero source polarization, and some degree of iteration will be required to isolate the gain calibration if the source polarization is not initially known. These complications can actually be used to advantage in solving for the instrumental calibration; in can be shown, for example, that a significantly linearly polarized calibrator enables a better instrumental polarization solution than an unpolarized calibrator.
In the following example, we show the processing steps for calibrating the instrumental polarization using a strongly (> 5%) polarized point-source calibrator (which is also the time-dependent gain calibrator) that has been observed over a range of parallactic angle (a single scan is not sufficient). We assume that we have calibrated the gain, bandpass, and cross-hand delay as described above, and that the gain calibration (polcal.gcal) was obtained assuming the calibrator was unpolarized.
First, we import some utility functions from the CASA recipes area:
from recipes.almapolhelpers import *
Since the gain calibrator was assumed unpolarized, the time-dependent gain solutions contain information about the source polarization. This can be seen by plotting the amp vs. time for this table using poln=’/’. The antenna-based polarization amplitude ratios will reveal the sinusoidal (in parallactic angle) of the source polarization. Run a utility method (qufromgain()) to extract the apparent source polarization estimates for each spw:
qu=qufromgain('polcal.gcal')
The source polarization reported for all spws should be reasonably consistent. This estimate is not as good as can be obtained from the cross-hands (see below) since it relies on the gain amplitude polarization ratio being stable which may not be precisely true. However, this estimate will be useful in resolving an ambiguity that occurs in the cross-hand estimates.
Next we estimate both the XY-phase offset and source polarization from the cross-hands. The XY-phase offset is a spectral phase-only bandpass relating the X and Y systems of the reference antenna. The cross-hand delay solved for above represents a systematic component (linear phase in frequency). If the XY-phase is solved for in a channel-dependent manner (as below), it is strictly not necessary to have solved for the cross-hand delay above, but it does not hurt (at it allows reasonably coherent channel averages for data examination). The source polarization occurs in the cross-hands as a sinusoidal function of parallactic angle that is common to both cross-hands on all baselines (for a point-source). If the XY-phase bandpass is uniformly zero, then the source linear polarization function will occur entirely in the real part of the cross-hand visibilities. Non-zero XY-phase has the effect of rotating the source linear polarization signature partially into the imaginary part, where circular (and instrumental) polarization occur (cf. the circular feed basis where the cross-hand phase merely rotates the position angle of linear polarization). The following solve averages all baselines together and first solves for a channelized XY-phase (the slope of the source polarization function in the complex plane), then corrects the slope and solves for a channel-averaged source polarization. This calibration is obtained using gaintype=’XYf+QU’ in gaincal:
default('gaincal')
vis = 'polcal_linfeed.ms'
caltable = 'polcal.xy0amb' # possibly with 180deg ambiguity
field = '1' # the calibrator
solint = 'inf'
combine = 'scan'
preavg = 200.0 # minimal parang change
smodel = [1,0,1,0] # non-zero U assumed
gaintype = 'XYf+QU'
gaintable = ['polcal.gcal','polcal.bcal','polcal.xdelcal]
gaincal()
Note that we imply non-zero Stokes U in smodel; this is to enforce the assumption of non-zero source polarization signature in the cross-hands in the ratio of data and model. This solve will report the center-channel XY-phase and apparent Q,U for each spw. The Q,U results should be recognizable in comparison to that reported by qufromgain() above. However, since the XY-phase has a 180 degree ambiguity (you can rotate the source polarization signature to lie entirely in the visibility real part by rotating clockwise or counter-clockwise), some or all spw QU estimates may have the wrong sign. We correct this using the xyamb() utility method, using the qu obtained from qufromgain() above (which is not ambiguous):
S=xyamb(xy='polcal.xy0amb',qu=qu,xyout='polcal.xy0')
The python variable S now contains the mean source model (Stokes I = 1; fractional Q,U; V=0) that can be used in a revision of the gain calibration and instrumental polarization calibration.
Next we revise the gain calibration using the full polarization source model:
default('gaincal')
vis = 'polcal_linfeed.ms'
caltable = 'polcal.gcal1'
field = '1'
solint = 'int' # or whatever was used previously
smodel = S # obtained from xyamb
gaintype = 'G'
gaintable = ['polcal.bcal']
parang = T # so source poln properly rotated
gaincal()
Note that parang=T so that the supplied source linear polarization is properly rotated in the parallel-hand visibility model. This new gain solution can be plotted with poln=’/’ as above to show that the source polarization is no longer distorting it. Also, if qufromgain is run on this new gain table, the reported source polarization should be statistically indistinguishable from zero.
Finally, we can now solve for the instrumental polarization:
default('polcal')
vis = 'polcal_linfeed.ms'
caltable = 'polcal.dcal'
field = '1'
solint = 'inf'
combine = 'scan'
preavg = 200
poltype = 'Dflls' # freq-dep LLS solver
refant = '' # no reference antenna
smodel = S
gaintable = ['polcal.gcal1','polcal.bcal','polcal.xdelcal','polcal.xy0']
polcal()
Note that no reference antenna is used since this solve will produce an absolute instrumental polarization solution that is registered to the assumed source polarization (S) and prior calibrations. Applying a refant (referring all instrumental polarization terms to a reference antenna’s X feed, which would then be assumed perfect) would, in fact, discard valid information about the imperfections in the reference antenna’s X feed. (Had we used an unpolarized calibrator, we would not have a valid xy-phase solution, nor would we have had access to the absolute instrumental polarization solution demonstrated here.)
A few points:
A full processing example for linear feed basis polarimetry is under development and will be distributed with an upcoming CASA release.
You can use the blcal task to solve for baseline-dependent (non-closing) errors. WARNING: this is in general a very dangerous thing to do, since baseline-dependent errors once introduced are difficult to remove. You must be sure you have an excellent model for the source (better than the magnitude of the baseline-dependent errors).
The inputs are (note that blcal does not yet use the docallib parameter:
# blcal :: Calculate a baseline-based calibration solution (gain or bandpass)
vis = '' # Name of input visibility file
caltable = '' # Name of output gain calibration table
field = '' # Select field using field id(s) or field name(s)
spw = '' # Select spectral window/channels
intent = '' # Select observing intent
selectdata = False # Other data selection parameters
solint = 'inf' # Solution interval
combine = 'scan' # Data axes which to combine for solve (scan, spw,
# and/or field)
freqdep = False # Solve for frequency dependent solutions
calmode = 'ap' # Type of solution" ('ap', 'p', 'a')
solnorm = False # Normalize average solution amplitudes to 1.0
gaintable = [''] # Gain calibration table(s) to apply on the fly
gainfield = [''] # Select a subset of
# calibrators from gaintable(s)
interp = [''] # Interpolation mode (in
# time) to use for each gaintable
spwmap = [] # Spectral windows combinations to form for
# gaintable(s)
gaincurve = False # Apply internal VLA antenna
# gain curve correction
opacity = [] # Opacity correction to apply (nepers), per spw
parang = False # Apply parallactic angle correction
The freqdep parameter controls whether blcal solves for “gain” (freqdep=False) or “bandpass” (freqdep=True) style non-closing calibration.
Other parameters are the same as in other calibration tasks. These common calibration parameters are described in § 4.4.1.
At some point, the user should examine (plotting or listing) the calibration solutions. Calibration tables can also be manipulated in various ways, such as by interpolating between times (and sources), smoothing of solutions, and accumulating various separate calibrations into a single table.
The plotcal task is available for examining solutions of all of the basic solvable types (G, T, B, D, M, MF, K). The inputs are:
# plotcal :: An all-purpose plotter for calibration results:
caltable = '' # Name of input calibration table
xaxis = '' # Value to plot along x axis (time,chan,amp,phase,real,imag,snr)
yaxis = '' # Value to plot along y axis (amp,phase,real,imag,snr)
poln = '' # Polarization to plot (RL,R,L,XY,X,Y,/)
field = '' # Field names or index: ''=all, '3C286,P1321*', '0~3'
antenna = '' # Antenna selection. E.g., antenna='3~5'
spw = '' # Spectral window: ''=all, '0,1' means spw 0 and 1
timerange = '' # Time selection ''=all
subplot = 111 # Panel number on display screen (yxn)
overplot = False # Overplot solutions on existing display
clearpanel = 'Auto' # Specify if old plots are cleared or not
iteration = '' # Iterate on antenna,time,spw,field
plotrange = [] # plot axes ranges: [xmin,xmax,ymin,ymax]
showflags = False # If true, show flags
plotsymbol = '.' # pylab plot symbol
plotcolor = 'blue' # initial plotting color
markersize = 5.0 # size of plot symbols
fontsize = 10.0 # size of label font
showgui = True # Show plot on gui
figfile = '' # ''= no plot hardcopy, otherwise supply name
ALERT: Currently, plotcal needs to know the MS from which caltable was derived to get indexing information. It does this using the name stored inside the table, which does not include the full path, but assumes the MS is in the cwd. Thus if you are using a MS in a directory other than the current one, it will not find it. You need to change directories using cd in IPython (or os.chdir() inside a script) to the MS location.
The controls for the plotcal window are the same as for plotxy (see § 3.3.2.1).
The xaxis and yaxis plot options available are:
of the calibration solutions that are in the caltable. The xaxis choices also include ’time’ and ’channel’ which will be used as the sensible defaults (if xaxis=’’) for gain and bandpass solutions respectively.
The poln parameter determines what polarization or combination of polarization is being plotted. The poln=’RL’ plots both R and L polarizations on the same plot. The respective XY options do equivalent things. The poln=’/’ option plots amplitude ratios or phase differences between whatever polarizations are in the MS (R and L. or X and Y).
The field, spw, and antenna selection parameters are available to obtain plots of subsets of solutions. The syntax for selection is given in § 2.3.
The subplot parameter is particularly helpful in making multi-panel plots. The format is subplot=yxn where yxn is an integer with digit y representing the number of plots in the y-axis, digit x the number of panels along the x-axis, and digit n giving the location of the plot in the panel array (where n = 1, ..., xy, in order upper left to right, then down). See § 3.3.2.8 for more details on this option.
The iteration parameter allows you to select an identifier to iterate over when producing multi-panel plots. The choices for iteration are: ’antenna’, ’time’, ’spw’, ’field’. For example, if per-antenna solution plots are desired, use iteration=’antenna’. You can then use subplot to specify the number of plots to appear on each page. In this case, set the n to 1 for subplot=yxn. Use the Next button on the plotcal window to advance to the next set of plots. Note that if there is more than one timestamp in a ’B’ table, the user will be queried to interactively advance the plot to each timestamp, or if multiplot=True, the antennas plots will be cycled through for each timestamp in turn. Note that iteration can take more than one iteration choice (as a single string containing a comma-separated list of the options). ALERT: the iteration order is fixed (independent of the order specified in the iteration string), for example:
iteration = 'antenna, time, field'
iteration = 'time, antenna, field'
will both iterate over each field (fastest) then time (next) and antenna (slowest). The order is:
iteration = 'antenna, time, field, spw'
from the slowest (outer loop) to fastest (inner loop).
The markersize and fontsize parameters are especially helpful in making the dot and label sizes appropriate for the plot being made. The screen shots in this section used this feature to make the plots more readable in the cookbook. Adjusting the fontsize can be tricky on multi-panel plots, as the labels can run together if too large. You can also help yourself by manually resizing the Plotter window to get better aspect ratios on the plots.
ALERT: Unfortunately, plotcal has many of the same problems that plotxy does, as they use similar code underneath. An overhaul is underway, so stay tuned.
For example, to plot amplitude or phase as a function of time for ’G’ solutions (after rescaling by fluxscale can look like
default('plotcal')
fontsize = 14.0 # Make labels larger
markersize = 10.0 # Make dots bigger
caltable = 'ngc5921.usecase.fluxscale'
yaxis = 'amp'
subplot = 211
plotcal()
yaxis = 'phase'
subplot = 212
plotcal()
The results are shown in Figure 4.4. This makes use of the subplot option to make multi-panel displays.
Figure 4.4: Display of the amplitude (upper) and phase (lower) gain solutions for all antennas and polarizations in the ngc5921 post-fluxscale table.
Similarly, to plot amplitude or phase as a function of channel for ’B’ solutions for NGC5921:
default('plotcal')
fontsize = 14.0 # Make labels larger
markersize = 10.0 # Make dots bigger
caltable = 'ngc5921.usecase.bcal'
antenna = '1'
yaxis = 'amp'
subplot = 311
plotcal()
yaxis = 'phase'
subplot = 312
plotcal()
yaxis = 'snr'
subplot = 313
plotcal()
The results are shown in Figure 4.5. This stacks three panels with amplitude, phase, and signal-to-noise ratio. We have picked antenna=’1’ to show.
For example, to show 6 plots per page of ’B’ amplitudes on a 3 × 2 grid:
default('plotcal')
fontsize = 12.0 # Make labels just large enough
markersize = 10.0 # Make dots bigger
caltable = 'ngc5921.usecase.bcal'
yaxis = 'amp'
subplot = 231
iteration = 'antenna'
plotcal()
See Figure 4.6 for this example. This uses the iteration parameter.
Developed at the NAASC, this is a generic task to display CASA Tsys and bandpass solution tables with options to overlay them in various combinations, and/or with an atmospheric transmission or sky temperature model. It works with both the ’new’ (CASA 3.4) and ’old’ calibration table formats, and allows for mixed mode spws (e.g. TDM and FDM for ALMA). It uses the new msmd tool to access the information about an MS. This task is still being developed as new ALMA observing modes are commissioned. So if you encounter problems, please report them.
The parameters of plotbandpass are as follows:
# plotbandpass :: Makes detailed plots of Tsys and bandpass solutions.
caltable = '' # Input table name, either a bandpass solution or a Tsys
# solution
antenna = '' # A comma-delimited string list of antennas (either names
# or integer indices) for which to display solutions.
# Default = all antennas.
field = '' # A comma-delimited string list of fields (either names
# or integer indices) for which to display solutions.
# Default = all fields.
spw = '' # A comma-delimited string list of spws for which to
# display solutions. Default = all spws.
yaxis = 'amp' # The quantity to plot on the y-axis ("amp", "phase",
# "both", "tsys", append "db" for dB).
xaxis = 'chan' # The quantity to plot on the x-axis ("chan" or "freq").
figfile = '' # The name of the plot file to produce.
plotrange = [0, 0, 0, 0] # The axes limits to use [x0,x1,y0,y1].
caltable2 = '' # A second cal table, of type BPOLY or B, to overlay on a
# B table
overlay = '' # Show multiple solutions in same frame in different
# colors (time, antenna, spw, baseband, or time,antenna)
showflagged = False # Show the values of the solution, even if flagged
timeranges = '' # Show only these timeranges, the first timerange being 0
markersize = 3 # Size of points
interactive = True # if False, then run to completion automatically without
# pause
showpoints = 'auto' # Draw points for the data (default=F for amp, T for
# phase)
showlines = 'auto' # Draw lines connecting the data (default=T for amp, F
# for phase)
subplot = '22' # 11..81,22,32 or 42 for RowsxColumns (default=22), any
# 3rd digit is ignored
poln = '' # Polarizations to plot: "" = all, or
# "RR","RL","LR","LL","XX","XY","YX","YY","RR,LL","XX,YY
# "
showatm = False # Compute and overlay the atmospheric transmission curve
solutionTimeThresholdSeconds = 30.0 # Consider 2 solutions simultaneous if within this
# interval in seconds
debug = False # Print verbose messages for debugging purposes
vis = '' # name of the ms for this table, in case it does not
# match the string in the caltable
showtsky = False # Compute and overlay the sky temperature curve instead
# of transmission
channeldiff = False # Set to a value > 0 (sigma) to plot derivatives of the
# solutions
basebands = '' # A baseband number or list of baseband numbers for which
# to display solutions. Default = all.
showBasebandNumber = False # Put the baseband converter number (BBC_NO) in the title
# of each plot
scans = '' # A scan or list of scans for which to display solutions.
# Default = all. Does not work with overlay="time".
figfileSequential = False # naming scheme for pngs: False: name by spw/antenna
# (default), True: figfile.000.png, figfile.001.png,
# etc.
The listcal task will list the solutions in a specified calibration table.
The inputs are:
# listcal :: List data set summary in the logger:
vis = '' # Name of input visibility file (MS)
caltable = '' # Input calibration table to list
field = '' # Select data based on field name or index
antenna = '' # Select data based on antenna name or index
spw = '' # Spectral window, channel to list
listfile = '' # Disk file to write, else to terminal
pagerows = 50 # Rows listed per page
An example listing is:
Listing CalTable: jupiter6cm.usecase.split.ms.smoothcal2 (G Jones)
---------------------------------------------------------------
SpwId = 0, channel = 0.
Time Field Ant : Amp Phase Amp Phase
--------------------- ---------- -------- --------------- ---------------
1999/04/16/14:10:43.5 'JUPITER' '1' : 1.016 -11.5 1.016 -9.2
'2' : 1.013 -5.3 0.993 -3.1
'3' : 0.993 -0.8 0.990 -5.1
'4' : 0.997 -10.7 0.999 -8.3
'5' : 0.985 -2.7 0.988 -4.0
'6' : 1.005 -8.4 1.009 -5.3
'7' : 0.894 -8.7 0.897 -6.8
'8' : 1.001 -0.1 0.992 -0.7
'9' : 0.989 -12.4 0.992 -13.5
'10' : 1.000F -4.2F 1.000F -3.2F
'11' : 0.896 -0.0 0.890 -0.0
'12' : 0.996 -10.6 0.996 -4.2
'13' : 1.009 -8.4 1.011 -6.1
'14' : 0.993 -17.6 0.994 -16.1
'15' : 1.002 -0.8 1.002 -1.1
'16' : 1.010 -9.9 1.012 -8.6
'17' : 1.014 -8.0 1.017 -7.1
'18' : 0.998 -3.0 1.005 -1.0
'19' : 0.997 -39.1 0.994 -38.9
'20' : 0.984 -5.7 0.986 3.0
'21' : 1.000F -4.2F 1.000F -3.2F
'22' : 1.003 -11.8 1.004 -10.4
'23' : 1.007 -13.8 1.009 -11.7
'24' : 1.000F -4.2F 1.000F -3.2F
'25' : 1.000F -4.2F 1.000F -3.2F
'26' : 0.992 3.7 1.000 -0.2
'27' : 0.994 -5.6 0.991 -4.3
'28' : 0.993 -10.7 0.997 -3.8
The calstat task will print the statistics of solutions in a specified calibration table.
The inputs are:
# calstat :: Displays statistical information on a calibration table
caltable = '' # Name of input calibration table
axis = 'amp' # Which values to use
datacolumn = 'gain' # Which data column to use
useflags = True # Take flagging into account?
# (not implemented)
For example:
CASA <3>: calstat('ngc5921.demo.gcal',axis='amp',datacolumn='gain')
Out[3]:
{'GAIN': {'max': 1.6031942367553711,
'mean': 1.4448433067117419,
'medabsdevmed': 0.0086394548416137695,
'median': 1.5732669830322266,
'min': 0.99916577339172363,
'npts': 280.0,
'quartile': 0.020265340805053711,
'rms': 1.4650156497955322,
'stddev': 0.24271160321065546,
'sum': 404.55612587928772,
'sumsq': 600.95579999685287,
'var': 0.058908922333086665}}
CASA <4>: calstat('ngc5921.demo.gcal',axis='phase',datacolumn='gain')
Out[4]:
{'GAIN': {'max': 0.091214209794998169,
'mean': -0.015221830284565011,
'medabsdevmed': 0.012778861448168755,
'median': -0.012778861448168755,
'min': -0.15903720259666443,
'npts': 280.0,
'quartile': 0.02537553571164608,
'rms': 0.031241731718182564,
'stddev': 0.027331476552707856,
'sum': -4.2621124796782031,
'sumsq': 0.27329283416317834,
'var': 0.00074700961055121926}}
The statistics can be captured as return variables from the task:
CASA <7>: mystat = calstat('ngc5921.demo.gcal',axis='amp',datacolumn='gain')
CASA <8>: print 'Gain Amp = ',mystat['GAIN']['mean'],'+/-',mystat['GAIN']['stddev']
Gain Amp = 1.44484330671 +/- 0.242711603211
ALERT: This task is still under development and currently offers no selection (e.g. by antenna) for the statistics.
The smoothcal task will smooth calibration solutions (most usefully G or T) over a longer time interval to reduce noise and outliers. The inputs are:
# smoothcal :: Smooth calibration solution(s) derived from one or more sources:
vis = '' # Name of input visibility file
tablein = '' # Input calibration table
caltable = '' # Output calibration table
field = '' # Field name list
smoothtype = 'median' # Smoothing filter to use
smoothtime = 60.0 # Smoothing time (sec)
Note that if no caltable is specified as output, smoothcal will overwrite the input tablein calibration table.
The smoothing will use the smoothtime and smoothtype parameters to determine the new data points which will replace the previous points on the same time sampling grid as for the tablein solutions. The currently supported smoothtype options:
Note that smoothtime defines the width of the time window that is used for the smoothing.
ALERT: Note that smoothcal currently smooths by field and spw, and thus you cannot smooth solutions from different sources or bands together into one solution.
An example using the smoothcal task to smooth an existing table:
smoothcal('n4826_16apr.ms',
tablein='n4826_16apr.gcal',
caltable='n4826_16apr.smoothcal',
smoothtime=7200.,
smoothtype='mean')
# Plot up before and after tables
plotcal('n4826_16apr.gcal','','amp',antenna='1',subplot=211)
plotcal('n4826_16apr.smoothcal','','amp',antenna='1',subplot=212)
This example uses 2 hours (7200 sec) for the smoothing time and smoothtype=’mean’. The plotcal results are shown in Figure 4.7.
ALERT: The accum task is generally no longer recommended for most calibration scenarios. Please write to the NRAO CASA helpdesk if you need support using accum.
The accum task is used to interpolate calibration solutions onto a different time grid, and to accumulate incremental calibrations into a cumulative calibration table. The manual accumulation of calibration is rarely required and can usually be achieved implicitly simply by running applycal with all the calibration tables given as a list in the gaintable parameter (and using gainfield, spwmap, and interp appropriately. However, sometimes it is desirable to see the interpolated calibration prior to application, and this section describes how this can be done.
Its inputs are:
# accum :: Accumulate incremental calibration solutions
vis = '' # Name of input visibility file
tablein = '' # Input (cumulative) calibration table; use '' on first run
accumtime = 1.0 # Timescale on which to create cumulative table
incrtable = '' # Input incremental calibration table to add
caltable = '' # Output (cumulative) calibration table
field = '' # List of field names to process from tablein.
calfield = '' # List of field names to use from incrtable.
interp = 'linear' # Interpolation mode to use for resampling incrtable solutions
spwmap = [-1] # Spectral window combinations to apply
The mapping implied here is
tablein + incrtable => caltable
(mathematically the cal solutions are multiplied as complex numbers as per the Measurement Equation). The tablein is optional (see below). You must specify an incrtable and a caltable.
The tablein parameter is used to specify the existing cumulative calibration table to which an incremental table is to be applied. Initially, no such table exists, and if tablein=’’ then accumulate will generate one from scratch (on-the-fly), using the timescale (in seconds) specified by the sub-parameter accumtime. These nominal solutions will be unit-amplitude, zero-phase calibration, ready to be adjusted by accumulation according to the settings of other parameters. When accumtime is negative (the default), the table name specified in tablein must exist and will be used. If tablein is specified, then the entries in that table will be used.
The incrtable parameter is used to specify the incremental table that should be applied to tablein. The calibration type of incrtable sets the type assumed in the operation, so tablein (if specified) must be of the same type. If it is not, accum will exit with an error message. (Certain combinations of types and subtypes will be supported by accum in the future.)
The caltable parameter is used to specify the name of the output table to write. If un-specified (’’), then tablein will be overwritten. Use this feature with care, since an error here will require building up the cumulative table from the most recent distinct version (if any).
The field parameter specifies those field names in tablein to which the incremental solution should be applied. The solutions for other fields will be passed to caltable unaltered. If the cumulative table was created from scratch in this run of accumulate, then the solutions for these other fields will be unit-amplitude, zero-phase, as described above.
The calfield parameter is used to specify the fields to select from incrtable to use when applying to tablein. Together, use of field and calfield permit completely flexible combinations of calibration accumulation with respect to fields. Multiple runs of accum can be used to generate a single table with many combinations. In future, a ’self’ mode will be enabled that will simplify the accumulation of field-specific solutions.
The spwmap parameter gives the mapping of the spectral windows in the incrtable onto those in tablein and caltable. The syntax is described in § 4.4.1.4.
The interp parameter controls the method used for interpolation. The options are (currently): ’nearest’ and ’linear’ for time-dependent interpolation, and ’nearest’, ’linear’, cubic, and spline for (optional) frequency-dependent interpolation. These are described in § 4.4.1.4. For most purposes, the ’linear’ option should suffice.
We now describe the two uses of accum.
ALERT: The accum task is generally no longer recommended for most calibration scenarios. Please write to the NRAO CASA helpdesk if you need support using accum.
Calibration solutions (most notably G or T) can be interpolated onto the timestamps of the science target observations using accum.
The following example uses accum to interpolate an existing table onto a new time grid:
accum(vis='n4826_16apr.ms',
tablein='',
accumtime=20.0,
incrtable='n4826_16apr.gcal',
caltable='n4826_16apr.20s.gcal',
interp='linear',
spwmap=[0,1,1,1,1,1])
plotcal('n4826_16apr.gcal','','phase',antenna='1',subplot=211)
plotcal('n4826_16apr.20s.gcal','','phase',antenna='1',subplot=212)
See Figure 4.8 for the plotcal results. The data used in this example is BIMA data (single polarization YY) where the calibrators were observed in single continuum spectral windows (spw=’0,1’) and the target NGC4826 was observed in 64-channel line windows (spw=’2,3,4,5’). Thus, it is necessary to use spwmap=[0,1,1,1,1,1] to map the bandpass calibrator in spw=’0’ onto itself, and the phase calibrator in spw=’1’ onto the target source in spw=’2,3,4,5’.
It is occasionally desirable to solve for and apply calibration incrementally. This is the case when a calibration table of a certain type already exists (from a previous solve), a solution of the same type and incremental relative to the first is required, and it is not possible or convenient to recover the cumulative solution by a single solve.
Much of the time, it is, in fact, possible to recover the cumulative solution. This is because the equation describing the solution for the incremental solution (using the original solution), and that describing the solution for their product are fundamentally the same equation—the cumulative solution, if unique, must always be the same no matter what initial solution is. One circumstance where an incremental solution is necessary is the case of phase-only self-calibration relative to a full amplitude and phase calibration already obtained (from a different field).
For example, a phase-only ’G’ self-calibration on a target source may be desired to tweak the full amplitude and phase ’G’ calibration already obtained from a calibrator. The initial calibration (from the calibrator) contains amplitude information, and so must be carried forward, yet the phase-only solution itself cannot (by definition) recover this information, as a full amplitude and phase self-calibration would. In this case, the initial solution must be applied while solving for the phase-only solution, then the two solutions combined to form a cumulative calibration embodying the net effect of both. In terms of the Measurement Equation, the net calibration is the product of the initial and incremental solutions.
Cumulative calibration tables also provide a means of generating carefully interpolated calibration, on variable user-defined timescales, that can be examined prior to application to the data with applycal. The solutions for different fields and/or spectral windows can be interpolated in different ways, with all solutions stored in the same table.
The analog of accum in classic AIPS is the use of CLCAL to combine a series of (incremental) SN calibration tables to form successive (cumulative) CL calibration tables. AIPS SN/CL tables are the analog of ’G’ tables in CASA.
The only difference between incremental and cumulative calibration tables is that incremental tables are generated directly from the calibration solving tasks (gaincal, bandpass, etc.), and cumulative tables are generated from other cumulative and incremental tables via accum. In all other respects (internal format, application to data with applycal, plotting with plotcal, etc.), they are the same, and therefore interchangeable. Thus, accumulate and cumulative calibration tables need only be used when circumstances require it.
The accum task represents a generalization on the classic AIPS CLCAL (see sidebox) model of cumulative calibration in that its application is not limited to accumulation of ’G’ solutions. In principle, any basic calibration type can be accumulated (onto itself), as long as the result of the accumulation (matrix product) is of the same type. This is true of all the basic types, except ’D’. Accumulation is currently supported for ’B’, ’G’, and ’T’, and, in future, ’F’ (ionospheric Faraday rotation), delay-rate, and perhaps others. Accumulation of certain specialized types (e.g., ’GSPLINE’, ’TOPAC’, etc.) onto the basic types will be supported in the near future. The treatment of various calibration from ancillary data (e.g., system temperatures, weather data, WVR, etc.), as they become available, will also make use of accumulate to achieve the net calibration.
Note that accumulation only makes sense if treatment of a uniquely incremental solution is required (as described above), or if a careful interpolation or sampling of a solution is desired. In all other cases, re-solving for the type in question will suffice to form the net calibration of that type. For example, the product of an existing ’G’ solution and an amplitude and phase ’G’ self-cal (solved with the existing solution applied), is equivalent to full amplitude and phase ’G’ self-cal (with no prior solution applied), as long as the timescale of this solution is at least as short as that of the existing solution.
One obvious application is to calibrate the amplitudes and phases on different timescales during self-calibration. Here is an example:
# Add clean model
ft(vis='jupiter6cm.usecase.split.ms',
model='jupiter6cm.usecase.clean1.model')
# Phase only self-cal on 10s timescales
gaincal(vis='jupiter6cm.usecase.split.ms',
caltable='jupiter6cm.usecase.phasecal1',
gaintype='G',
calmode='p',
refant='6',
solint=10.0,
minsnr=1.0)
# Plot up solution phase and SNR
plotcal('jupiter6cm.usecase.phasecal1','','phase',antenna='1',subplot=211)
plotcal('jupiter6cm.usecase.phasecal1','','snr',antenna='1',subplot=212)
# Amplitude and phase self-cal on scans
gaincal(vis='jupiter6cm.usecase.split.ms',
caltable='jupiter6cm.usecase.scancal1',
gaintable='jupiter6cm.usecase.phasecal1',
gaintype='G',
calmode='ap',
refant='6',
solint='inf',
minsnr=1.0)
# Plot up solution amp and SNR
plotcal('jupiter6cm.usecase.scancal1','','amp',antenna='1',subplot=211)
plotcal('jupiter6cm.usecase.scancal1','','snr',antenna='1',subplot=212)
# Now accumulate these - they will be on the 10s grid
accum(vis='jupiter6cm.usecase.split.ms',
tablein='jupiter6cm.usecase.phasecal1',
incrtable='jupiter6cm.usecase.scancal1',
caltable='jupiter6cm.usecase.selfcal1',
interp='linear')
# Plot this up
plotcal('jupiter6cm.usecase.selfcal1','','amp',antenna='1',subplot=211)
plotcal('jupiter6cm.usecase.selfcal1','','phase',antenna='1',subplot=212)
The final plot is shown in Figure 4.9
ALERT: Only interpolation is offered in accum, no smoothing (as in smoothcal).
After the calibration solutions are computed and written to one or more calibration tables, one then needs to apply them to the data.
Alert: This section is written using the traditional applycal parameters. Users are encouraged to consult Appendix G for information on how to use the “Cal Library” to manage and apply calibration, which will ultimately provide more flexibility in calibration application, including in a wider variety of applications.
After all relevant calibration types have been determined, they must be applied to the target source(s) before splitting off to a new MS or before imaging. This is currently done by explicitly taking the data in the DATA column in the MAIN table of the MS, applying the relevant calibration tables, and creating the CORRECTED_DATA scratch column. The original DATA column is untouched.
The applycal task does this. The inputs are:
# applycal :: Apply calibrations solutions(s) to data
vis = '' # Name of input visibility file
field = '' # Select field using field id(s) or field name(s)
spw = '' # Select spectral window/channels
intent = '' # Select observing intent
selectdata = True # Other data selection parameters
timerange = '' # Select data based on time range
uvrange = '' # Select data within uvrange (default units meters)
antenna = '' # Select data based on antenna/baseline
scan = '' # Scan number range
observation = '' # Select by observation ID(s)
msselect = '' # Optional complex data selection (ignore for now)
docallib = False # Use callib or traditional cal apply parameters
gaintable = [] # Gain calibration table(s) to apply on the fly
gainfield = [] # Select a subset of calibrators from gaintable(s)
interp = [] # Interp type in time[,freq],
# per gaintable. default=linear,linear
spwmap = [] # Spectral windows
# combinations to form for gaintable(s)
calwt = [True] # Calibrate data weights per gaintable.
parang = False # Apply parallactic angle correction
applymode = '' # Calibration mode:
# ""="calflag","trial","flagonly", or "calonly"
flagbackup = True # Automatically back up the
# state of flags before the run?
As in other tasks, setting selectdata=True will open up the other selection sub-parameters (see § 2.3). In addition, you can also select data based on the scan intents that were set during the observations (find them through listobs). Many of the other parameters are the common calibration parameters that are described in § 4.4.1.
The single non-standard parameter is the calwt option to toggle the ability to scale the visibility weights by the inverse of the products of the scale factors applied to the amplitude of the antenna gains (for the pair of antennas of a given visibility). This should in almost all cases be set to its default (True). The weights should reflect the inverse noise variance of the visibility, and errors in amplitude are usually also in the weights.
Alert: Current (as of February 2014) Jansky VLA data has no calibrated weights to the data (unless they are created from switched power). To avoid trouble, calwt=False should be set for those data sets. Older, pre-Jansky VLA data should still be calibrated with calwt=True.
For applycal, the list of final cumulative tables is given in gaintable. In this case you will have run accum if you have done incremental calibration for any of the types, such as ’G’. You can also feed gaintable the full sets and rely on use of gainfield, interp and spwmap to do the correct interpolation and transfer. In particular, for frequency interpolation, the interpolation methods ending in ’PD’, nearestPD and linearPD also scale the phase by the frequency ratio between the measured and interpolated values. It is often more convenient to go through accumulation of each type with accum as described above (see § 4.5.6.2), as this makes it easier to keep track of the sequence of incremental calibration as it is solved and applied. You can also do any required smoothing of tables using smoothcal (§ 4.5.5), as this is not yet available in accum or applycal.
applycal has different applymodes: ’calflag’ will apply all flags from a calibration table to the data and apply the calibration itself to the remaining visibilities. ’trial’ will only report on the calibration table flags but not manipulate the data, ’flagonly’ applies the flags but not the calibration itself, and ’calonly’ will apply the calibration and but not the solution table flags. Data that would ’calflag’ would flag are thus passed through uncalibrated. This option can be useful when applycal is executed in consecutive steps, one calibration table at a time. Portions of the data that were not calibrated in the first run can then be calibrated in a second run with a different calibration table. This option should be used with care such that no uncalibrated data remains in the final data product.
applycal will flag all data that have no calibration solution. Flags will distribute into all of your scratch columns, i.e. it will affect your uncalibrated visibilities, too. To be able to restore the flags to the state before applycal is starting its duty, the task will make a backup of your current flags by default ( flagbackup=True). Restore them with flagmanager, if you are not happy with the applycal results.
If you are not doing polarization calibration or imaging, then you can set parang=False to make the calculations faster. If you are applying polarization calibration, or wish to make polarization images, then set parang=True so that the parallactic angle rotation is applied to the appropriate correlations. Currently, you must do this in applycal as this cannot be done on-the-fly in clean or mosaic. See § 4.4.1.3 for more on parang.
For example, to apply the final bandpass and flux-scaled gain calibration tables solutions to the NGC5921 data:
default('applycal')
vis='ngc5921.usecase.ms'
# We want to correct the calibrators using themselves
# and transfer from 1445+099 to itself and the target N5921
# Start with the fluxscale/gain and bandpass tables
gaintable=['ngc5921.usecase.fluxscale','ngc5921.usecase.bcal']
# pick the 1445+099 (field 1) out of the gain table for transfer
# use all of the bandpass table
gainfield = ['1','*']
# interpolation using linear for gain, nearest for bandpass
interp = ['linear','nearest']
# only one spw, do not need mapping
spwmap = []
# all channels, no other selection
spw = ''
selectdata = False
# no prior calibration
gaincurve = False
opacity = 0.0
# select the fields for 1445+099 and N5921 (fields 1 and 2)
field = '1,2'
applycal()
# Now for completeness apply 1331+305 (field 0) to itself
field = '0'
gainfield = ['0','*']
applycal()
# The CORRECTED_DATA column now contains the calibrated visibilities
In another example, we apply the final cumulative self-calibration of the Jupiter continuum data obtained in the example of § 4.5.6.2:
applycal(vis='jupiter6cm.usecase.split.ms',
gaintable='jupiter6cm.usecase.selfcal1',
selectdata=False)
Again, it is important to remember the relative nature of each calibration term. A term solved for in the presence of others is, in effect, residual to the others, and so must be used in combination with them (or new versions of them) in subsequent processing. At the same time, it is important to avoid isolating the same calibration effects in more than one term, e.g., by solving for both ’G’ and ’T’ separately (without applying the other), and then using them together.
It is always a good idea to examine the corrected data after calibration (using plotxy to compare the raw (’data’) and corrected (’corrected’) visibilities), as we describe next.
Once the source data is calibrated using applycal, you should examine the uv data and flag anything that looks bad. If you find source data that has not been flanked by calibration scans, delete it (it will not be calibrated).
For example, to look at the calibrated Jupiter data in the last example given in the previous section:
plotxy('jupiter6cm.usecase.split.ms','uvdist','amp','corrected',
selectdata=True,correlation='RR LL',fontsize = 14.0)
will show the CORRECTED_DATA column. See Figure 4.10.
Figure 4.10: The final ’amp’ versus ’uvdist’ plot of the self-calibrated Jupiter data, as shown in plotxy. The ’RR LL’ correlations are selected. No outliers that need flagging are seen.
See § 3.3 for a description of how to display and edit data using plotms or plotxy, and § 7.5 for use of the viewer to visualize and edit a Measurement Set.
Whenever calibration tasks are run, the models associated with the MS will be overwritten. Sometimes, however, one would like to completely remove the model and the task delmod can perform this functionality:
# delmod :: Deletes model representations in the MS
vis = '' # Name of input visibility file (MS)
otf = True # Delete the on-the-fly model data keywords
field = '' # Select field using field id(s) or field name(s)
scr = False # Delete the MODEL_DATA scr col (if it exists)
To do so, the parameter otf should be set to True. delmod can also be used if for any reason a MODEL_column was created and should be removed to avoid confusion between the on-the-fly model and the MODEL column (the MODEL_DATA column was required in CASA 3.3 and earlier). This can be achieved with scr=T.
delmod generally replaces the functionality of the older clearcal task. If one still decides to use the MODEL_DATA columns, however, clearcal is still useful and will reset both the MODEL_DATA and CORRECTED_DATA columns to DATA:
CASA <11>: inp clearcal
# clearcal :: Re-initializes the calibration for a visibility data set
vis = '' # Name of input visibility file (MS)
field = '' # Select field using field id(s) or field name(s)
spw = '' # Select spectral window/channel.
intent = '' # Select observing intent
addmodel = False # Add MODEL_DATA scratch column
with field, spw, and intent being data selection parameters. addmodel can be used to opt in/out of formation of the MODEL_DATA column.
With the introduction of the on-the-fly calculation of the MODEL visibilities, and the fact that applycal overwrites any previously existing CORRECTED_DATA column, clearcal is not required anymore unless usescratch=True is chosen in calibration tasks, and it is also not recommended to use clearcal to create the scratch columns at the beginning of data calibration; all benefits from the on-the-fly model would be made obsolete.
Starting in CASA 4.6.0, split2 has been renamed to split and became the default in CASA. The old implementation of split is still available under the name oldsplit. The interface of both implementations is the same, but the new split uses the MSTransform framework underneath. See § 4.7.4 for detailed information on mstransform.
The split task will apply calibration and output a new MS containing a specified list of sources (usually a single source). The inputs are:
# split :: Create a visibility subset from an existing visibility set
vis = '' # Name of input Measurement set or Multi-MS
outputvis = '' # Name of output Measurement set or Multi-MS
keepmms = True # If the input is a Multi-MS the output will also be a Multi-MS
field = '' # Select field using ID(s) or name(s)
spw = '' # Select spectral window/channels
scan = '' # Select data by scan numbers
antenna = '' # Select data based on antenna/baseline
correlation = '' # Correlation: '' ==> all, correlation='XX,YY'
timerange = '' # Select data by time range
intent = '' # Select data by scan intent.
array = '' # Select (sub)array by array ID number(s)
uvrange = '' # Select data bt baseline length
observation = '' # Select data by observation ID(s).
feed = '' # Multi-feed numbers: Not yet implemented.
datacolumn = 'corrected' # Which data column(s) to process
keepflags = True # Keep *completely flagged rows* instead of dropping them
width = 1 # Number of channels to average to form one output channel
timebin = '0s' # Bin width for time averaging
Usually you will run split with datacolumn=’corrected’ as previous operations (e.g. applycal) will have placed the calibrated data in the CORRECTED_DATA column of the MS. This will produce a new MS with this corrected data in its DATA column. The modes available in datacolumn are:
corrected
data
model
data,model,corrected
float_data
lag_data
float_data,data
lag_data,data
all
float_data is used for single dish processing.
For example, to split out 46 channels (5-50) from spw 1 of our NGC5921 calibrated dataset:
split(vis='ngc5921.usecase.ms',
outputvis='ngc5921.split.ms',
field='2', # Output NGC5921 data (field 2)
spw='0:5~50', # Select 46 chans from spw 0
datacolumn='corrected') # Take the calibrated data column
Time and channel averaging are available using the timebin and width parameters.
The timebin parameter gives the averaging interval. It takes a quantity, e.g.
timebin = '30s'
combine = ''
and by default will combine scans and states during averaging. For more details on the time average implementation and how it handles the weights, see § 4.7.4.4
The width parameter defines the number of channels to average to form a given output channel. This can be specified globally for all spw, e.g.
width = 5
or specified per spw, e.g.
width = [2,3]
to average 2 channels of 1st spectral window selected and 3 in the second one. See further details in § 4.7.4.5
ALERT: When averaging channels split will produce negative channel widths (as reported by listobs) if frequency goes down with increasing channel number, whether or not the input channel widths are negative. The bandwidths and channel resolutions will still be positive."
Sometimes the u,v,w coordinates of a Measurement Set are not recorded correctly by the correlator. In those cases, it may be necessary to recalculate them based on the antenna positions. fixvis will perform this task.
# fixvis :: Recalculates (u, v, w) and/or changes Phase Center
vis = '' # Name of the input visibility set.
outputvis = '' # Name of the output visibility set. (Can be the same
# as vis.)
field = '' # Fields to operate on. = all.
refcode = '' # reference frame to convert UVW coordinates to
reuse = True # base UVW calculation on the old values?
phasecenter = '' # use this direction as phase center
A useful feature of fixvis is that it can also change the phase center of a Measurement Set. This can be done with absolute coordinates or using offsets. An example is:
fixvis(vis='Moon.ms',outpuvis='Moon-fixed.ms',field='Moon', phasedir='J2000 9h25m00s 05d12m00s')
that will recalculate the u,v,w coordinates relative to the new phase center for the field ’Moon’.
For strong spectral line sources (like RFI sources), the Gibbs phenomenon may cause ringing across the frequency channels of an observation. This is called the Gibbs phenomenon and a proven remedy is the Hanning smoothing algorithm. Hanning smoothing is a running mean across the spectral axis with a triangle as a smoothing kernel. The central channel is weighted by 0.5 and the two adjacent channels by 0.25 to preserve the flux. Hanning smoothing significantly reduces Gibbs ringing but there’s no gain without a penalty and here it is the loss of a factor of two in spectral resolution.
The new hanningsmooth task (based on mstransform) has slightly changed in CASA 4.6. The main difference with oldhanningsmooth, the previous hanningsmooth implementation, is that the new task no longer writes to the input MS, but it always creates an output MS. The new task can also handle a Multi-MS and process it in parallel (see more information in the parallelization chapter (§ 10).
In CASA, the hanningsmooth task will apply Hanning smoothing to a spectral line uv data Measurement Set. The inputs are:
# hanningsmooth :: Hanning smooth frequency channel data to remove Gibbs ringing
vis = '' # Name of input Measurement set or Multi-MS.
outputvis = '' # Name of output Measurement set or Multi-MS.
keepmms = True # If the input is a Multi-MS the output will also
be a Multi-MS.
field = '' # Select field using ID(s) or name(s).
spw = '' # Select spectral window/channels.
scan = '' # Select data by scan numbers.
antenna = '' # Select data based on antenna/baseline.
correlation = '' # Correlation: '' ==> all, correlation='XX,YY'.
timerange = '' # Select data by time range.
intent = '' # Select data by scan intent.
array = '' # Select (sub)array(s) by array ID number.
uvrange = '' # Select data by baseline length.
observation = '' # Select by observation ID(s).
feed = '' # Multi-feed numbers: Not yet implemented.
datacolumn = 'all' # Input data column(s) to process.
The datacolumn parameter determines which of the data columns is to be Hanning smoothed: ’all’, ’model’, ’corrected’, ’data’, ’float_data’ or ’lag_data’. ’all’ will use whichever of the visibility data columns that are present in the input MS. If ’corrected’ is specified, the task will smooth the input CORRECTED_DATA column and save the smoothed data in DATA of the output MS.
mstransform is a multipurpose task that provides all the functionality of split, partition, cvel, hanningsmooth and applycal with the possibility of applying each of these transformations separately or together in an in-memory pipeline, thus avoiding unnecessary I/O steps. The list of transformations that mstransform can apply is as follows:
Notice that the order in the above list is not arbitrary. When various transformations are applied on the data using mstransform the order in which the transformations are pipe one after the other is the one shown in the above list.
Besides mstransform in itself, there are a series of tasks that mimic the old interfaces and are based on the mstransform framework: split, hanningsmooth and cvel2. Notice that as of CASA 4.6 the mstransform based versions of split and hanningsmooth are the default ones, whereas cvel still is based on the old implementation by default, and the cvel2 interface points to the mstransform implementation. For backwards compatibility safety the old implementation of split and hanningsmooth are available under the names of oldsplit and oldhanningsmooth.
A complete functionality description is available at the following web site:
http://www.eso.org/~scastro/ALMA/casa/MST/MSTransformDocs/MSTransformDocs.html
mstransform / split are able to create a new MS with a specific data selection, for instance splitting a science target. The new MS contains only the selected data and also the subtables are re-generated to contain only the metadata matching the data selection. The details about pure split operation are described in section (§ 4.7.1)
vis = '' # Name of input Measurement set or Multi-MS.
outputvis = '' # Name of output Measurement Set or Multi-MS.
tileshape = [0] # List with 1 or 3 elements giving the
tile shape of the disk data columns.
field = '' # Select field using ID(s) or name(s).
spw = '' # Select spectral window/channels.
scan = '' # Select data by scan numbers.
antenna = '' # Select data based on antenna/baseline.
correlation = '' # Correlation: '' ==> all, correlation='XX,YY'.
timerange = '' # Select data by time range.
intent = '' # Select data by scan intent.
array = '' # Select (sub)array(s) by array ID number.
uvrange = '' # Select data by baseline length.
observation = '' # Select by observation ID(s).
feed = '' # Multi-feed numbers: Not yet implemented.
datacolumn = 'corrected' # Which data column(s) to process.
keepflags = True # Keep *completely flagged rows* or
drop them from the output.
usewtspectrum = False # Create a WEIGHT_SPECTRUM column in the output MS.
The new features related with data election and re-indexing contained in mstransform / split but not in oldsplit are the following:
MSTransform is the framework used by the partition task, and even though it can be used directly to produce an MMS by specifying createmms = True it is highly recommended to use directly the partition task as explained in chapter (§ 10). Nevertheless, for the sake of completeness this is the list of expandable parameters shown when createmms = True
createmms = True # Create a multi-MS output
from an input MS.
separationaxis = 'auto' # Axis to do parallelization across
(scan,spw,baseline,auto)
numsubms = 9 # The number of Sub-MSs to create
(auto or any number)
As of CASA 4.5 mstransform incorporates the possibility of applying on the-the-fly (OTF) calibration by specifying docallib = True, which in turns allows to specify the “Cal Library” filename (callib parameter) whose format is described in Appendix G. This transformation is the first one applied to the data, producing effectively a corrected data column on-the-fly, which can be further transformed.
docallib = True # Enable OTF calibration
callib = '' # Cal Library filename
mstransform / split are both able to perform a regular (weighted) time average like oldsplit (in mstransform by specifying timeaverage = True and in split by default). However, there are some differences listed below. Additionally, mstransform it is able to perform a baseline dependent timeaverage as described in the paper Effects of Baseline Dependent Time Averaging of UV Data by W.D. Cotton (OBIT No. 13, 2008).
timeaverage = True # Average data in time.
timebin = '0s' # Bin width for time averaging.
timespan = '' # Span the timebin across scan, state or both.
maxuvwdistance = 0.0 # Maximum separation of start-to-end baselines
that can be included in an average(meters)
Both mstransform / split are able to perform a regular (weighted) channel average (in mstransform by specifying chanaverage =T rue and in split by default). The differences w.r.t. the channel average algorithm of oldsplit listed in the list below.
chanaverage = True # Average data in channels.
chanbin = 1 # Width (bin) of input channels
to average to form an output channel.
Both mstransform / hanningsmooth are able to perform hannining smooth (in mstransform by specifying hanning = True and in hanningsmooth by default). The only difference w.r.t. oldhanningsmooth is that a new MS is created.
Both mstransform / cvel2 are able combine SPWs (in mstransform by specifying combinespws = True and in cvel2 by default). The algorithm is in general the same as the old cvel, however, there are two significant differences in the new framework:
Both mstransform / cvel2 are able to perform spectral regridding / reference frame transformation (in mstransform by specifying regridms = True and in cvel2 by default). However mstransform is able to perform spectral regridding / reference frame transformation on each selected SPW separately, that is w/o combining the selected SPWs. As of CASA 4.5 both algorithms are fully aligned including the latest developments to take into account ephemerides and radial velocity correction.
regridms = True # Regrid the MS to a new spw,
channel structure or frame.
mode = 'channel' # Regridding mode
(channel/velocity/frequency/channel_b).
nchan = -1 # Number of channels in the output spw (-1=all).
start = 0 # First channel to use in the output spw
(mode-dependant)
width = 1 # Number of input channels that are
used to create an output channel.
nspw = 1 # Number of output spws to create in output MS.
interpolation = 'linear' # Spectral interpolation method.
phasecenter = '' # Image phase center: position or field index.
restfreq = '' # Rest frequency to use for output.
outframe = '' # Output reference frame (''=keep input frame).
veltype = 'radio' # Velocity definition.
A completely new feature in mstransform is the ability to separate an input SPW into several ones, or to combine various input SPWs into a single one with a uniform grid (resolving overlaps/gaps) to then separate it in several output SPWs. This option is activated under the regridding section (therefore by specifying regridms = True), together with the nspw) parameter which when bigger than 1 implies that the input SPW / combination of input SPWs should be separated:
regridms = True # Regrid the MS to a new spw,
nspw = 1 # Number of output spws to create in output MS.
The uvsub task will subtract the Fourier transform of the associated model of the MS (added to the MS with the tasks ft or setjy) from that in the CORRECTED_DATA column in the input MS and store the result in that same CORRECTED_DATA column.
The reverse operation is achieved by specifying reverse = True: in that case uvsub will add the value of the Fourier transform of the associated model to that in the CORRECTED_DATA column in the input MS and store the result in that same CORRECTED_DATA column.
The inputs are:
# uvsub :: Subtract/add model from/to the corrected visibility data.
vis = '' # Name of input visibility file (MS)
reverse = False # reverse the operation (add rather than subtract)
For example:
uvsub('ngc5921.split.ms')
ALERT: Currently, uvsub operates on the CORRECTED_DATA column in the MS vis. Eventually we will provide the option to write out a new MS.
At this point, consider whether you are likely to need continuum subtraction. If there is significant continuum emission present in what is intended as a spectral line observation, continuum subtraction may be desirable. You can estimate and subtract continuum emission in the uv-plane prior to imaging or wait and subtract an estimate of it in the image-plane. Note that neither method is ideal, and the choice depends primarily upon the distribution and strength of the continuum emission. Subtraction in the uv-plane is desirable if continuum emission dominates the source, since deconvolution of the line emission will be more robust if it not subject to the deconvolution errors of the brighter continuum. There is also a performance benefit since the continuum is nearly the same in each channel of the observation, and it is desirable to avoid repeating its deconvolution in each channel. However, doing the continuum estimation in the uv-plane has the serious drawback that interpolating visibilities between channels is only a good approximation for emission from near the phase center. Thus, uvcontsub will do an increasingly poor job for emission distributed further from the phase center. If the continuum emission is relatively weak, it is usually adequate to subtract it in the image plane; this is described in the Image Analysis section of this cookbook. Here, we describe how to do continuum subtraction in the uv-plane.
The uv-plane continuum subtraction is performed by the uvcontsub task. First, determine which channels in your data cube do not have line emission, perhaps by forming a preliminary image as described in the next chapter. This image will also help you decide whether or not you need to come back and do uv-plane continuum subtraction at all.
The inputs to uvcontsub are:
# uvcontsub :: Continuum fitting and subtraction in the uv plane
vis = '' # Name of input MS. Output goes to vis + ".contsub"
# (will be overwritten if already exists)
field = '' # Select field(s) using id(s) or name(s)
fitspw = '' # Spectral window:channel selection for fitting the
# continuum
combine = '' # Data axes to combine for the continuum estimation
# (none, or spw and/or scan)
solint = 'int' # Continuum fit timescale (int recommended!)
fitorder = 0 # Polynomial order for the fits
spw = '' # Spectral window selection for output
want_cont = False # Create vis + ".cont" to hold the continuum estimate.
For each baseline, and over the timescale specified in solint, uvcontsub will provide a polynomial fit to the real and imaginary parts of the (continuum-only) channels specified in fitspw (using the standard spw selection syntax), and then subtract this model from all channels specified in spw, or from all channels in spectral windows of fitspw if spw=’’. By setting the subparameter excludechannels=True, the channel selection in fitspw will be inverted. In that case one can select the line channels themselves and/or corrupted channels that are not used in the continuum fit to the data. fitspw can also take frequency ranges, e.g.
fitspw='*:113.767~114.528GHz;114.744~115.447GHz'
where ’*’ indicates to go across all spws.
Typically, low orders for the polynomial work best, like 0th (a constant), or 1st order (a linear fit). Use higher orders with caution and check your results carefully.
Usually, one should set solint=’int’ which does no averaging and fits each integration. However, if the continuum emission comes from a small region around the phase center and fitorder = 0, then you can set solint larger (as long as it is shorter than the timescale for changes in the visibility function of the continuum). If your scans are short enough you can also use scan averaging with combine=’scan’ and solint=’inf’. Be warned, setting solint too large will introduce “time smearing” in the estimated continuum and thus not properly subtract emission not at the phase center. Increasing solint speeds up the calculation but it does not improve the overall result quality of uvcontsub - although the continuum estimates of each baseline may be noisy (just like each visibility in a continuum MS may be noisy), it is better to use the ensemble of individual fits than to average the ensemble before fitting. Note that plotms can do time and baseline averaging on the fly to help you examine noisy data.
So, the recommended procedure is as follows:
For example, we perform uv-plane continuum subtraction on our NGC5921 dataset:
# Want to use channels 4-6 and 50-59 for continuum
uvcontsub(vis='ngc5921.usecase.ms',
field='N5921',
spw='', # all spw (only 0 in this data)
fitspw='0:4~7;50~59' # channels 4-6 and 50-59
solint='int', # leave it at the default
fitorder=0) # mean only
# You will see it made a new MS:
# ngc5921.usecase.ms.contsub"
Although not strictly a calibration operation, spectral regridding of a MS is available to aid in calibration operations (e.g. continuum subtraction) and preparation for imaging. For this purpose, the cvel task has been developed.
The inputs are:
# cvel :: regrid an MS to a new spectral window / channel structure or frame
vis = '' # Name of input measurement set
outputvis = '' # Name of output measurement set
passall = False # Pass through (write to output MS) non-selected data with
# no change
field = '' # Select field using field id(s) or field name(s)
spw = '' # Select spectral window/channels
selectdata = True # Other data selection parameters
timerange = '' # Range of time to select from data
array = '' # (sub)array indices
antenna = '' # Select data based on antenna/baseline
scan = '' # scan number range
mode = 'channel' # Regridding mode
nchan = -1 # Number of channels in output spw (-1=all)
start = 0 # first input channel to use
width = 1 # Number of input channels to average
interpolation = 'linear' # Spectral interpolation method
phasecenter = '' # Image phase center: position or field index
restfreq = '' # rest frequency (see help)
outframe = '' # Output frame (not case-sensitive, ''=keep input frame)
veltype = 'radio' # velocity definition
hanning = False # If true, Hanning smooth data before regridding to remove
# Gibbs ringing.
The key parameters for the operation of cvel are the regridding mode, the output reference outframe, veltype, restfreq (which may be a list of rest frequencies to match the different spws) and the standard selection parameters (in particular spw and field).
The syntax for mode options (’channel’,’velocity’,’frequency’,’channel_b’) has been made compatible with the respective modes of clean (§ 5.2.5). The combination of selected spw and mode will determine the output channels and spw(s):
spw = '0,1'; mode = 'channel'
# will produce a single spw containing all channels in spw 0 and 1
spw='0:5~28^2'; mode = 'channel'
# will produce a single spw made with channels (5,7,9,...,25,27)
spw = '0'; mode = 'channel': nchan=3; start=5; width=4
# will produce an spw with 3 output channels
# new channel 1 contains data from channels (5+6+7+8)
# new channel 2 contains data from channels (9+10+11+12)
# new channel 3 contains data from channels (13+14+15+16)
spw = '0:0~63^3'; mode='channel'; nchan=21; start = 0; width = 1
# will produce an spw with 21 channels
# new channel 1 contains data from channel 0
# new channel 2 contains data from channel 2
# new channel 21 contains data from channel 61
spw = '0:0~40^2'; mode = 'channel'; nchan = 3; start = 5; width = 4
# will produce an spw with three output channels
# new channel 1 contains channels (5,7)
# new channel 2 contains channels (13,15)
# new channel 3 contains channels (21,23)
The simplest use of cvel is to shift a single spectral window into an output frame without regridding. This is done with mode=’channel’. For example:
cvel(vis='test_w3oh_nohann.ms',
outputvis ='test_w3oh_nohann_chanbary.ms',
mode='channel',nchan=-1,start=0,width=1,
interpolation='linear',
phasecenter='',
spw='',
restfreq='1665.4018MHz',
outframe='BARY')
does this for an observation of the OH line.
There is also a special mode=’channel_b’ that does not force a linear output frequency grid, e.g. for irregularly spaced/overlapping spectral windows), but is nominally faster. This is not equivalent to a clean output gridding mode, although clean will work on this spectral lattice.
Mode channel is intended to not interpolate between channels. It will perform binning if needed. For most scientific applications we recommend using the mode=’velocity’’ and mode=’frequency’ options, as it is easiest to determine what the resulting channelization will be. For example:
cvel(vis='test_w3oh_nohann.ms',
outputvis ='test_w3oh_nohann_cvellsrk.ms',
mode='velocity',nchan=45,start='-35.0km/s',width='-0.55km/s',
interpolation='linear',
phasecenter='',
spw='',
restfreq='1665.4018MHz',
outframe='LSRK')
cvel(vis='test_w3oh_nohann.ms',
outputvis ='test_w3oh_nohann_cvelbary.ms',
mode='velocity',nchan=45,start='-35.0km/s',width='-0.55km/s',
interpolation='linear',
phasecenter='',
spw='',
restfreq='1665.4018MHz',
outframe='BARY')
will transform a MS into the LSRK and BARYcenter frames respectively.
The sign of the width parameter determines whether the channels run along increasing or decreasing values of frequency or velocity (i.e. if the cube is reversed or not).
The intent of cvel regridding is to transform channel labels and the visibilities to a spectral reference frame which is appropriate for the science analysis, e.g. from TOPO to LSRK, e.g. to correct for Doppler shifts throughout the time of the observation. Naturally, this will change the shape of the spectral features to some extent. According to the Nyquist theorem you should oversample a spectrum with twice the numbers of channels to retain the shape. Based on some tests, however, we recommend to observe with at least 3-4 times the number of channels for each significant spectral feature (like 3-4 channels per linewidth). This will minimize regridding artifacts in cvel.
If cvel has already established the grid that is desired for the imaging, clean should be run with the default channel mode ( > width=1) or with exactly the same frequency/velocity parameters as was used in cvel. This will avoid additional regridding in clean. Hanning smoothing is optionally offered in cvel, but tests have shown that already the regridding process itself, if it involved a transformation from TOPO to a non-terrestrial reference frame, implies some smoothing (due to channel interpolation) such that Hanning smoothing may not be necessary.
The interpolation method fftshift calculates the transformed visibilities by applying a FFT, then a phase ramp, and then an inverse FFT. It will also perform pre-averaging, if necessary (this will increas the S/N). Note that if you want to use this interpolation method, your frequency grid needs to be equidistant, i.e. it only works in mode velocity with veltype=radio, in mode frequency, and in mode channel (in the latter only if the input grid is itself equidistant in frequency). Note also that, as opposed to all other interpolation methods, this method will apply a constant (frequency-independent) shift in frequency which is not fully correct in the case of large fractional bandwidth of the given spectral window.
It is often desirable to fit simple analytic source component models directly to visibility data. Such fitting has its origins in early interferometry, especially VLBI, where arrays consisted of only a few antennas and the calibration and deconvolution problems were poorly constrained. These methods overcame the calibration uncertainties by fitting the models to calibration-independent closure quantities and the deconvolution problem by drastically limiting the number of free parameters required to describe the visibilities. Today, even with larger and better calibrated arrays, it is still desirable to use visibility model fitting in order to extract geometric properties such as the positions and sizes of discrete components in radio sources. Fits for physically meaningful component shapes such as disks, rings, and optically thin spheres, though idealized, enable connecting source geometry directly to the physics of the emission regions.
Visibility model fitting is carried out by the uvmodelfit task. The inputs are:
# uvmodelfit :: Fit a single component source model to the uv data:
vis = '' # Name of input visibility file
field = '' # field name or index
spw = '' # spectral window
selectdata = False # Activate data selection details
niter = 5 # Number of fitting iterations to execute
comptype = 'P' # Component type (P=pt source,G=ell. gauss,D=ell. disk)
sourcepar = [1, 0, 0] # Starting guess (flux,xoff,yoff,bmajaxrat,bpa)
varypar = [] # Which parameters can vary in fit
outfile = '' # Optional output component list table
ALERT: This task currently only fits a single component.
The user specifies the number of non-linear solution iterations (niter), the component type (comptype), an initial guess for the component parameters (sourcepar), and optionally, a vector of Booleans selecting which component parameters should be allowed to vary (varypar), and a filename in which to store a CASA componentlist for use in other applications (file). Allowed comptypes are currently point ’P’ or Gaussian ’G’.
The function returns a vector containing the resulting parameter list. This vector can be edited at the command line, and specified as input (sourcepar) for another round of fitting.
The sourcepar parameter is currently the only way to specify the starting parameters for the fit. For points, there are three parameters: I (total flux density), and relative direction (RA, Dec) offsets (in arcsec) from the observation’s phase center. For Gaussians, there are three additional parameters: the Gaussian’s semi-major axis width (arcsec), the aspect ratio, and position angle (degrees). It should be understood that the quality of the result is very sensitive to the starting parameters provided by the user. If this first guess is not sufficiently close to the global χ2 minimum, the algorithm will happily converge to an incorrect local minimum. In fact, the χ2 surface, as a function of the component’s relative direction parameters, has a shape very much like the inverse of the absolute value of the dirty image of the field. Any peak in this image (positive or negative) corresponds to a local χ2 minimum that could conceivable capture the fit. It is the user’s responsibility to ensure that the correct minimum does the capturing.
Currently, uvmodelfit relies on the likelihood that the source is very near the phase center (within a beamwidth) and/or the user’s savvy in specifying the starting parameters. This fairly serious constraint will soon be relieved somewhat by enabling a rudimentary form of uv-plane weighting to increase the likelihood that the starting guess is on a slope in the correct χ2 valley.
Improvements in the works for visibility model fitting include:
Example (see Figure 4.11):
#
# Note: It's best to channel average the data if many channels
# before running a modelfit
#
split('ngc5921.ms','1445+099_avg.ms',
datacolumn='corrected',field='1445*',width='63')
# Initial guess is that it's close to the phase center
# and has a flux of 2.0 (a priori we know it's 2.47)
uvmodelfit('1445+099_avg.ms', # use averaged data
niter=5, # Do 5 iterations
comptype='P', # P=Point source, G=Gaussian, D=Disk
sourcepar=[2.0,.1,.1], # Source parameters for a point source
spw='0', #
outfile='gcal.cl') # Output component list file
# Output looks like:
There are 19656 - 3 = 19653 degrees of freedom.
iter=0: reduced chi2=0.0418509: I=2, dir=[0.1, 0.1] arcsec
iter=1: reduced chi2=0.003382: I=2.48562, dir=[-0.020069, -0.0268826] arcsec
iter=2: reduced chi2=0.00338012: I=2.48614, dir=[0.00323428, -0.00232235] arcsec
iter=3: reduced chi2=0.00338012: I=2.48614, dir=[0.00325324, -0.00228963] arcsec
iter=4: reduced chi2=0.00338012: I=2.48614, dir=[0.00325324, -0.00228963] arcsec
iter=5: reduced chi2=0.00338012: I=2.48614, dir=[0.00325324, -0.00228963] arcsec
If data weights are arbitrarily scaled, the following formal errors
will be underestimated by at least a factor sqrt(reduced chi2). If
the fit is systematically poor, the errors are much worse.
I = 2.48614 +/- 0.0176859
x = 0.00325324 +/- 0.163019 arcsec
y = -0.00228963 +/- 0.174458 arcsec
Writing componentlist to file: /home/sandrock/smyers/Testing/Patch2/N5921/gcal.cl
# Fourier transform the component list to a model of the MS
ft('1445+099_avg.ms', complist='gcal.cl')
# Plot data versus uv distance
plotxy('1445+099_avg.ms', xaxis='uvdist', datacolumn='corrected')
# Specify green circles for model data (overplotted)
plotxy('1445+099_avg.ms', xaxis='uvdist', datacolumn='model',
overplot=True, plotsymbol='go')
Figure 4.11: Use of plotxy to display corrected data (red and blue points) and uv model fit data (green circles).
Alert: statwt is still an experimental task. Please check the results carefully and report any problems to the NRAO CASA helpdesk.
In most cases, the data that comes from the telescopes have the correct absolute or relative weights associated (absolute weights will be supplied once the VLA switched power application becomes standard; for ALMA the Tsys application is already in place). However, there are data sets where one would like to adjust the weights based on the scatter of the visibilities (typically as a function of time, antenna, and/or baseline). This calculation is performed by the task statwt that updates the WEIGHT and SIGMA columns of the Measurement Set. statwt inputs are:
# statwt :: Reweight visibilities according to their scatter
vis = '' # Name of measurement set
dorms = False # Use rms instead of stddev?
byantenna = False # Estimate the noise per antenna -not
# implemented (vs. per baseline)
fitspw = '' # The signal-free spectral window:channels
# to estimate the scatter from
fitcorr = '' # The signal-free correlation(s) to estimate
# the scatter from (not implemented)
combine = '' # Let estimates span changes in spw, corr,
# scan and/or state
timebin = '0s' # Bin length for estimates (not implemented)
minsamp = 2 # Minimum number of unflagged visibilities
# for estimating the scatter
field = '' # Select field using ID(s) or name(s)
spw = '' # Select spectral window/channels
antenna = '' # Select data based on antenna/baseline
timerange = '' # Select data by time range
scan = '' # Select data by scan numbers
intent = '' # Select data by scan intents
array = '' # Select (sub)array(s) by array ID number
correlation = '' # Select correlations to reweight (DEPRECATED in CASA v4.5)
observation = '' # Select by observation ID(s)
datacolumn = 'corrected' # Which data column to calculate the scatter
# from
statwt should only be run after all calibration steps have been performed. The parameter dorms switches from a scatter standard deviation to a root mean square scatter estimator. datacolumn specifies the column on which the task operates and the usual data selection parameters apply. Channels with strong RFI or a spectral line should be avoided for the calculation and good channel range should be specified via fitspw. In its current implementation, statwt uses data samples of an integration time interval but eventually wider sample intervals can be specified by the timebin parameter. Those samples are contained within a scan, spw, and polarization product but using the combine can relax this restriction. minsamp sets the minimum number of unflagged visibilities used for the calculation.
Alert: As of CASA 4.5, selection using correlation has been deprecated in statwt; in prior versions, this was not working correctly, and it is unlikely setting weights in a correlation-dependent manner is advisable.
conjugatevis is an easy task to flip the signs of the phases of visibilities, thus creating the complex conjugate numbers. The inputs are like:
# conjugatevis :: Change the sign of the phases in all visibility columns.
vis = '' # Name of input visibility file.
spwlist = '' # Spectral window selection
outputvis = '' # Name of output visibility file
overwrite = False # Overwrite the outputvis if it exists.
The task works on all scratch columns.
(This task was added early in JVLA commissioning, when some data suffered from a phase sign error.)
When an astronomical object has a proper motion, in particular objects in our solar system, a static (RA,Dec) position in the FIELD table of the MeasurementSet will not accurately describe the time-dependent position. Prior to CASA 4.2, there was no support for ephemeris objects other than the built-in reference frames for the Sun, the Moon, and the planets out to PLUTO. With CASA 4.2, several new features were introduced which help the user to attach an arbitrary ephemeris to a given field and work with the object from calibration to imaging.
The CASA format for ephemeris tables has been defined was introduced in the early development stages of CASA in connection with the Measures module. The me tool permits position calculations based on ephemerides in this format. Two examples for such tables can be found in the distribution directory in subdirectory data/ephemerides: VGEO is an ephemeris of Venus in the geocentric reference frame while VTOP is an ephemeris for the same object in the TOPO reference fame for the observatory location of the VLA. With the introduction of solar system source models (Butler) in the setjy task, a nearly complete set of ephemerides for the larger bodies in our solar system had to be made available. These are stored in nearly the same format as the above examples VGEO and VTOP (but with a few enhancements) in directory data/ephemerides/JPL-Horizons. If your object’s ephemeris is among those stored in data/ephemerides/JPL-Horizons, you can simply copy the ephemeris from there. Otherwise, you can request the ephemeris from the JPL-Horizons using the CASA commands (for example)
import recipes.ephemerides.request as jplreq
jplreq.request_from_JPL(objnam='Mars',startdate='2012-01-01',enddate='2013-12-31',
date_incr='0.1 d', get_axis_orientation=False,
get_axis_ang_orientation=True,
get_sub_long=True, use_apparent=False, get_sep=False,
return_address='YOUR_EMAIL_ADDESS',
mailserver='YOUR_MAIL_SERVER_ADDRESS')
where you need to fill in the parameters objnam, startdate, enddate,date_incr (the time interval between individual ephemeris table entries), return_address (your email address where you want to receive the ephemeris), and mailserver (the smtp server through which you want to send the request email). The other parameters should be set as shown. Within a short time, you should receive the requested ephemeris as an email from NASA’s JPL-Horizonssystem. Save the email into a file with the “save as” function of your mail client. See the next section on how to attach it to your dataset.
As of CASA 4.6.0, importasdm fills the SOURCE coodinates with the correct postions based on the ephemerides table. Alternatively, one can use the task fixplanets to set the ephemeris of a given field in a Measurement Set. Her is an example:
fixplanets(vis='uid___A002_X1c6e54_X223.ms',
field='Titan', fixuvw=True, direction='mytitanephemeris')
where you need to set the parameters vis to the name of your MS and the parameter field to the name of the field to which you want to attach the ephemeris. The parameter direction must be set to the name of your ephemeris table. Accepted formats are (a) the CASA format (as in VGEO or the ephemerides in data/ephemerides/JPL-Horizons as described above) and (b) the JPL-Horizons mail format which you obtain by saving an ephemeris email you received from JPL-Horizons. The parameter fixuvw should be set to True in order to trigger a recalculation of the UVW coordinates in your MS based on the new ephemeris. The task fixplanets can also be used for other field direction modifications. Please refer to the help text of the task.
Note that among the ephemerides in the directory data /ephemerides/JPL-Horizons/ you should only use those ending in ’_J2000.tab’. They are the ones in J2000 coordinates.
Once you have attached the ephemeris to a field of an MS, it will automatically be handled in tasks like split and concat which need to hand on the ephemeris to their output MSs. In particular concat recognizes when fields of the MSs to be concatenated use the same ephemeris and merges these fields if the time coverage of the provided ephemeris in the first MS also covers the observation time of the second MS. The ephemeris of the field in the first MS will then be used for the merged field. In order to inspect the ephemeris attached to a field, you can find it inside the FIELD subdirectory of your MS. The optional column EPHEMERIS_ID in the FIELD table points to the running number of the ephemeris table. A value of −1 indicates that no ephemeris is attached. Note that in case an ephemeris is attached to a field, the direction column entries for that field in the FIELD table will be interpreted as an offset to the ephemeris direction and are therefore set to (0.,0.) by default. This offset feature is used in mosaic observations where several fields share the same ephemeris with different offsets. The TIME column in the FIELD table should be set to the beginning of the observation for that field and serves as the nominal time for ephemeris queries.
The ephemerides contain radial velocity information. The task cvel can be used to transform spectral windows into the rest frame of the ephemeris object by setting the parameter outframe to “SOURCE” as in the following example:
cvel(vis='europa.ms',
outputvis='cvel_europa.ms', outframe='SOURCE', mode = 'velocity',
width = '0.3km/s', restfreq = '354.50547GHz')
This will make cvel perform a transformation to the GEO reference frame followed by an additional Doppler correction for the radial velocity given by the ephemeris for the each field. (Typically, this should happen after calibration and after splitting out the spectral widows and the target of interest). The result is an MS with a single combined spectral window in reference frame REST. From this frame, further transformations to other reference frames are not possible.
The ALMA Science Data Model (the raw data format for ALMA data) now foresees an Ephemeris table. This feature has been in use since the beginning of ALMA Cycle 3 both for science targets and calibrator objects. With ALMA Cycle 3 (or later) datasets, the task importasdm will automatically translate the contents of the ASDM Ephemeris table into separate ephemeris tables in CASA format and attach them to the respective fields.
In the case of mosaics, all fields belonging to a mosaic on an ephemeris object will share the same ephemeris. The individual mosaic pointings will use the offset mechanism described above to define the position of each field.
The data reduction tutorials on casaguides.nrao.edu provide walkthroughs for high and low frequency, spectral line and polarization calibration techniques.