Antennae Band7 - Calibration for CASA 3.3
- This script assumes that you have downloaded Antennae_Band7_UnCalibratedMSAndTablesForReduction.tgz from AntennaeBand7#Obtaining_the_Data
- Details of the ALMA observations are provided at AntennaeBand7
- This portion of the guide covers calibration of the raw visibility data. To skip to the imaging portion of the guide, see: Antennae Band7 - Imaging.
Unpack the Data
Once you have downloaded the Antennae_Band7_UnCalibratedMSandTablesForReduction.tgz, unpack the file in a terminal outside CASA using
>tar -xvzf Antennae_Band7_UnCalibratedMSandTablesForReduction.tgz
then change directory (cd) to the directory Antennae_Band7_UnCalibratedMSandTablesForReduction.
You may wish to type
>ls
to look at the files present. You should see a bunch of files with extension ".ms" indicating that these are CASA measurement set (MS) files. The data have already been converted to MS format using the CASA task importasdm. Accompanying the data are some basic calibration tables holding system temperature (Tsys) and water vapor radiometer (WVR) information that we have generated outside of CASA (for Early Science CASA will be able to generate these).
To begin, start CASA by typing
>casapy
Initial Inspection
First we will take stock of what we have. If you have not already done so, begin by reviewing the description of the observations here AntennaeBand7 . The 10 data sets each target either the northern or the southern mosaic, as follows:
Northern Mosaic:
- uid___A002_X1ff7b0_Xb.ms
- uid___A002_X207fe4_X3a.ms
- uid___A002_X207fe4_X3b9.ms
- uid___A002_X2181fb_X49.ms
Southern Mosaic:
- uid___A002_X1ff7b0_X1c8.ms
- uid___A002_X207fe4_X1f7.ms
- uid___A002_X207fe4_X4d7.ms
- uid___A002_X215db8_X18.ms
- uid___A002_X215db8_X1d5.ms
- uid___A002_X215db8_X392.ms
The first step is to get basic information about the data: targets observed, time range, spectral setup, and so on. We do this using the task listobs, which will output a detailed summary of each dataset.
# In CASA
# Define a python list holding the names of all of our data sets
basename_all=["uid___A002_X1ff7b0_Xb","uid___A002_X207fe4_X3a","uid___A002_X207fe4_X3b9",
"uid___A002_X2181fb_X49","uid___A002_X1ff7b0_X1c8","uid___A002_X207fe4_X1f7",
"uid___A002_X207fe4_X4d7","uid___A002_X215db8_X18","uid___A002_X215db8_X1d5",
"uid___A002_X215db8_X392"]
# Loop over each element in the list and
for asdm in basename_all:
listobs(vis=asdm+'.ms', listfile=asdm+'.listobs.txt', verbose=True)
These commands define a python list called "basename_all", which contains the name of all 10 MS files. The for loop executes for each element in basename_all, calling listobs and directing the output a file called, e.g., "uid___A002_X1ff7b0_Xb.listobs.txt" for the first measurement set.
You can browse through the listobs output as you would normally look at a text file (use emacs, vi, or another editor). You can also send the output to the terminal from inside of CASA. To do so type:
# In CASA
cat uid___A002_X1ff7b0_Xb.listobs.txt
or
# In CASA
os.system('more uid___A002_X1ff7b0_Xb.listobs.txt')
CASA knows a few basic shell commands like 'cat', 'ls', and 'rm' but for more complex commands you will need to run them inside 'os.system("command")'. For more information see http://casa.nrao.edu/ .
Here is an example of the (abridged) output from listobs for the first dataset in the list, uid___A002_X1ff7b0_Xb.ms, which targets the Northern Mosaic:
============================================================================= MeasurementSet Name:/Users/despada/Desktop/Imaging/Antennae/Datasets/band7/uid___A002_X1ff7b0_Xb.ms ============================================================================= Observer: Unknown Project: T.B.D. Observation: ALMA(11 antennas) Data records: 181357 Total integration time = 4931.71 seconds Observed from 28-May-2011/01:25:27.6 to 28-May-2011/02:47:39.3 (UTC) Fields: 26 ID Code Name RA Decl Epoch SrcId 0 none 3c279 12:56:11.1666 -05.47.21.5247 J2000 0 1 none Titan 12:42:43.9481 -01.43.38.3190 J2000 1 2 none NGC4038 - A* 12:01:53.1701 -18.52.37.9200 J2000 2 3 none NGC4038 - A* 12:01:51.9030 -18.51.49.9437 J2000 2 4 none NGC4038 - A* 12:01:52.4309 -18.51.49.9437 J2000 2 5 none NGC4038 - A* 12:01:52.9587 -18.51.49.9437 J2000 2 6 none NGC4038 - A* 12:01:53.4866 -18.51.49.9436 J2000 2 7 none NGC4038 - A* 12:01:54.0144 -18.51.49.9436 J2000 2 8 none NGC4038 - A* 12:01:52.1669 -18.51.56.4319 J2000 2 9 none NGC4038 - A* 12:01:52.6948 -18.51.56.4318 J2000 2 10 none NGC4038 - A* 12:01:53.2226 -18.51.56.4318 J2000 2 11 none NGC4038 - A* 12:01:53.7505 -18.51.56.4318 J2000 2 12 none NGC4038 - A* 12:01:51.9030 -18.52.02.9201 J2000 2 13 none NGC4038 - A* 12:01:52.4309 -18.52.02.9200 J2000 2 14 none NGC4038 - A* 12:01:52.9587 -18.52.02.9200 J2000 2 15 none NGC4038 - A* 12:01:53.4866 -18.52.02.9200 J2000 2 16 none NGC4038 - A* 12:01:54.0144 -18.52.02.9199 J2000 2 17 none NGC4038 - A* 12:01:52.1669 -18.52.09.4082 J2000 2 18 none NGC4038 - A* 12:01:52.6948 -18.52.09.4082 J2000 2 19 none NGC4038 - A* 12:01:53.2226 -18.52.09.4082 J2000 2 20 none NGC4038 - A* 12:01:53.7505 -18.52.09.4081 J2000 2 21 none NGC4038 - A* 12:01:51.9030 -18.52.15.8964 J2000 2 22 none NGC4038 - A* 12:01:52.4309 -18.52.15.8964 J2000 2 23 none NGC4038 - A* 12:01:52.9587 -18.52.15.8963 J2000 2 24 none NGC4038 - A* 12:01:53.4866 -18.52.15.8963 J2000 2 25 none NGC4038 - A* 12:01:54.0144 -18.52.15.8963 J2000 2 (nVis = Total number of time/baseline visibilities per field) Spectral Windows: (9 unique spectral windows and 2 unique polarization setups) SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 0 4 TOPO 184550 1500000 7500000 183300 I 1 3840 TOPO 344845.586 488.28125 1875000 344908.33 XX YY 2 1 TOPO 343908.086 1875000 1875000 344908.33 XX YY 3 3840 TOPO 356845.586 488.28125 1875000 344908.33 XX YY 4 1 TOPO 343908.086 1875000 1875000 344908.33 XX YY 5 128 TOPO 344900.518 15625 2000000 344908.33 XX YY 6 1 TOPO 343892.705 1796875 1796875 344908.33 XX YY 7 128 TOPO 356900.518 15625 2000000 344908.33 XX YY 8 1 TOPO 343892.705 1796875 1796875 344908.33 XX YY Antennas: 11 'name'='station' ID= 0-3: 'DV02'='A015', 'DV04'='J505', 'DV06'='T704', 'DV07'='A004', ID= 4-7: 'DV08'='A072', 'DV09'='A008', 'DV10'='A009', 'DV11'='A016', ID= 8-10: 'PM01'='T702', 'PM02'='A017', 'PM03'='J504' ================================================================================
And here is an example of the listobs for uid___A002_X1ff7b0_X1c8.ms, which targets the Southern Mosaic:
Insert listobs from a Southern mosaic data set here.
This output shows that three sources were observed in each data set: 3c279, Titan, and the Antennae.
- The Antennae are our science target. Note that the source name changes between the Northern and Southern Mosaic and that the source corresponds to a number of individual fields (the field id column). These are the individual mosaic pointings. There are 23 for the Northern Mosaic and 29 for the Southern Mosaic.
- Titan will be used to set the absolute flux scale of the data.
- 3c279 plays two roles: it will serve as our bandpass calibrator, used to characterize the frequency response of the telescopes, and as our phase calibrator, which we will use to track changes in the phase and amplitude response of the telescopes over time. Observations of 3c279 are interleaved with observations of the Antennae.
The output also shows that the data contain many spectral windows, using the labeling scheme in the listobs above. These are:
- spw 0 targets ~185 GHz and holds water vapor radiometer data
- spw 1 and spw 3 hold our science data. These are "Frequency Domain Mode" (FDM) data with small (0.49 kHz) channel width and wide (1.875 GHz) total bandwidth. As a result these have a lot of channels (3840). spw 1 holds the lower sideband (LSB) data and includes the CO(3-2) line. We will focus on these data. For the CO(3-2) line the channel width corresponds to 0.426 km/s and the bandwidth of spw 1 to 1634 km/s.
- spw 2 and spw 4 hold frequency-averaged versions of spw 1 and 3 ("Channel 0" for those familiar with AIPS). These are used for quick inspection. We will not use them here.
- spw 5 and spw 7 hold lower resolution processing ("Time Domain Mode", TDM) data for the same part of the spectrum (baseband) as spws 1 and 3. These data have only 128 channels across 2 GHz bandwidth and so have a much coarser channel spacing than the FDM data. These were used to generate the calibration tables that we include in the tarball but will not otherwise appear in this guide.
The final column of the listobs output in the logger (not shown above) gives the scan intent. Later we will use this information to flag the pointing scans and the hot and ambient load calibration scans.
We'll now have a look at the configuration of the antennas used to take the data using the task plotants (Figure 1).
# In CASA
basename_all=["uid___A002_X1ff7b0_Xb","uid___A002_X207fe4_X3a","uid___A002_X207fe4_X3b9",
"uid___A002_X2181fb_X49","uid___A002_X1ff7b0_X1c8","uid___A002_X207fe4_X1f7",
"uid___A002_X207fe4_X4d7","uid___A002_X215db8_X18","uid___A002_X215db8_X1d5",
"uid___A002_X215db8_X392"]
for asdm in basename_all:
print "Antenna configuration for : "+asdm
plotants(vis=asdm+'.ms', figfile=asdm+'.plotants.png')
dummy_string = raw_input("Hit <Enter> to see the antenna configuration for the next data set.")
This will loop through all 10 data sets, show you the antenna position for each, and save that as a file named, e.g., "uid___A002_X1ff7b0_Xb.plotants.png" for the first data set. The "raw_input" command asks CASA to wait for your input before proceeding.
How to Deal With 10 Measurement Sets
It should already be clear from the initial inspection that dealing with 10 data sets at the same time can be a bit problematic. This is especially tricky in our case because the Antennae data contain two distinct sets of observations: the Northern and Southern Mosaics. The source name changes between these two scripts and there are different numbers of fields in the mosaic.
As a general rule one would reduce each individual observation separately or at the very least only group data observed in a uniform way and very close in time.
Unfortunately, a CASA Guide stepping through the reduction for each of 10 data sets would quickly become unwieldy. Therefore we will use a few tricks to reduce the Antennae data in a kind of batch mode. You have already seen the first trick: we can define a python list holding the names of each data set and then loop over this list to execute the same command on each data set. For example:
# In CASA
basename_all=["uid___A002_X1ff7b0_Xb","uid___A002_X207fe4_X3a","uid___A002_X207fe4_X3b9",
"uid___A002_X2181fb_X49","uid___A002_X1ff7b0_X1c8","uid___A002_X207fe4_X1f7",
"uid___A002_X207fe4_X4d7","uid___A002_X215db8_X18","uid___A002_X215db8_X1d5",
"uid___A002_X215db8_X392"]
for asdm in basename_all:
print asdm
Note that after cutting and pasting a 'for' loop like this you often have to press return twice to execute. You may also want to take care to paste a line at a time if you are having trouble copy and pasting. Alternatively, you can try "cpaste" to paste blocks of code. To do so type "cpaste" at the CASA prompt, paste your commands, and then type "--" and hit return on the final (otherwise empty) line. This should look something like this:
CASA <8>: cpaste Pasting code; enter '--' alone on the line to stop. :basename_all=["uid___A002_X1ff7b0_Xb","uid___A002_X207fe4_X3a","uid___A002_X207fe4_X3b9", : "uid___A002_X2181fb_X49","uid___A002_X1ff7b0_X1c8","uid___A002_X207fe4_X1f7", : "uid___A002_X207fe4_X4d7","uid___A002_X215db8_X18","uid___A002_X215db8_X1d5", : "uid___A002_X215db8_X392"] : :for asdm in basename_all: : print asdm :-- uid___A002_X1ff7b0_Xb uid___A002_X207fe4_X3a uid___A002_X207fe4_X3b9 uid___A002_X2181fb_X49 uid___A002_X1ff7b0_X1c8 uid___A002_X207fe4_X1f7 uid___A002_X207fe4_X4d7 uid___A002_X215db8_X18 uid___A002_X215db8_X1d5 uid___A002_X215db8_X392 CASA <9>:
if you have trouble, just carefully paste one line at a time directly into CASA and hit return until the desired command executes.
You only need to define your list of MS files once per CASA session. Then "basename_all" will be a variable in the casapy shell. You can check if it exists by typing "print basename_all". In the interests of allowing you to easily exit and restart CASA and pick this guide up at any point we will redefine "basename_all" in each section of the guide. Feel free to skip this step if you've already defined it in your session.
This page will step you through the reduction of the whole Antennae Band 7 SV data set using these for loops. We will not be able to show every diagnostic plot but we give an example of each and the syntax to generate the rest.
A Priori Flagging
Even before we look in detail, we know that there are some data that we wish to exclude. We will start by flagging "shadowed" data where one antenna blocks the line of sight of another. We will also flag scans that were used to carry out pointing and atmospheric calibration, identified by their scan intent. Finally, we'll flag the autocorrelation data (the correlation of the signal from an antenna with itself) as are only interested in cross-correlation data to make an interferometric image.
Start by defining our list of MS files:
# In CASA
basename_all=["uid___A002_X1ff7b0_Xb","uid___A002_X207fe4_X3a","uid___A002_X207fe4_X3b9",
"uid___A002_X2181fb_X49","uid___A002_X1ff7b0_X1c8","uid___A002_X207fe4_X1f7",
"uid___A002_X207fe4_X4d7","uid___A002_X215db8_X18","uid___A002_X215db8_X1d5",
"uid___A002_X215db8_X392"]
Then flag shadowed data using the command flagdata:
# In CASA
for asdm in basename_all:
flagdata(vis=asdm+'.ms',mode = 'shadow', diameter=12.0, flagbackup = F,)
In the flagdata task we choose:
- vis = asdm+'.ms' : each measurement set
- mode='shadow',diameter=12.0: flag shadowed data, taking into account that antennas are 12m diameter
- flagbackup = F: Do not automatically back up the flag files. We will save all of the a priori flags together using flagmanager at the end of this subsection.
Now flag the pointing and atmospheric calibration scans using flagdata in 'manualflag' mode and selecting on 'intent':
# In CASA
for asdm in basename_all:
flagdata(vis=asdm+'.ms', mode='manualflag', intent='*POINTING*',flagbackup = F)
flagdata(vis=name+'.ms', mode='manualflag', intent='*ATMOSPHERE*',flagbackup = F)
Now flag the autocorrelation data with flagautocorr.
# In CASA
for asdm in basename_all:
flagautocorr(vis=asdm+'.ms')
Finally store the current flagging information using flagmanager:
# In CASA
for asdm in basename_all:
flagmanager(vis = asdm+'.ms', mode = 'save', versionname = 'Apriori')
We can now roll back the flags to match the current version, called 'Apriori', whenever we want. It would have been possible to set flagdata to flagbackup=T so that it stores the flags at each of the flagging step automatically, but this way is a bit more efficient.
Examine Tsys and WVR Calibration Tables, Apply, and Split
The Antennae_Band7_UnCalibratedMSandTablesForReduction directory includes system temperature (Tsys) and water vapor radiometer (WVR) calibration tables, which appear as files with extensions '.tsys.cal.fdm' and '.wvr.cal'. These tables have been built from the spw 0 (WVR) and spws 5 and 7 (Tsys) and provided to you because CASA does not generate them at the moment (this situation will change soon). The Tsys calibration corrects (to first-order) for the atmospheric opacity as a function of time and frequency and associates weights with each visibility that persist through imaging. The Tsys calibrations were derived from the TDM data and then interpolated to the FDM frequency coverage. The WVR calibration uses observations of atmospheric water lines to correct for phase variations as a function of time.
We inspect the Tsys tables for the spectral window spw=1 with the task plotcal. We want to check that Tsys data have reasonable values. Estimated Tsys under average atmospheric observations for a given band can be found for example in the sensitivity calculator.
We plot the Tsys for all the antennas and polarizations (XX and YY) as a function of frequency for the first dataset.
#In CASA
asdm=basename_all[0]
plotcal(caltable=asdm+'.tsys.cal.fdm',
xaxis="freq",yaxis="amp",
timerange='<1e8',spw='1',
plotsymbol=".")
Here spw='1' means the spectral window where CO(3-2) is. A large value for the time range (timerange='<1e8') is chosen to force averaging within each scan in the dataset. Not averaging all scans allow us to see any large deviation as a function of time as well.
In Figure 2, the Tsys solutions look reliable, with typical values ~150 K, except for some large values of Tsys at ~300 and 400 K (green and blue lines). By plotting each antenna individually, one finds that this corresponds to the two polarizations of DV04. Therefore, the data for that antenna will need to be flagged. There is also a mesospheric absorption line at about 345.2 GHz that makes Tsys larger, as seen in all of the other antennas. By applying the Tsys calibration tables in the data we will minimize the contribution of these absorption lines.
We now apply the Tsys and the WVR calibration tables to the data with applycal. Again we loop through all the datasets and sources to apply the Tsys and WVR calibration tables. It is important to only apply Tsys and WVR corrections obtained close in time to the data being corrected, so in addition to looping over data sets we define the list of unique source names and loop over these. Then by setting gainfield and field to the same value we ensure that Tsys and WVR calibrations are only applied to the source for which they area measured.
# In CASA
field_names = ['Titan','3c279','NGC4038*','Antenna*']
for asdm in basename_all:
for field in field_names:
applycal(vis=asdm+".ms", spw='1',
field=field, gainfield=field,
interp='nearest',
gaintable=[asdm+".tsys.cal.fdm",asdm+'.wvr.cal'],
flagbackup=F)
where:
- field: the field to which we will apply the calibration,
- gainfield: the field from which we wish to take the calibration table
- spw='1' : select only spectral window 1
- interp='nearest': use the interpolation mode to the 'nearest' solution.
Now you can use plotms to show some of the before-and-after effects of applying the calibration. Just run plotms to plot the amplitude versus channel in order to check if atmospheric absorption lines have been effectively corrected.
# In CASA
asdm=basename[0]
plotms(vis=asdm+'.ms',
field='3c279',
xaxis='channel', yaxis='amp',
selectdata=T, spw='1', correlation='XX',
avgtime='1e8',avgscan=T,
coloraxis='baseline' )
We have chosen field='3c279', correlation='XX', and averaged over time (avgtime='1e8', avgscan=T). We plot the different baselines with a different color. If the window is still open, check 'force reload'. This will display the uncorrected phases across the band. To display the corrected phases, you will need to select 'corrected' in the 'Data Column' field of the 'Axes' tab in plotms and re-plot. The amplitude as a function of channel should look flatter with the corrected data. You can also check that the amplitude of the corrected data as a function of time (xaxis='time') is flatter as well, as we have removed part of the atmospheric contribution on the amplitude using the Tsys calibration (for example, due to different elevation of the source with time).
We also recommend to plot plotms with keywords xaxis='channel', yaxis='amp' (this time with avgchannel=3840 instead of avgtime='1e8',avgscan=T) to check the effect of the WVR corrections. We expect smaller phase scatter in the corrected data, as the effect of atmospheric fluctuations on phase (time scales of ~1 sec) is calibrated with the help of WVR.
With the Tsys and WVR calibrations applied successfully and the a priori flagging taken care of we will now split out the corrected data. We will keep only the corrected data, specified via 'datacolumn', and only spectral window 1, which contains FDM (high spectral resolution) observations of the CO(3-2) line. We give the new MS files the extension ".wvrtsys.ms" to indicate that they have been corrected for WVR and Tsys effects. Because split will not overwrite existing files, we remove any previous versions of the new MS before beginning.
# In CASA
for asdm in basename_all:
os.system('rm -rf '+asdm+'.wvrtsys.ms')
split(vis=directory+asdm+'.ms', outputvis=asdm+'.wvrtsys.ms',
datacolumn='corrected', spw='1')
The WVR and Tsys-corrected data now sit in the DATA column of the new measurement sets, which only one spectral window (now labeled spectral window 0 though it was spectral window 1 in the original data). You may wish to run listobs to illustrate the changes:
# In CASA
for asdm in basename_all:
listobs(vis=asdm+'.wvrtsys.ms', listfile=asdm+'.wvrtsys.listobs.txt', verbose=True)
Note the new spectral window information:
Insert clipped spw information from listobs here.
Inspect and Flag Data
Next we do some additional inspection using plotms in order to flag data. We will do this in the time (continuum plot) and spectral domains.
First we will plot amplitude versus time (see Figure 3), averaging over all channels of spectral window 1 (spw = 1, where the CO(3-2) line is) and colorizing by field, for example for the first dataset. Scans on Titan are colored red, the bandpass and phase calibrator 3c279 is colored black, and the different pointings of the Antennae mosaic in different colors. Figure 4 shows the phase versus time plot, with the same color.
Check carefully that the amplitudes and phases vary smoothly with time. The Tsys corrected data should have approximately constant amplitudes and the WVR corrected data should usually have lower phase scatter. Note that the amplitudes in Figure 3 are decreasing, as a result of the decreasing elevation of the source. The Tsys corrected data (choose 'corrected' in Axes > Data Column in plotms) show constant amplitudes.
These are the plotms instances to produce the continuum plots, amplitude and phase versus time:
# In CASA
asdm=basename_all[0]
plotms(vis=asdm+'.ms',
xaxis='time', yaxis='amp',
selectdata=True, spw='1', correlation='XX',antenna='*&*',
avgchannel='3840', avgscan=T,
iteraxis='baseline',coloraxis='field')
plotms(vis=asdm+'.ms',
xaxis='time', yaxis='phase',
selectdata=True, spw='1', correlation='XX',antenna='*&*',
avgchannel='3840', avgscan=T,
iteraxis='baseline',coloraxis='field')
where:
- xaxis='time', yaxis='X' : a plot of X (amplitude or phase) versus time.
- avgchannel='3840' : average over all the channels in the spectral window.
- spw='1', correlation='XX',antenna='*&*': Select only the spectral window 1, polarization XX and cross-correlation data.
- iteraxis='baseline',coloraxis='field': Iterate over baseline, and colorize by different fields.
Select correlation='YY' to inspect the other polarization.
Especially important is to look that there are no sudden phase jumps or decrease in amplitude, larger scatter in phase for the gain calibrator (probably as a result of decorrelation due to weather conditions, which will be more relevant for the longest baselines), and that enough phase calibrator visibilities are available to calibrate the source data. In case you find any of these problems, you may want to flag the data.
Second we plot the amplitude and phase versus frequency for the two correlations, XX and YY. Figure 5 and 6 show an example of spectral plot for the bandpass and phase calibrator 3c279, for both correlations. Again, since this source is a quasar, the amplitudes should be constant (Tsys corrected) and phases varying smoothly with frequency.
These are the task instances to obtain these plots:
# In CASA
plotms(vis = asdm+'.split.ms',
xaxis = 'frequency',yaxis = 'phase',
field='3c279'
avgtime = '1e8',avgscan = T,
selectdata=True, antenna = '*&*',
iteraxis='baseline')
plotms(vis = asdm+'.split.ms',
xaxis = 'frequency',yaxis = 'amp',
field='3c279'
avgtime = '1e8',avgscan = T,
selectdata=True, antenna = '*&*',
iteraxis='baseline')
where:
- xaxis = 'frequency',yaxis = 'X': plot X (amplitude or phase) versus frequency
- field = '3c279': plot only our bandpass and phase calibrator
- avgtime = '1e8',avgscan = T: average all scans and integrations
- antenna = '*&*: plot only cross-correlation data
- iteraxis='baseline': iterate for each baseline
There are no large phase delays found in any of the datasets (i.e., usually less than one wrap over the bandpass), so bandpass calibration should remove this effect properly.
Using these plots we look for channels that would need to be flagged, for example, in the edges of the band, or as a result of spurious interferences.
You can plot other sources as well. By selecting any pointing of Antennae, you should be able to see clearly the (still uncalibrated) CO(3-2) line. The flux calibrator also present some emission line in this spectral window, that will need to be flagged (check the flux calibration subsection).
First we use flagdata to remove the edge channels from both sides of the bandpass:
# In CASA
flagdata(flagbackup = F,vis = asdm+'.split.ms',spw = '1:0~7,1:3831~3839')
Continue to inspect the data with plotms, plotting different axes and colorizing by the different parameters. Don't forget to average the data if possible to speed the plotting process. The time ranges to insert in flagdata can be obtained using plotms Tools Hover/Display. Instead of using the following flagdata commands, you can also flag by hand in plotms. To do this, select your bad data by clicking on the 'Mark Regions" button, then on 'Flag".
File:AntennaeDataInspection-Band7.txt contains the different problems that have been identified for all the datasets. We indicated how to flag the bad data in different instances of the flagdata command. For example, for the first dataset we find that DV04 Tsys is too large in comparison with the other antennas. Also, in the continuum phase plot of all baselines * & PM03, we find that corr=YY for the spw=3 need to be flagged. In the spectral amplitude plot, it is possible to discern some spurious interferences in a few channels. Finally we flag the data with mode='manualflag' and selecting the data to flag. Finally, we save the flag version using flagmanager.
asdm="uid___A002_X1ff7b0_Xb"
flagdata(vis = asdm+'.ms',mode='manualflag',flagbackup = F,antenna='DV04')
flagdata(vis = asdm+'.ms',mode='manualflag',flagbackup = F,
antenna='PM03',correlation='YY',spw='3')
flagdata(vis=asdm+'.ms', flagbackup=F, antenna=['DV12','PM03'],
spw=['0:639~640;1663~1664;2431~2432','1:1146~1147;2182~2184'])
flagmanager(vis =asdm+'.ms',mode = 'save',versionname = 'FlagFinal')
Other minor problems that are reported for this dataset are: 1) in the spectrum phase plot: DV02 10 degrees peak to peak noise in one of the correlations and 2) in the continuum phase plot: *&DV09, there is a sudden change of ~ 200 deg. at 2:03:20, but these can be calibrated.
Bandpass Calibration
Next we plot the phase as a function of time and frequency for the bandpass calibrator, 3c279. For the first plot, Figure 7, we use avgscan=T and avgtime='1e8' to average in time over all scans and integrations, and we specify coloraxis='baseline' to colorize by baseline. For the second, Figure 8, we use spw='0:40~3800' and avgchannel='3840' to average over the central channels of the first spectral window. For both plots we will iterate on antenna (interaxis='antenna'). Use the green arrows of the plotms GUI to view the plots for different antennas.
# In CASA
plotms(vis= asdm+'.split.ms',
xaxis='freq', yaxis='phase',
selectdata=True, field='3c279', corr='XX', antenna='*&*',
avgtime='1e8', avgscan=T,
coloraxis='baseline', iteraxis='antenna')
# In CASA
plotms(vis= asdm+'.split.ms',
xaxis='time', yaxis='phase',
selectdata=True, field='3c279',
spw='0:40~3800', antenna='*&*',corr='XX',
avgchannel='3840', avgscan=T,
coloraxis='baseline', iteraxis='antenna')
Figure 7 shows that the phase varies with time, thus it will need to be gain calibrated. In Figure 8 we see that phase variations as a function of frequency in dataset uid_A002_X1ff7b0_Xb are small, typically ~ 30 degrees, thus bandpass calibration will suffice (i.e. no delay calibration is needed).
We issue gaincal on 3c279 to determine phase(-only) gain solutions. We use solint='int' for the solution interval, which means that one gain solution will be determined for every integration time to prevent de-correlation of the signal. Once phase is corrected, we can determine the bandpass solutions with bandpass. We apply the phase calibration table on-the-fly with the parameter "gaintable". Bandpass response can vary from day to day, therefore we calculate independent bandpass solutions for each dataset.
#In CASA
for asdm in basename_all:
os.system('rm -rf '+asdm+'.bpphase.gcal,'+asdm+'.bandpass.bcal')
gaincal(vis = asdm+'.split.ms',
selectdata=T,field = '3c279',spw = '0:40~3800',
caltable = asdm+'.bpphase.gcal',
solint = 'int',refant = 'DV09',calmode='p')
bandpass(vis = asdm+'.split.ms',
field = '3c279',
gaintable = asdm+'.bpphase.gcal',
caltable = asdm+'.bandpass.bcal',
bandtype='B',
solint = 'inf',combine = 'scan', solnorm=T,refant = 'DV09',
minblperant=3,minsnr=2,fillgaps=1)
where:
- gaintable = asdm+'.bpphase.gcal', caltable = asdm+'.b1.cal': Gain calibration table, and bandpass calibration table
- solint='int' or 'inf': The former is to consider integration by integration. The latter, combined with the default combine='scan', sets the solution interval to the entire observation
- refant = 'DV09': Set the reference antenna to DV06
- calmode='p': Gain cal calibration only phase
- minblperant=3: Minimum number of baselines required per antenna for each solve
- minsnr=2: Minimum SNR for solutions
- bandtype='B': Channel by channel solution for each specified spw
- fillgaps=1: Interpolate channel gaps 1 channel wide
- solnorm=T: Normalize the bandpass amplitudes and phases of the corrections to unity
Do not worry about the message "Insufficient unflagged antennas" when running the bandpass task, which relates to the flagged edge channels.
We plot the time variation of the phase solutions (asdm+'.bpphase.gcal') using plotcal, to check that they vary smoothly with time.
#In CASA
plotcal(caltable = asdm+'.bpphase.gcal',
xaxis = 'time',yaxis = 'phase',
iteration = 'antenna',plotrange=[0,0,-180,180])
We also plot the bandpass solutions with plotcal, and we see that the solutions seem reasonable, with amplitudes close to 1 (Figure 9), and phases that does not vary much over the spectral window.
#In CASA
plotcal(caltable = asdm+'.bandpass.bcal',
xaxis = 'freq',yaxis = 'amp',
plotrange = [0,0,0.8,1.2])
plotcal(caltable = asdm+'.bandpass.bcal',
xaxis = 'freq',yaxis = 'phase',
plotrange = [0,0,-100,100])
Gain (Phase and Amplitude) Calibration
We set the flux for our flux calibrator, Titan, using the task setjy, by applying the Butler-JPL-Horizons 2010 model. First, we check the flux calibrator data. We plot the spectrum and find a bright line emission (Figure 10) in the spectral window. We will flag it using flagdata and update the flag version using flagmanager.
# in CASA
for asdm in basename_all:
flagdata(vis=asdm+'.split.ms',flagbackup=F,
field=['Titan'],
spw=['0:1100~1700'])
flagmanager(vis =asdm+'.split.ms',mode = 'save',versionname = 'FlagFlux')
Second, we do a new gain calibration applying the bandpass calibration solutions on-the-fly.
We solve for amplitude and phase simultaneously and determine average solutions per scan.
# in CASA
for asdm in basename_all:
os.system('rm -rf '+asdm+'.amp.gcal,'+asdm+'.flux.cal')
setjy(vis = asdm+'.split.ms',field = 'Titan',
standard = 'Butler-JPL-Horizons 2010')
gaincal(vis=asdm+'.split.ms',
gaintable=asdm+'.bandpass.bcal',
caltable=asdm+'.intphase.gcal',
calmode='p',
field='Titan,3c279',
spw='0~1:40~3800',
refant='DV09', solint='int',minsnr=2.0,minblperant=4)
gaincal(vis=asdm+'.split.ms',
gaintable=asdm+'.bandpass.bcal',
caltable=asdm+'.scanphase.gcal',
calmode='p',
field='Titan,3c279',
spw='0~1:40~3800',
refant='DV09', solint='inf',minsnr=2.0,minblperant=4)
gaincal(vis = asdm+'.split.ms',
gaintable =[asdm+'.bandpass.bcal',asdm+'.intphase.gcal'],
caltable = asdm+'.amp.cal',
calmode='ap'
field = 'Titan, 3c279',
refant = 'DV09',solint = 'int')
Finally, we will bootstrap the flux density of the secondary calibrator from that of Titan using the task fluxscale.
# in CASA
fluxscale(vis = asdm+'.split.ms',
caltable = asdm+'.amp.gcal',
fluxtable = asdm+'.flux.cal',
reference = 'Titan',
transfer = '3c279')
The flux of Titan at these frequencies is about 2.9 Jy. For example, for dataset uid___A002_X1ff7b0_Xb.f1.cal:
#2011-07-13 07:31:04 INFO setjy Titan spwid= 0 [I=2.846, Q=0, U=0, V=0] Jy
The new flux table asdm+'.flux.gcal' replaces the previous asdm+'.amp.gcal table in future application of the calibration to the data, i.e. the new flux table contains both asdm+'.amp.gcal and the newly acquired flux scaling. Unlike the gain calibration steps, this is not an incremental table.
- gaintable = asdm+'.bandpass.bcal': We apply the bandpass calibration on-the-fly
- caltable = 'asdm+'.amp.cal: the output gain calibration table
- calmode = 'ap': To solve for amplitude and phase
We find that the flux of 3c279 is 10.45 Jy, by averaging the fluxes obtained from the ten available datasets. This flux agree within 10% with the most recent 0.850 mm measurements from the SMA calibrator list [1] : (01 Jul 2011, SMA 9.75 ± 0.49).
Apply the Calibrations and Inspect
Now we will use applycal to apply the bandpass and gaincal tables that we generated in the previous sections. First, we will apply the solutions to the bandpass and secondary calibrator (in this case, it is the same, 3c279). Second to the flux calibrator from the bandpass solutions and its own phase and flux solutions. Finally we will apply the bandpass, phase and amplitude solutions to the science target, the flux calibrator (Titan) and the bandpass calibrator (3c279).
#In CASA
for asdm in basename_all:
applycal(vis=asdm+'.split.ms',field='3c279',
gaintable=[asdm+'bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
interp=['nearest','nearest','nearest'],
gainfield=['3c279','3c279','3c279'],flagbackup=T)
applycal(vis=asdm+'.split.ms',field='Titan',
gaintable=[asdm+'bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
interp=['nearest','nearest','nearest'],
gainfield=['3c279','titan','titan'],flagbackup=T)
applycal(vis=asdm+'.split.ms',field=['Ant*','NGC*']
interp=['nearest','linear','linear'],
gaintable=[asdm+'bandpass.bcal',asdm+'.scanphase.gcal',asdm+'.flux.cal'],
gainfield=['3c279','3c279','3c279'],flagbackup=T)
We plot the corrected amplitudes and phases of 3c279 as a function of time and frequency, to check that the phases are close to zero and the amplitudes are constant.
# In CASA
asdm=basename_all[0]
plotms(vis = asdm+'.split.ms', xaxis='time', yaxis='amp',
ydatacolumn='corrected', selectdata=True, field='3c279',
averagedata=True, avgchannel='3840', avgtime='',
avgscan=F, avgbaseline=F, coloraxis='spw')
plotms(vis = asdm+'.split.ms', xaxis='time', yaxis='pha',
ydatacolumn='corrected', selectdata=True, field='3c279',
averagedata=True, avgchannel='3840', avgtime='',
avgscan=F, avgbaseline=F, coloraxis='spw')
plotms(vis = asdm+'.split.ms', xaxis='time', yaxis='pha',
ydatacolumn='corrected', selectdata=True, field='3c279',
averagedata=True, avgchannel='', avgtime='1e6',
avgscan=T, avgbaseline=F, coloraxis='baseline')
plotms(vis = asdm+'.split.ms', xaxis='time', yaxis='pha',
ydatacolumn='corrected', selectdata=True, field='3c279',
averagedata=True, avgchannel='', avgtime='1e6',
avgscan=T, avgbaseline=F, coloraxis='baseline')
In Fig. 11 and 12 we plot phase vs. channel and amp vs. time for 3c279 for the uid___A002_X1ff7b0_Xb dataset.
Finally we can use plotms to examine the corrected amplitude and phase of Antennae galaxies as a function of frequency, just by changing the field keyword to field='NGC*','Antennae*'.
Split and Concatenate Data for Northern and Southern Mosaics
The individual data sets are now calibrated. We can safely split out the calibrated data for our science target and drop the calibrators. As we do so, we will smooth the data in frequency by setting, averaging together groups of 10 channels by setting width=10 in split. The new data will have a channel width corresponding to about ~4.5 km/s. This will make the imaging steps much more tractable.
#In CASA
basename_all=["uid___A002_X1ff7b0_Xb","uid___A002_X207fe4_X3a","uid___A002_X207fe4_X3b9",
"uid___A002_X2181fb_X49","uid___A002_X1ff7b0_X1c8","uid___A002_X207fe4_X1f7",
"uid___A002_X207fe4_X4d7","uid___A002_X215db8_X18","uid___A002_X215db8_X1d5",
"uid___A002_X215db8_X392"]
for asdm in basename_all:
os.system('rm -rf '+asdm+'.cal.ms')
split(vis = asdm+'.wvrtsys.ms',outputvis = asdm+'.cal.ms',
field = ['NGC*','Antennae*'],spw='0',width=10)
listobs(asdm+'.cal.ms',listfile=asdm+'.cal.listobs.txt')
For convenience we concatenate all data for the Northern Mosaic into a single big MS and place all data for the Southern Mosaic into another file. To do this, we construct a list that holds the names of all the Southern Mosaic MS files and another that holds the name of all the Northern Mosaic MS files then feed these into the concat task.
# In CASA
basename_all_North=["uid___A002_X1ff7b0_Xb","uid___A002_X207fe4_X3a","uid___A002_X207fe4_X3b9","uid___A002_X2181fb_X49"]
basename_all_South=["uid___A002_X1ff7b0_X1c8","uid___A002_X207fe4_X1f7","uid___A002_X207fe4_X4d7",
"uid___A002_X215db8_X18","uid___A002_X215db8_X1d5","uid___A002_X215db8_X392"]
cal_south_vis = [vis+'.cal.ms' for vis in basename_south]
cal_north_vis = [vis+'.cal.ms' for vis in basename_north]
os.system('rm -rf Antennae_South.cal.ms')
concat(vis=cal_south_vis, concatvis='Antennae_South.cal.ms', timesort=T)
os.system('rm -rf Antennae_North.cal.ms')
concat(vis=cal_north_vis, concatvis='Antennae_North.cal.ms', timesort=T)
The syntax used to construct the 'cal_south_vis' variable loops over basename_south and makes a list after adding '.cal.ms' to each member. To see the list 'print cal_south_vis'.
Continue on to Imaging of the Science Target
Now you can continue on to the imaging guide.
Daniel Espada 12:00 UT, 27 July 2011
