NGC3256 Band3 Calibration for CASA 3.4

Overview

This portion of the NGC3256Band3 CASA Guide will cover the calibration of the raw visibility data. To skip to the imaging portion of the guide, see: NGC3256 Band3 - Imaging.

If you haven't already downloaded the raw data, you may do that now by clicking on the region closest to your location and downloading the file named 'NGC3256_Band3_UnCalibratedMSandTablesForReduction.tgz':

North America

Europe

East Asia

Once the download has finished, unpack the file:

# In a terminal outside CASA
tar -xvzf NGC3256_Band3_UnCalibratedMSandTablesForReduction.tgz

cd NGC3256_Band3_UnCalibratedMSandTablesForReduction

# Start CASA
casapy

The data have already been converted to CASA Measurement Set (MS) format using the CASA task importasdm. Accompanying the data are some basic calibration tables you will need for the following reduction, as well as the *.ms.flagversions files that are automatically generated by importasdm.

Initial Inspection and A priori Flagging

We will eventually concatenate the six datasets used here into one large dataset. However, we will keep them separate for now, as some of the steps to follow require individual datasets (specifically, the application of the Tsys and WVR tables). We therefore start by defining an array "basename" that includes the names of the six files in chronological order. This will simplify the following steps by allowing us to loop through the files using a simple for-loop in python. Remember that if you log out of CASA, you will have to re-issue this command. We will remind you of this in the relevant sections by repeating the command at the start.

# In CASA
basename=['uid___A002_X1d54a1_X5','uid___A002_X1d54a1_X174','uid___A002_X1d54a1_X2e3',
'uid___A002_X1d5a20_X5','uid___A002_X1d5a20_X174','uid___A002_X1d5a20_X330']

The usual first step is then to get some basic information about the data. We do this using the task listobs, which will output a detailed summary of each dataset supplied.

# In CASA
for name in basename:
        listobs(vis=name+'.ms')

Note that after cutting and pasting a for-loop you often have to press return several times to execute. The output will be sent to the CASA logger. You will have to scroll up to see the individual output for each of the six datasets. Here is an example of the most relevant output for the first file in the list.

	Fields: 3
	  ID   Code Name                RA              Decl          Epoch   SrcId nRows  
	  0    none 1037-295            10:37:16.07900 -29.34.02.8130 J2000   0     38759  
	  1    none Titan               00:00:00.00000 +00.00.00.0000 J2000   1     16016  
	  2    none NGC3256             10:27:51.60000 -43.54.18.0000 J2000   2     151249 
	Spectral Windows:  (9 unique spectral windows and 2 unique polarization setups)
	  SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)  TotBW(kHz)  Corrs  
	  0           4 TOPO  184550      1500000       7500000     I   
	  1         128 TOPO  113211.988  15625         2000000     XX  YY  
	  2           1 TOPO  114188.55   1796875       1796875     XX  YY  
	  3         128 TOPO  111450.813  15625         2000000     XX  YY  
	  4           1 TOPO  112427.375  1796875       1796875     XX  YY  
	  5         128 TOPO  101506.187  15625         2000000     XX  YY  
	  6           1 TOPO  100498.375  1796875       1796875     XX  YY  
	  7         128 TOPO  103050.863  15625         2000000     XX  YY  
	  8           1 TOPO  102043.05   1796875       1796875     XX  YY  
	Sources: 48
	  ID   Name                SpwId RestFreq(MHz)  SysVel(km/s) 
	  0    1037-295            0     -              -            
	  0    1037-295            9     -              -            
	  0    1037-295            10    -              -            
	  0    1037-295            11    -              -            
	  0    1037-295            12    -              -            
	  0    1037-295            13    -              -            
	  0    1037-295            14    -              -            
	  0    1037-295            15    -              -            
	  0    1037-295            1     -              -            
	  0    1037-295            2     -              -            
	  0    1037-295            3     -              -            
	  0    1037-295            4     -              -            
	  0    1037-295            5     -              -            
	  0    1037-295            6     -              -            
	  0    1037-295            7     -              -            
	  0    1037-295            8     -              -            
	  1    Titan               0     -              -            
	  1    Titan               9     -              -            
	  1    Titan               10    -              -            
	  1    Titan               11    -              -            
	  1    Titan               12    -              -            
	  1    Titan               13    -              -            
	  1    Titan               14    -              -            
	  1    Titan               15    -              -            
	  1    Titan               1     -              -            
	  1    Titan               2     -              -            
	  1    Titan               3     -              -            
	  1    Titan               4     -              -            
	  1    Titan               5     -              -            
	  1    Titan               6     -              -            
	  1    Titan               7     -              -            
	  1    Titan               8     -              -            
	  2    NGC3256             0     -              -            
	  2    NGC3256             9     -              -            
	  2    NGC3256             10    -              -            
	  2    NGC3256             11    -              -            
	  2    NGC3256             12    -              -            
	  2    NGC3256             13    -              -            
	  2    NGC3256             14    -              -            
	  2    NGC3256             15    -              -            
	  2    NGC3256             1     -              -            
	  2    NGC3256             2     -              -            
	  2    NGC3256             3     -              -            
	  2    NGC3256             4     -              -            
	  2    NGC3256             5     -              -            
	  2    NGC3256             6     -              -            
	  2    NGC3256             7     -              -            
	  2    NGC3256             8     -              -            
	Antennas: 7:
	  ID   Name  Station   Diam.    Long.         Lat.         
	  0    DV04  J505      12.0 m   -067.45.18.0  -22.53.22.8  
	  1    DV06  T704      12.0 m   -067.45.16.2  -22.53.22.1  
	  2    DV07  J510      12.0 m   -067.45.17.8  -22.53.23.5  
	  3    DV08  T703      12.0 m   -067.45.16.2  -22.53.23.9  
	  4    DV09  N602      12.0 m   -067.45.17.4  -22.53.22.3  
	  5    PM02  T701      12.0 m   -067.45.18.8  -22.53.22.2  
	  6    PM03  J504      12.0 m   -067.45.17.0  -22.53.23.0  

This output shows that three fields were observed: 1037-295, Titan, and NGC3256. Field 0 (1037-295) will serve as the gain calibrator and bandpass calibrator; field 1 (Titan) will serve as the flux calibrator; and field 2 (NGC3256) is, of course, the science target.

Note that there are more than four SpwIDs even though the observations were set up to have four spectral windows. The spectral line data themselves are found in spectral windows 1,3,5,7, which have 128 channels each. The first one (spw 1) is centered on the CO(1-0) emission line in the galaxy NGC 3256 and is our highest frequency spectral window. There is one additional spectral window (spw 3) in the Upper Side Band (USB), and there are two spectral windows (spw 5 and 7) in the Lower Side Band (LSB). These additional spectral windows are used to measure the continuum emission in the galaxy, and may contain other emission lines as well.

Spectral windows 2,4,6,8 contain channel averages of the data in spectral windows 1,3,5,7, respectively. These are not useful for the offline data reduction. Spectral window 0 contains the WVR data. You may notice that there are additional SpwIDs listed in the "Sources" section which are not listed in the "Spectral Windows" section. These spectral windows are reserved for the WVRs of each antenna (seven in our case). At the moment, all WVRs point to spw 0, which contains nominal frequencies. The additional spectral windows (spw 9-15) are therefore not used and can be ignored.

Another important thing to note is that the position of Titan is listed as 00:00:00.0000 +00.00.00.0000. This is due to the fact that for ephemeris objects, the positions are currently not stored in the asdm. This will be handled correctly in the near future, but at present, we have to fix this offline. We will correct the coordinates below by running the procedure fixplanets, which takes the position from the pointing table.

The final column of the listobs output in the logger (not shown above) gives the scan intent. This information is used later to flag the pointing scans and the hot and ambient load calibration scans, using scan intent as a selection option. Also these intents will be used in the future for pipeline processing.

Seven antennas were used for the dataset listed above. Note that numbering in python always begins with "0", so the antennas have IDs 0-6. To see what the antenna configuration looked like at the time of the this observation, we will use the task plotants.

plotants output

# In CASA
plotants(vis=basename[0]+'.ms', figfile=basename[0]+'_plotants.png')

This will plot the antenna configuration on your screen as well as save it under the specified filename for future reference. This will be important later on when we need to choose a reference antenna, since the reference antenna should be close to the center of the array (as well as stable and present for the entire observation).

If you repeat the plotants command for the other five datasets, you will see that there is an additional antenna (DV10) present on the second day of observations. Other than that, the configuration stays constant during the course of the observations.

Flagging

The first editing we will do is some a priori flagging with tflagdata. We will start by flagging the shadowed data and the autocorrelation data. ALMA data contains both the cross correlation and autocorrelation data, but here we are only interested in the cross-correlation data. Additionally, for compact configurations of the array, one antenna can shadow another, blocking its view. These data also need to be flagged.

There are a number of scans in the data that were used by the online system for pointing calibration. These scans are no longer needed, and we can flag them easily with tflagdata by selecting on 'intent'.

Remember that you first need to redefine the "basename" array if you logged out of CASA prior to starting this section.

# In CASA
basename=['uid___A002_X1d54a1_X5','uid___A002_X1d54a1_X174','uid___A002_X1d54a1_X2e3',
'uid___A002_X1d5a20_X5','uid___A002_X1d5a20_X174','uid___A002_X1d5a20_X330']

Now we will loop over the datasets, running the flagging commands:

# In CASA
for name in basename:
	tflagdata(vis=name+'.ms', mode='shadow', applyflags=T, flagbackup = F, async=False)

for name in basename:
        default(tflagdata)
        vis=name+'.ms'
	applyflags=T
	flagbackup = T
	async=False
        mode='manual'
	intent='*POINTING*'
        tflagdata()
	intent='*ATMOSPHERE*'
        tflagdata()
	intent=''
        antenna='*&&&'
        tflagdata()

We will then store the current flagging state for each dataset using the flagmanager:

# In CASA
for name in basename:
        flagmanager(vis = name+'.ms', mode = 'save', versionname = 'Apriori')

Delay Correction (mainly for Antenna DV07)

Due to an issue with antenna DV07 during the commissioning period when these data were taken, it shows large delays in phase for the first three datasets. While the bandpass calibration will attempt to fit and remove small phase delays (i.e., less than one wrap over the bandpass), large delays like those seen here will result in failed solutions. If we want to salvage the data for this antenna, we therefore need to correct the delays by calculating a K-type delay calibration table with gencal. We emphasize that this is not usually a part of the typical calibration procedure, but it may be useful to the reader to see how such a correction is made.

Remember that you first need to redefine the "basename" array if you logged out of CASA prior to starting this section.

# In CASA
basename=['uid___A002_X1d54a1_X5','uid___A002_X1d54a1_X174','uid___A002_X1d54a1_X2e3',
'uid___A002_X1d5a20_X5','uid___A002_X1d5a20_X174','uid___A002_X1d5a20_X330']

DV07: Phase vs Channel

You can see the phase delays by plotting phase versus channel in plotms, as we do here for a single spw, correlation, and baseline with DV07:

# In CASA
plotms(vis=basename[0]+'.ms', xaxis='channel', yaxis='phase',
        spw='3', antenna='PM03&DV07', correlation='XX', avgtime='1e8')

We run gaincal on all datasets, which yields a delay correction for all antennas. This will remove the phase wrapping and seen on baselines incolving DV07. and will aslo serve as a first phase calibration for all baselines. As the gaincal task will not overwrite existing tables, the script starts by deleting any existing versions of the calibration tables with the same name.

DV07: Phase vs frequency before delay correction

DV07: Phase vs frequency after delay correction

# In CASA
for name in basename:
	os.system('rm -rf cal-'+name+'_del.K')
	gaincal(vis=name+'.ms', caltable='cal-'+name+'_del.K',
		field="1037*",spw="1:10~120,3:10~120,5:10~120,7:10~120",
		solint="inf",combine="scan",refant="DV04",
		gaintype="K")

We will apply these K tables to the data in the next section together with the WVR and Tsys correction. To visualize the effect that the delay correction has you may apply the calibration, and then plot first the uncorrected data, and then the corrected data (by selecting 'corrected' in the axes - Y axis - data column field in the plotms GUI). The phase wrapping apparent on a number of baselines should have dissapeared, and the phases should not show large drifts anymore:

# In CASA
for name in basename:
         applycal(vis=name+'.ms', flagbackup=F, spw='1,3,5,7',
         interp='nearest', gaintable='cal-'+name+'_del.K')

plotms(vis=basename[0]+'.ms', xaxis='freq', yaxis='phase',
         spw='1,3,5,7', antenna='', correlation='XX', avgtime='1e8',
         coloraxis='baseline', avgscan=T, selectdata=T, field='1037*')

Tsys calibration and WVR Correction

In the next step the data will be put on a proper flux scale using Tsys measurements and corrected for path length variations due to water vapour fluctuations. The Tsys measurements correct for the atmospheric opacity (to first-order) and allow the calibration sources to be measured at elevations that differ from the science target.

While a set of Tsys and water vapour calibration tables are delivered together with the data, these will only work in CASA version 3.3 and earlier. For CASA 3.4 and higher, these tables have to be recreated.

Create (and inspect) Tsys tables

We now create and inspect the Tsys tables:

# In CASA
for name in basename:
     os.system('rm -rf cal-tsys_'+name)
     gencal(
         vis=name+'.ms',
         caltype='tsys', 
         caltable='cal-tsys_'+name)

Inspect the tables:

# In CASA
for spw in ['1','3','5','7']:
    for name in basename:
        plotcal(caltable='cal-tsys_'+name, xaxis='freq', yaxis='amp',
                spw=spw, subplot=421, overplot=False,
                iteration='antenna', plotrange=[0, 0, 40, 180], plotsymbol='.',
                figfile='cal-tsys_per_spw_'+spw+'_'+name+'.png')

Note that we only plot the spectral windows that contain the spectral line data. In addition to plotting on your screen, the above command will also produce a plot file (png) for each of the datasets and spectral windows. An example plot is shown to the right for uid___A002_X1d54a1_X174.ms.

Upon examination of these plots, we see that for uid___A002_X1d54a1_X174.ms (measurement set 1), there is an outlying feature in spw=7, antenna DV04. This feature occurs during scans 5 and 9, so we flag those data with the following command:

# In CASA
tflagdata(vis='uid___A002_X1d54a1_X174.ms',
	  applyflags=T, flagbackup = T,
	  mode='manual', antenna='DV04', scan='5,9', spw='7')

Aside from the large amplitudes in the edge channels (which we will handle below), the plots look acceptable. Note that in the lowest spectral window ID (spw 1), Tsys rises toward higher frequencies. This is due to a spectral line from O2 at 117 GHz.

Create water vapour correction tables

To create the wvr calibration tables:

# In CASA
for name in basename:
	os.system('rm -rf cal-wvr_'+name)
	wvrgcal(vis=name+'.ms', caltable='cal-wvr-'+name,  toffset=-1, segsource=True,
		tie=["Titan,1037-295,NGC3256"], statsource="1037-295")

Apply Tsys, water vapour, and delay correction

We will apply the Tsys and wvr tables with applycal, together with the delay correction. We do this for each field separately so that the appropriate calibration data are applied to the right fields. The "field" parameter specifies the field to which we will apply the calibration, and the "gainfield" parameter specifies the field from which we wish to take the calibration solutions from the gaintable. In the call to applycal, we will specify interpolation="nearest". This is important particularly for the WVR corrections; it doesn't make a difference for the delay corrections because they have no time dependence.

# In CASA
for name in basename:
	for field in ['Titan','1037*','NGC*']:
		applycal(vis=name+'.ms', spw='1,3,5,7', flagbackup=F, 
			 field=field, gainfield=field,
			 interp=['nearest','nearest','nearest'], 
			 gaintable=['cal-tsys_'+name,'cal-wvr-'+name,'cal-'+name+'_del.K'])

Without giving the full recipe here, we suggest at this point that you use plotms to plot channel-averaged amplitudes as a function of time, comparing the DATA and CORRECTED columns after applying the Tsys correction. This way you can check that calibration has done what was expected, which is put the data onto the Kelvin temperature scale.

Now we split out the CORRECTED data column of the datasets with the task split, retaining only spectral windows 1,3,5,7. This will get rid of the extraneous spectral windows, including the channel averaged spectral windows and spw 0, which is the one that contained the WVR data. We give the resulting datasets the extension "_line" to indicate that they only contain the spectral windows with the spectral line data.

# In CASA
for name in basename:
	os.system('rm -rf '+name+'_line.ms*')
	split(vis=name+'.ms', outputvis=name+'_line.ms',
		datacolumn='corrected', spw='1,3,5,7')

Additional Data Inspection

Flagging: spectral window edges

We will now do some additional inspection with plotms. First we will plot amplitude versus channel, averaging over time and baselines in order to speed up the plotting process; do this plot for all 6 measurement sets:

(Remember that you first need to redefine the "basename" array if you logged out of CASA prior to starting this section.)

# In CASA
basename=['uid___A002_X1d54a1_X5','uid___A002_X1d54a1_X174','uid___A002_X1d54a1_X2e3',
'uid___A002_X1d5a20_X5','uid___A002_X1d5a20_X174','uid___A002_X1d5a20_X330']

Amplitude vs. channel, averaged over all spw's and all baselines, for measurement set 0 (zoomed to show the low amplitude channels as well)

# In CASA
plotms(vis=basename[0]+'_line.ms', xaxis='channel', yaxis='amp',
       averagedata=T, avgbaseline=T, avgtime='1e8', avgscan=T)

From these plots we see that the edge channels have abnormally high or low amplitudes (you may actually have to zoom in to see the normal amplitudes against the extreme high amplitude outliers). We will use tflagdata to remove the edge channels from both sides of the bandpass:

# In CASA
for name in basename:
	tflagdata(name+'_line.ms', flagbackup=T,  mode='manual', spw='*:0~16')
	tflagdata(name+'_line.ms', flagbackup=T,  mode='manual', spw='*:125~127')

Next, we will look at amplitude versus time, averaging over all channels and colorizing by field (do the following plot for all 6 measurement sets). The first thing to notice is the difference in Titan's amplitude between the two days and the large range the temporal change in amplitude during the second day. The plot on the right shows the first measurement set of the first day of observations, only showing spw 1 in this case. Scans on Titan are colored red, NGC3256 is orange, and the calibrator 1037-295 is colored black. If you select other spws, you can see some outlying points, which will be flagged later on.

Amplitude vs. time for spw 0, averaged over all channels, for measurement set 0

# In CASA
plotms(vis=basename[0]+'_line.ms', xaxis='time', yaxis='amp',
       averagedata=T, avgchannel='128', coloraxis='field',
       iteraxis='spw')

Titan: contamination by Saturns rings, fix position

Titan is our primary flux calibrator. However, for the second day of observations, Titan had moved too close to Saturn, and Saturn's rings moved into the primary beam. Another way to see this is to plot amplitude versus uv-distance:

# In CASA
plotms(vis=basename[0]+'_line.ms',  xaxis='uvdist', yaxis='amp',
       averagedata=T, avgchannel='128', field='Titan',
       iteraxis='spw')

If you do this for all measurement sets, you will notice that for the first three the amplitude is roughly constant over the total range of uv distances, indicating a compact, unresolved source. For the last three measurement sets, the amplitudes are much larger for shorter baselines, indicating the presence of large scale emission, due to Saturns rings. We will therefore need to flag the Titan scans for the second day:

# In CASA
for i in range(3,6): # loop over the last three data sets
	name=basename[i]
	tflagdata(vis = name+'_line.ms', flagbackup = T,
		  mode='manual', field='Titan')

We also find that during the first day, the Titan observations in spw 2 and 3 are also affected by Saturn. These spectral windows are at lower frequencies and therefore correspond to slightly larger primary beams. This results in Saturn just being picked up in spw 2 and 3, but not in spw 0 and 1. We do not flag these data here, but have to take this effect into account when we do the flux calibration later.

Next, we will fix the position of Titan in the combined dataset. Recall that the position of the Titan field is currently set to 00:00:00.0000 +00.00.00.0000. The following procedure will replace this with the actual mean position observed by the telescopes and, at the same time, it will recalculate the uvw coordinates.

If you are running casa revision r15777 or newer, fixplanets() is already defined as a task. However, if you are using an older revision, then you must first initialize this script with this command:

execfile(casadef.python_library_directory+'/fixplanets.py')

Now fix the position of Titan:

# In CASA
for i in range(3): # loop over the first three data sets
	name=basename[i]
	fixplanets(vis=name+'_line.ms', field='Titan', fixuvw=True)

The third parameter in fixplanets, set to True, indicates that the uvw-coordinates for Titan are recalculated. Note that for Cycle 0 data, the coordinates of ephemeris objects will be treated correctly in the data.

Now check to see that the coordinates for Titan have been corrected (e.g., doing a listob for the first measurement set):

2012-04-27 11:19:55 INFO listobs	Fields: 3
2012-04-27 11:19:55 INFO listobs	  ID   Code Name                RA              Decl          Epoch   SrcId nRows  
2012-04-27 11:19:55 INFO listobs	  0    none 1037-295            10:37:16.07900 -29.34.02.8130 J2000   0     8960   
2012-04-27 11:19:55 INFO listobs	  1    none Titan               12:51:24.87318 -02.32.21.3997 J2000   1     3696   
2012-04-27 11:19:55 INFO listobs	  2    none NGC3256             10:27:51.60000 -43.54.18.0000 J2000   2     34944  

some more flagging

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. You will find the following:

• Baselines with DV07 have very high amplitudes in spw 3, correlation YY
• Baselines with DV08 have very low amplitudes in spw 3, correlation YY, but only for the last observation
• Baselines with PM03 have low amplitudes at 2011/04/17/02:15:00 for spw 0
• Baselines with PM03 have low amplitudes at 2011/04/16/04:15:15 for a number of reads in correlation XX

The times to insert in tflagdata can be obtained using plotms Tools Hover/Display. Instead of using the following tflagdata 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".

We flag the bad data with the following commands:

# In CASA
for name in basename:
	tflagdata(vis=name+'_line.ms', flagbackup=T, 
		  mode='manual', spw='3', correlation='YY', antenna='DV07')

tflagdata(vis=basename[5]+'_line.ms', flagbackup=T, 
	  mode='manual', spw='3', correlation='YY', antenna='DV08')

tflagdata(vis=basename[4]+'_line.ms', flagbackup=T, 
        mode='manual', spw='0',
        correlation='', antenna='PM03', timerange='2011/04/17/02:15:00~02:15:32')

tflagdata(vis=basename[1]+'_line.ms', flagbackup=T, 
	  mode='manual', spw='', 
	  correlation='XX', antenna='PM03', timerange='2011/04/16/04:13:50~04:18:00')

Bandpass Calibration

We are now ready to begin the bandpass calibration. First, we will inspect the bandpass calibrator by plotting the phase as a function of frequency and time. For the first plot we use avgscan=T and avgtime='1E6' to average in time over all scans, and we specify coloraxis='baseline' to colorize by baseline. For the second, we use spw='0:30~90' and avgchannel='128' to average over the central 61 channels of the first spectral window. For both plots we will iterate on antenna, so you will need to use the green arrows at the bottom of the plotms GUI to view different antennas.

(Remember that you first need to redefine the "basename" array if you logged out of CASA prior to starting this section (see above).)

Phase vs. frequency for the phase calibrator, averaged over time, for measurement set 1

Phase vs. time for the phase calibrator, all baselines with antenna PM03, measurement set 1

# In CASA
plotms(vis=basename[1]+'_line.ms', xaxis='freq', yaxis='phase', selectdata=True,
	field='1037*', avgtime='1E6', avgscan=T, coloraxis='baseline', iteraxis='antenna')

# In CASA
plotms(vis=basename[1]+'_line.ms', xaxis='time', yaxis='phase', selectdata=True,
	field='1037*', spw='0:30~90', avgchannel='128', avgscan=T, 
        coloraxis='baseline', iteraxis='antenna')

The top plot on the right shows the time-averaged phase as a function of frequency for the calibrator 1037-295, with all baselines shown at once (i.e., without iteraxis='antenna') for presentation purposes. We see that the phase variations across the spectral windows are modest, as expected, as a first order bandpass correction was already done with the delay correction. The second plot shows how the phase varies as function of time. For clarity, we only show baselines with one antenna. There are clearly phase variations on short time scales that we wish to correct for before calculating the bandpass solutions.

Hence, we run gaincal on the bandpass calibrator to determine phase-only gain solutions. We will use solint='int' for the solution interval, which means that one gain solution will be determined for every integration time. This short integration time is possible because the bandpass calibrator is a very bright point source, so we have very high signal-to-noise and a perfect model. This will correct for any phase variations in the bandpass calibrator as a function of time, a step which will prevent decorrelation of the vector-averaged bandpass solutions. We will then apply these solutions on-the-fly when we run bandpass.

We will use the average of channels 40 to 80 to increase our signal-to-noise in the determination of the antenna-based phase solutions. Averaging over a subset of channels near the center of the bandpass is acceptable when the phase variation as a function of channel is small, which it is here. For our reference antenna, we choose DV07. We call the output calibration tables "cal-basename[i].BP.int.p".

# In CASA
for name in basename:
	gaincal(vis=name+'_line.ms', caltable='cal-'+name+'.BP.int.p', 
		spw='*:40~80', field='1037*',
		selectdata=T, solint='int', refant='DV07', calmode='p')

We then check the time variations of the phase solutions with plotcal. We will plot the XX and YY polarization products separately and make different subplots for each of the spectral windows. This is done by setting the "iteration" parameter to "spw" and specifying subplot=221. By setting the parameter "figfile" to a non-blank value, it will also generate png files of the plots.

Phase-only gaincal solutions vs. time for correlation XX, measurement set 0

# In CASA
for name in basename:
	plotcal(caltable='cal-'+name+'.BP.int.p', xaxis = 'time', yaxis = 'phase',
		poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw',
		figfile='cal-'+name+'-phase_vs_time_XX.BP.int.p.png', subplot = 221)

# In CASA
for name in basename:
	plotcal(caltable='cal-'+name+'.BP.int.p', xaxis = 'time', yaxis = 'phase',
		poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw',
		figfile='cal-'+name+'-phase_vs_time_YY.BP.int.p.png', subplot = 221)

If you examine the solutions you will notice that they look very reasonable; in particular, there are no apparent large phase jumps. You may also apply the calibration tables to the datasets and plot phase vs. time for the calibrator and see how the 'data' and 'corrected' columns compare ('corrected' should scatter around 0).

Now that we have a first measurement of the phase variations as a function of time, we can determine the bandpass solutions with bandpass. We will apply the phase calibration table on-the-fly with the parameter "gaintable". Now that the phases are corrected, the data can be time-averaged over longer intervals to maximise SNR in each individual channel. Do not worry about the message "Insufficient unflagged antennas", which relates to the flagged edge channels.

# In CASA
for name in basename:
	bandpass(vis=name+'_line.ms', caltable='cal-'+name+'.BP.inf.ap', 
		 gaintable='cal-'+name+'.BP.int.p', interp='nearest',
		 field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan',
		 bandtype='B', fillgaps=1, refant = 'DV07', solnorm = T)

• caltable='cal-'+name+'.BP.inf.ap' : Output bandpass calibration tables
• gaintable='cal-'+name+'.BP.int.p' : Gain calibration tables to apply (on-the-fly, preliminary phase calibration)
• minblperant=3 : Minimum number of baselines required per antenna for each solve
• minsnr=2 : Minimum SNR for solutions
• solint='inf' : This setting, combined with the default combine='scan', sets the solution interval to the entire observation
• combine='scan' : The solutions cross scans
• bandtype='B' : The default type of bandpass solution, which does a 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

Bandpass phase and amplitude solutions, measurement set 0

We then plot the bandpass solutions with the following commands:

# In CASA
for name in basename:
	plotcal(caltable = 'cal-'+name+'.BP.inf.ap', xaxis='freq', yaxis='phase', spw='',
		subplot=111, overplot=False, plotrange = [0,0,-45,45],
		plotsymbol='.', timerange='',
		figfile='cal-'+name+'-bandpass-phase.BP.inf.ap.png')

# In CASA
for name in basename:
	plotcal(caltable = 'cal-'+name+'.BP.inf.ap', xaxis='freq', yaxis='amp', spw='',
		subplot=111, overplot=False,
		plotsymbol='.', timerange='',
		figfile='cal-'+name+'-bandpass-amp.BP.inf.ap.png')

The solutions seem reasonable, so we will in the following apply them on-the-fly during gain calibration.

Gain Calibration

The first step is to set the flux density for Titan using the task setjy. We will use the Butler-JPL-Horizons 2010 model:

(Remember that you first need to redefine the "basename" array if you logged out of CASA prior to starting this section (see above).)

# In CASA
for name in basename:
	setjy(vis=name+'_line.ms', 
	      field='Titan', standard='Butler-JPL-Horizons 2010', 
	      spw='0,1,2,3')

The flux density of Titan is 370 mJy at 114 GHz (spw 0) and drops to 288 mJy at 101 GHz (spw 2):

2011-12-08 13:33:30     INFO    setjy::imager::setjy()         Titan (fld ind 1) spw 0  [I=0.36975, Q=0, U=0, V=0] Jy, (JPL-Butler Solar System Object)
2011-12-08 13:33:30     INFO    setjy::imager::setjy()         Titan (fld ind 1) spw 1  [I=0.35864, Q=0, U=0, V=0] Jy, (JPL-Butler Solar System Object)
2011-12-08 13:33:30     INFO    setjy::imager::setjy()         Titan (fld ind 1) spw 2  [I=0.28766, Q=0, U=0, V=0] Jy, (JPL-Butler Solar System Object)
2011-12-08 13:33:30     INFO    setjy::imager::setjy()         Titan (fld ind 1) spw 3  [I=0.29642, Q=0, U=0, V=0] Jy, (JPL-Butler Solar System Object)

The gain calibration will be done in three steps. First, short term phase variations will be accounted for by solving for the phases on an integration by integration basis. This is done to minimize the effect of decorrelation in the amplitude calibration due to short term phase fluctuations (which will be particularly important when working at higher frequencies and/or longer baselines).

# In CASA
for name in basename:
	gaincal(vis=name+'_line.ms' , caltable='cal-'+name+'.int.p', 
		spw ='*:16~112', field = '1037*,Titan', minsnr=1.0,
		solint= 'int', selectdata=T, solnorm=False, 
		refant = 'DV07',interp='nearest',
		gaintable = 'cal-'+name+'.BP.inf.ap', calmode = 'p')

• caltable='cal-'+name+'.int.p' : the output gain calibration tables
• minsnr=1.0: To reject solutions with a signal-to-noise less than 1.0
• calmode = 'p' : To solve for phase only
• spw='*:16~112' : to select all spectral windows, but only the inner 75% of the channels of each
• solint='int' : one solution per integration, to take out short term phase fluctuations
• solnorm=False : We do not want to normalize the solutions to unity since we wish to relate the measured amplitudes of the secondary calibrator (1037-295) to the flux calibrator (Titan)
• gaintable = 'cal-'+name+'.BP.inf.ap' : We apply the bandpass calibration on-the-fly

The next step will provide the amplitude calibration solution. We do a new gain calibration, this time applying the bandpass calibration solutions and the short-time phase solution on-the-fly. We solve for amplitude and phase simultaneously and determine average solutions per scan:

# In CASA
for name in basename:
        gaincal(vis=name+'_line.ms' , caltable='cal-'+name+'.inf.ap', 
                spw ='*:16~112', field = '1037*,Titan', minsnr=1.0,
                solint= 'inf', selectdata=T, solnorm=False, 
                refant = 'DV07',interp = ['nearest','nearest'],
                gaintable = ['cal-'+name+'.BP.inf.ap','cal-'+name+'.int.p'], calmode = 'ap')

• caltable='cal-'+name+'.inf.ap' : the output gain calibration tables
• minsnr=1.0: To reject solutions with a signal-to-noise less than 1.0
• calmode = 'ap' : To solve for amplitude and phase
• spw='*:16~112' : to select all spectral windows, but only the inner 75% of the channels of each
• solint='inf' : Together with the default for the "combine" parameter, this setting will solve for one solution per scan
• solnorm=False : We do not want to normalize the solutions to unity since we wish to relate the measured amplitudes of the secondary calibrator (1037-295) to the flux calibrator (Titan)
• gaintable = ['cal-'+name+'.BP.inf.ap','cal-'+name+'.int.p'] : We apply the bandpass and short-term phase calibration on-the-fly

Finally, in the third step an additional phase solution is derived to be used to correct the phases on the science target; this phase solution averages the calibrator phase over the entire calibration scan (using the short-term phase solution derived above, we would only interpolate between corrections derived for the last/first integration of the calibrator scans before/after the science target observation, and would thus be dominated by short term phase scatter instead of long-term phase drifts).

# In CASA
for name in basename:
	gaincal(vis=name+'_line.ms' , caltable='cal-'+name+'.inf.p', 
		spw ='*:16~112', field = '1037*,Titan', minsnr=1.0,
		solint= 'inf', selectdata=T, solnorm=False, 
		refant = 'DV07',interp='nearest',
		gaintable = 'cal-'+name+'.BP.inf.ap', calmode = 'p')

• caltable='cal-'+name+'.inf.p' : the output gain calibration tables
• minsnr=1.0: To reject solutions with a signal-to-noise less than 1.0
• calmode = 'p' : To solve for phase
• spw='*:16~112' : to select all spectral windows, but only the inner 75% of the channels of each
• solint='inf' : Together with the default for the "combine" parameter, this setting will solve for one solution per scan
• solnorm=False : We do not want to normalize the solutions to unity since we wish to relate the measured amplitudes of the secondary calibrator (1037-295) to the flux calibrator (Titan)
• gaintable = 'cal-'+name+'.BP.inf.ap' : We apply the bandpass calibration on-the-fly

Now we will examine the amplitude and phase solutions as a function of time, iterating on spectral window and plotting the XX and YY correlations separately for clarity. Zoom in on the GUI to examine the plots in more detail. You can also use the "Mark Region" and "Locate" buttons on the toolbar to identify points. Use the field parameter to select which calibrator you want to plot the solutions for.

Gain amplitude Solutions, X polarisation, measurement set 0

Gain phase Solutions, X polarisation, measurement set 0

# In CASA
for name in basename:
	plotcal(caltable='cal-'+name+'.inf.p', xaxis = 'time', yaxis = 'phase',
		poln='X', plotsymbol='o', iteration = 'spw',
		figfile='cal-'+name+'-phase_vs_time_XX.inf.p.png', subplot = 221)

# In CASA
for name in basename:
	plotcal(caltable='cal-'+name+'.inf.p', xaxis = 'time', yaxis = 'phase',
		poln='Y', plotsymbol='o', iteration = 'spw',
		figfile='cal-'+name+'-phase_vs_time_YY.inf.p.png', subplot = 221)

# In CASA
for name in basename:
	plotcal(caltable='cal-'+name+'.inf.ap', xaxis = 'time', yaxis = 'phase',
		poln='X', plotsymbol='o', iteration = 'spw',
		figfile='cal-'+name+'-phase_vs_time_XX.inf.ap.png', subplot = 221)

(phase correction for these solutions should be close to zero)

# In CASA
for name in basename:
	plotcal(caltable='cal-'+name+'.inf.ap', xaxis = 'time', yaxis = 'phase',
		poln='Y', plotsymbol='o', iteration = 'spw',
		figfile='cal-'+name+'-phase_vs_time_YY.inf.ap.png', subplot = 221)

(phase correction for these solutions should be close to zero)

# In CASA
for name in basename:
	plotcal(caltable='cal-'+name+'.inf.ap', xaxis = 'time', yaxis = 'amp',
		poln='X', plotsymbol='o', iteration = 'spw',
		figfile='cal-'+name+'-amp_vs_time_XX.inf.ap.png', subplot = 221)

# In CASA
for name in basename:
	plotcal(caltable='cal-'+name+'.inf.ap', xaxis = 'time', yaxis = 'amp',
		poln='Y', plotsymbol='o', iteration = 'spw',
		figfile='cal-'+name+'-amp_vs_time_YY.inf.ap.png', subplot = 221)

We find very stable amplitude solutions for the phase calibrator. The lower points for the first day are the solutions for Titan.

Flux Calibration

We will now bootstrap the flux density of the secondary calibrator from that of Titan using the task fluxscale. The new flux tables cal-+basename[i]+.inf.flux replace the previous cal-+basename[i]+.inf.ap tables in future application of the calibration to the data, i.e. the new flux table contains both cal-+basename[i]+.inf.ap and the newly acquired flux scaling. Unlike the gain calibration steps, this is not an incremental table.

As we previously flagged the Titan scans from the second day due to contamination from nearby Saturn's rings, we will first determine the flux of the secondary calibrator using the Titan scans from the first day, and use this flux value to flux calibrate the second days data. This would not always be an acceptable solution; however, in this case, examination of the measured amplitudes for the phase calibrator over both days indicate that the flux scale is relatively stable.

# In CASA
for i in range(3): # loop over the first three data sets
	name=basename[i]
	fluxscale(vis=name+'_line.ms', caltable='cal-'+name+'.inf.ap',
		  fluxtable='cal-'+name+'.inf.flux', reference="Titan",
		  transfer="1037*", refspwmap=[0,1,1,1])

Also note that we have used refspwmap=[0,1,1,1]. This means that the 1037-295 observations in spw 2 and 3 are referenced to Titan observations observed in spw 1. We have to do this here because the Titan data in spw 2 and 3 are contaminated by Saturn, as explained above.

The logger produces the following output for the three measurement sets:

2011-12-09 10:13:15	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=0 is: 1.85992 +/- 0.0604849 (SNR = 30.7502, N= 14)
2011-12-09 10:13:15	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=1 is: 1.87698 +/- 0.0573121 (SNR = 32.7501, N= 14)
2011-12-09 10:13:15	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=2 (ref SpW=1) is: 1.96112 +/- 0.0561983 (SNR = 34.8964, N= 14)
2011-12-09 10:13:15	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=3 (ref SpW=1) is: 1.87457 +/- 0.0545976 (SNR = 34.3342, N= 12)

2011-12-09 10:13:16	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=0 is: 1.81931 +/- 0.0550219 (SNR = 33.0651, N= 14)
2011-12-09 10:13:16	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=1 is: 1.79592 +/- 0.0478398 (SNR = 37.5404, N= 14)
2011-12-09 10:13:16	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=2 (ref SpW=1) is: 1.87941 +/- 0.0524405 (SNR = 35.839, N= 14)
2011-12-09 10:13:16	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=3 (ref SpW=1) is: 1.88375 +/- 0.05432 (SNR = 34.6788, N= 12)

2011-12-09 10:13:16	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=0 is: 1.86175 +/- 0.0393508 (SNR = 47.3117, N= 14)
2011-12-09 10:13:16	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=1 is: 1.85193 +/- 0.0494996 (SNR = 37.413, N= 14)
2011-12-09 10:13:16	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=2 (ref SpW=1) is: 1.94516 +/- 0.0508357 (SNR = 38.2637, N= 14)
2011-12-09 10:13:16	INFO	fluxscale::::	 Flux density for 1037-295 in SpW=3 (ref SpW=1) is: 1.90879 +/- 0.061465 (SNR = 31.055, N= 12)

We find that the flux density of 1037-295 is ~1.85 Jy in the higher frequency spectral windows and ~1.9 Jy in the lower frequency spectral windows. These fluxes agree very well with the 3mm measurements from the ATCA and OVRO of a few years ago (see SMA calibrator list).

We now set the flux of 1037-295 in the measurement sets of the second day. We use a negative spectral index to get a higher flux in the lower frequency spectral windows (the result of applying setjy can be checked by plotting amplitude vs. frequency and selecting the 'model' data column in the plotms gui). Then we derive a new, flux scaled gain calibration table (again correcting for short-term phase variations):

# In CASA
for i in range(3,6): # loop over the last three data sets
	name=basename[i]
	setjy(vis=name+'_line.ms', 
	      field='1037*', fluxdensity=[1.835], 
	      spix=-0.35, reffreq='115GHz',
	      spw='0,1,2,3')

# In CASA
for i in range(3,6): # loop over the last three data sets
	name=basename[i]
	gaincal(vis=name+'_line.ms' , caltable='cal-'+name+'.inf.flux', 
		spw ='*:16~112', field = '1037*', minsnr=1.0,
		solint= 'inf', selectdata=T, solnorm=False, 
		refant = 'DV07',interp=['nearest','nearest'],
		gaintable = ['cal-'+name+'.BP.inf.ap','cal-'+name+'.int.p'], calmode = 'ap')

Applying the calibrations

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 from the secondary calibrator to the science target and the secondary calibrator itself. We want to use the solutions from the secondary calibrator for the science target because, as the phase calibrator, it was observed throughout the observations. Its flux scale has also now been bootstrapped to that of the flux calibrator, so it can serve as an amplitude and phase calibrator. The application of the solutions to the secondary calibrator itself is a way for us to check the quality of the calibration. Below, we will do the same for the flux calibrator (Titan), applying its amplitude and phase solutions to itself to check the calibration.

# In CASA
for name in basename:
	applycal(vis=name+'_line.ms', flagbackup=F, field='1037*',
	interp=['nearest','nearest','nearest'], gainfield = ['1037*', '1037*', '1037*'],
	gaintable=['cal-'+name+'.BP.inf.ap','cal-'+name+'.int.p','cal-'+name+'.inf.flux'])

# In CASA
for name in basename:
	applycal(vis=name+'_line.ms', flagbackup=F, field='NGC*',
	interp=['nearest','linear','linear'], gainfield = ['1037*', '1037*', '1037*'],
	gaintable=['cal-'+name+'.BP.inf.ap','cal-'+name+'.inf.p','cal-'+name+'.inf.flux'])

• gaintable=['cal-'+name+'.BP.inf.ap','cal-'+name+'.inf.p','cal-'+name+'.inf.flux'] : We apply the bandpass tables, the (scan averaged) phase-gaincal tables, and the gaincal tables with the correct flux scaling
• gainfield = ['1037*', '1037*', '1037*'] : We use the solutions from the secondary calibrator (1037-295) in all three tables
• interp=['nearest','linear','linear']: We opt for interpolation mode "nearest" for the bandpass calibration and (for 1037-295) for phase and amplitude calibration; for NGC3256, we use a linear interpolation for the phase and amplitude calibration.
• flagbackup=F : We do not back up the state of the flags before applying calibration

Check calibrated data, additional flagging

Corrected Amp vs Freq, time average, measurement set 0

Let's now use plotms to examine the corrected amplitude of NGC3256 as a function of frequency:

# In CASA
plotms(vis=basename[0]+'_line.ms', xaxis='freq', yaxis='amp',
	ydatacolumn='corrected', selectdata=True, field='NGC*',
	averagedata=True, avgchannel='', avgtime='10000s',
	avgscan=True, avgbaseline=T, coloraxis='spw')

Use the "Zoom" tool to zoom in on the different spectral windows. (Note that spw 2-3 are on the left). There are three emission lines in spectral windows 0-1. There is no emission in spw 2 and 3. Also note that there are still some high channels on the edge of spw 2-3, particularly in measurement set 5; similarly, if you plot the corrected phase as a function of frequency for 1037-295 you will notice that a few channels at the spectral window edges show noisy phases. We will exclude these channels for our subsequent imaging.

# In CASA
plotms(vis=basename[0]+'_line.ms', xaxis='freq', yaxis='phase',
	ydatacolumn='corrected', selectdata=True, field='1037*',
	averagedata=True, avgchannel='', avgtime='1e8',
	avgscan=True, avgbaseline=F, coloraxis='baseline')

It also makes sense to plot the corrected amplitudes and phases of 1037-295 as a function of time and frequency, to check that the phases are close to zero and the amplitudes are constant. The result of the amplitude versus time plot is shown to the right for measurement set 0. Note that the amplitude differs between spectral windows; this is expected for a source with a spectral index different from 0.

Corrected amplitude vs. time for the calibrator, measurement set 0

Corrected phase vs. time for the calibrator, measurement set 1

# In CASA
plotms(vis=basename[0]+'_line.ms', xaxis='time', yaxis='amp',
	ydatacolumn='corrected', selectdata=True, field='1037*',
	averagedata=True, avgchannel='128', avgtime='',
	avgscan=False, avgbaseline=F, coloraxis='spw')

# In CASA
plotms(vis=basename[1]+'_line.ms', xaxis='time', yaxis='phase',
	ydatacolumn='corrected', selectdata=True, field='1037*',
	averagedata=True, avgchannel='128', avgtime='',
	avgscan=False, avgbaseline=F, coloraxis='spw')

Based on the inspection of these plots, we identify some bad data and flag them UPDATE!:

• in measurement set 1, phases are off for one read at 04:17:09
• in measurement set 2, phases are exceptionally noisy for one read at 05:21:16
• in measurement set 3, phases are very noisy for one read at 00:18:36 and for one read at 00:19:19

# In CASA
tflagdata(vis=basename[1]+'_line.ms',  flagbackup=T,
        mode='manual', selectdata=T, timerange='2011/04/16/04:17:06~04:17:12')

# In CASA
tflagdata(vis=basename[2]+'_line.ms',  flagbackup=T,
        mode='manual', selectdata=T, timerange='2011/04/16/05:21:13~05:21:19')

# In CASA
tflagdata(vis=basename[3]+'_line.ms',  flagbackup=T,
        mode='manual', selectdata=T, timerange='2011/04/17/00:18:33~00:18:39')

# In CASA
tflagdata(vis=basename[3]+'_line.ms',  flagbackup=T,
        mode='manual', selectdata=T, timerange='2011/04/17/00:19:16~00:19:21')