SMA CO Line Data 3.1

Overview

This script describes SMA data reduction using miriad to fill and apply Tsys correction, and then writing uvfits. CASA is then used to read uvfits, flag, calibrate, continuum subtract, and image the data. Plotting for datasets with many spw such as SMA data (plotxy) can be slow, where possible the new plotter plotms has been used. The data themselves are a six-pointing mosaic of the infrared dark cloud G19.3+0.07 with 12CO(2-1) in the USB from Brogan et al. in prep. Please do not use these data for scientific purposes.

A tarfile of the uvfits data can be found here. An integrated intensity map of the SMA data overlaid on a Spitzer image can be found here.

Creating the measurement set from uvfits data

# In CASA
# Remove any files from previous run-throughs
os.system('rm -rf spw_fits.files')
os.system('tar xvf spw_fits_files.tar')
os.system('rm -rf g19_d2usb*')

Read in each uvfits file (need to put upper range to one more than have, to satisfy ipython peculiarity). The myfiles makes a list of the filenames that can be used in a later call.

os.system('rm -rf spw_ms_files')
os.system('mkdir spw_ms_files')
myfiles = []
for i in range(1,25):
        msfile = "spw_ms_files/g19_d2usb.spw"+str(i)+".ms"
        importuvfits(fitsfile="spw_fits_files/g19_d2usb.spw"+str(i),
                     vis=msfile)
        myfiles.append(msfile)

Create a single ms from the filenames listed in 'myfile':

concat(vis=myfiles,concatvis='g19_d2usb.ms',timesort=True)

Finally, initialize the scratch columns using clearcal and use listobs to look at the properties in the concatenated file:

clearcal(vis='g19_d2usb.ms')
 
listobs(vis='g19_d2usb.ms',verbose=False)

Put screengrab here from listobs.

Brief description of the purposes of the sources in this project. () gives source id number.

Bandpass calibrator: 3c454.3 (9) Gain calibrators: 1743-038 (1) and 1908-201 (6) Flux calibrators: 3c454.3 (9) and Uranus (13) Science target: g19 (2,3,4,5,7,8) six-pointing mosaic Other sources in the ms were not part of this program.

Initial inspection and flagging

Use plotms to inspect the data. Screengrab here.

clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1~9,13',
       avgchannel='3072')

Now you can use the locate features in plotms to find the bad data points on field 8 and antenna ids 1-2 and 0-2 and either flag in plotms or for script purposes write a flagdata command. Here we will do all flagging within the script.

NOTE: locate gives antenna ids, but flag data looks first at the name column, since SMA antenna names are numbers this can be confusing. The corresponding NAMES are antennas 2-3 and 1-3, which you can check with listobs.

flagdata(vis='g19_d2usb.ms',mode='manualflag',field='8',spw='20~23',
	 timerange='07:52:45~07:52:47',antenna='2&3;1&3')

Plotting again shows that the bad data are gone: (screen grab)

clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='8',
       avgchannel='3072')

NOTE: to confirm the data were flagged in plotms, you need to change selection in plot (for example change field and then change back), and then change back to original to clear the cache and force a re-read of the data.

Make an inital amplitude vs. channel plot of the bandpass calibrator 3c454.3 averaging all the data together in time. crossscan=T needed because there are two scans on 3c454.3. Note the rolloff of 24 individual spw. Try zooming to see better. (Screengrab here)

NOTE: after you run plotxy, run clearstat() before trying to run plotms or you will get a table lock.

clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="data",
       iteration="baseline",
       antenna="",spw="",field="9",
       averagemode="vector",width="1",timebin="all",
       crossscans=T,crossbls=False,stackspw=F,
       subplot=331,plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

As you can see in this plot, SMA spw overlap by a fair amount. Here we will flag the first seven channels on either side of each spw to aid calibration.

default('flagdata')
flagdata(vis='g19_d2usb.ms', mode='manualflag',spw='0~23:0~6;121~127')

On many baselines the first integration of every scan is low. These can also be flagged with flagdata. The integration time of these data is 30.9s. Note that more recent SMA data rarely shows this problem.

flagdata(vis="g19_d2usb.ms",mode="quack",field='',
	 spw="0~23",quackinterval=40.0)

The following can be used to plot the antenna positions. For these observations, antenna 8 was in the barn for repair.

clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="x")

Bandpass Calibration

Next plot phase vs. channel. Note the large descrepancy in phase for spw=4-7 on some baselines, this will be important later.

clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="phase",datacolumn="data",
       iteration="baseline",
       antenna="",spw="",field="9",
       averagemode="vector",width="1",timebin="all",
       crossscans=T,crossbls=False,stackspw=F,
       subplot=331,plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

Now plot amplitude and phase as a function of time.

clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9',
       avgchannel='3072')

The gap in 3C454.3 observation is due to a pause to determine pointing solutions. The flags below catch times just before and after when the amplitudes and phase are off. This flagging could also have been accomplished interactively.

default('flagdata')
flagdata(vis="g19_d2usb.ms",field="9",timerange='12:16:58~12:19:00')
flagdata(vis="g19_d2usb.ms",field="9",timerange='12:08:00~12:09:30')

Here because msplot cannot do multipanel iteration yet we use the slower plotxy.

clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="data",
       iteration="baseline", antenna="",spw="0~23",field="9",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

Start with an initial phase only solution for 3c454.3 to take out major phase variations with time. For maximum sensitivity, we average all the spw together, with the exception of spw 4-7 which show large phase offsets from the other spw in phase vs. channel.

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal',
	field='9', spw='0~3,8~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='p',solint='int', combine='spw')

clearstat()
plotcal(caltable="g19_d2usb.ms.pcal",xaxis="time",yaxis="phase",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Note: A regular 'B' type solution can work well if you have lots of S/N on your bandpass calibrator. You may wish to make this table and compare to 'BPOLY' solution further down. **Can I make the code in grey here too?

bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpcal',
 	field='9', spw='0~23', bandtype='B', solint='inf', combine='scan', 
 	refant='3', gaincurve=False, opacity=0.0,
 	gaintable='g19_d2usb.ms.pcal',spwmap=[0]) 

plotcal(caltable="g19_d2usb.ms.bpcal",xaxis="freq",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=311,
	overplot=False,clearpanel="Auto",iteration="antenna",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

A polynomial solution is more sensitive than simply fitting a polynomial to the channel based solutions 'B' after the fact. For the mm/submm where the S/N is often low on the bandpass calibrator 'BPOLY' is often the best choice, but every dataset is different. You may especially need to play with degamp and degphase. The combine='scan' is needed to average together the two 3c454.3 scans. The spwmap=[0] is telling it to apply the phase solution that was found from the average of all spw to each spw before determining the bandpass solution.

bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpoly7',
	field='9', spw='0~23', bandtype='BPOLY', 
	solint='inf', combine='scan', 
	refant='3', degamp=7, degphase=7,
	gaintable='g19_d2usb.ms.pcal',spwmap=[0])

clearstat()
plotcal(caltable="g19_d2usb.ms.bpoly7",xaxis="freq",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=311,
	overplot=False,clearpanel="Auto",iteration="antenna",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Unfortunately, this kind of plot with many spw and colors takes a long time. Be patient for now its being worked on. You can safely ignore where the solution gets large at edges of spw. This is showing where we flagged channels - this is being worked on. Zoom in if necessary to see that inner portions are nice and flat.

The next step is to apply the bandpass solution while re-calculating the phase only solution for 3c454.3 on a short (integration time) time interval, but now including all spw (including 4~7 since antenna based solutions will have taken out phase vs.spw differences).

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal2',
	field='9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='p',solint='int', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7'])

clearstat()
plotcal(caltable="g19_d2usb.ms.pcal2",xaxis="time",yaxis="phase",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Next we apply the antenna based bandpass, and phase only solution to solve for amplitude solutions over a longer solint (to increase S/N).The spwmap=[[0],[]] tells it that the pcal solution is independent of spw while the bpoly7 solution is per spw, spwmap order must be in the same order as the tables in gaintable.

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.apcal',
	field='9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='a',solint='300s', combine='spw',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.bpoly7'],
	spwmap=[[0],[]])

clearstat()
plotcal(caltable="g19_d2usb.ms.apcal",xaxis="time",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Although in previous steps solutions were applied "on-the-fly", the data were not actually changed. Running applycal will fill the corrected data column in the ms. Again be careful with specification of gaintable, and keeping the same order in spwmap and gainfield. NOTE this applycal is not necessary for later steps -- it is being done in order to allow plotting of data calibration to this point.

applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7'],
	spwmap=[[0],[0],[]],gainfield=['9','9','9'])

clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
       datacolumn="corrected",iteration="baseline",
       antenna='',spw="0~23",field="9",
       averagemode="vector",width="1",timebin="all",crossscans=T,
       crossbls=False,stackspw=False,extendflag=F,
       extendchan="",extendspw="",extendant="",extendtime="",
       plotcolor='darkcyan',subplot=331,
       multicolor="none",figfile="")

As this plotxy shows, after antenna based bandpass, SMA data often show residual baseline based problems as a function of spw (i.e. a few abnormally low amplitudes for a few spw on a few baselines). To fix this problem we need to do a baseline based amplitude only solution per spw, but we'll need to remove an overall normalization factor first by combining all spw (this option will be incorporated into blcal one day so that only one run is necessary).

blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal1',
	field='9', spw='0~23',solint='inf', combine='spw,scan', 
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7'], 
	calmode='a',spwmap=[[0],[0],[]])

blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal2',
	field='9', spw='0~23',solint='inf', combine='scan', 
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.blcal1','g19_d2usb.ms.bpoly7',], 
	calmode='a',spwmap=[[0],[0],[0],[]])

clearstat()
plotcal(caltable="g19_d2usb.ms.blcal1",xaxis="freq",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=511,
	overplot=False,clearpanel="Auto",iteration="antenna",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Now we can do a new applycal including ONLY the second blcal2 solution. Note that applycal overwrites the corrected column so there is no need to clearcal. Like the antenna based bandpass the blcal solution is independent for each spw, so its spwmap entry is [].

applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[0],[0],[],[]],gainfield=['9','9','9','9'])

Note how flat the full bandpass is after all this calibration

clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
       datacolumn="corrected",iteration="baseline",
       antenna='',spw="0~23",field="9",
       averagemode="vector",width="1",timebin="all",crossscans=T,
       crossbls=False,stackspw=False,extendflag=F,
       extendchan="",extendspw="",extendant="",extendtime="",
       plotcolor='darkcyan',subplot=331,
       multicolor="none",figfile="")

Inspection of gain calibrators and setting the flux scale

Apply antenna and baseline based bandpass solutions to the gain calibrators:

clearstat()
applycal(vis='g19_d2usb.ms',spw='0~23', field='1,6',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]],gainfield=['9','9'])

Have a look at both calibrators with channel vs. amplitude:

# First field 1:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="",field="1",
       averagemode="vector",width="1",timebin="all",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Next field 6:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="",field="6",
       averagemode="vector",width="1",timebin="all",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

Have a look at both calibrators in phase and amplitude vs. time. For the moment, unfortunately multicolor='field' doesn't work, so this is easiest to do separately.

# First, amplitude vs. time for field 1: 
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="1",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Second, amplitude vs. time for field 6:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="6",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="red",
       overplot=F,showflags=False,interactive=True,figfile="")

# Third, phase vs. time for field 1:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="1",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Last, phase vs. time for field 6:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="6",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

At the time of these observations, the calibrator J1911 had the most recent and relatively stable monitor data (SMA website), so we will use it for absolute flux calibration. In future tools for using resolved planets will be implemented. ** How to set up nice table format? 1911-201 1mm 17 Jul 2005 11:38 SMA 274.5 1.88 +/- 0.13 mgurwell 1mm 05 Sep 2005 08:10 SMA 221.4 1.99 +/- 0.12 mgurwell

setjy(vis='g19_d2usb.ms', field='6', spw='0~23', fluxdensity=[2.0,0.,0.,0.])

Final gain and flux calibration for all calibrators

To avoid decorrelation of the amplitude, the phase solutions fed to the amplitude calibration step should be at the shortest solint your S/N can support, typically the integration time.

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcal',
	field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='p',solint='int', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]])

Unless you chose to do some more sophisticated smoothing, the phase solutions applied to your target might as well be on the scan time of your phase calibrator(s). Without combine='scan', the solution will not cross scan boundaries even though solint='inf'. This is an easy way to "manually" interpolate.

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcalscan',
	field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='p',solint='inf', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]])

Plot the phase solutions of the phase calibrators. Can repeat for both short and long solution intervals.

clearstat()
plotcal(caltable="g19_d2usb.ms.allpcalscan",xaxis="time",yaxis="phase",
	field="1,6",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")