EVLA high frequency Spectral Line tutorial - IRC+10216 CASA3.4
Please use CASA 3.4 for this tutorial
- 1 Overview
- 2 Getting the data
- 3 How to Use This casaguide
- 4 The Observing Log, Antenna Position Corrections, Opacities, Gaincurves, and other Calibration "Priors"
- 5 Initial Inspection and Flagging
- 6 Set Up the Model for the Flux Calibrator
- 7 Bandpass and Delay
- 8 Gain Calibration
- 9 Applycal and Inspect
- 10 Redo Calibration after more Flagging
- 11 Redo Applycal and Inspect
- 12 Split
- 13 UV Continuum Subtraction and Setting Up for Self-Calibration
- 14 Velocity Systems and Doppler corrections
- 15 Image the Spectral Line Data
- 16 Image the Continuum data
- 17 Image Analysis and Viewing
- 18 Spectrum Fitting
- 19 Self-Calibration
This tutorial describes the data reduction for two spectral lines observed toward the AGB star IRC+10216. In this EVLA "OSRO1" mode observation one subband was observed in each of two basebands, with the subbands centered on the HC3N and SiS lines near 36 GHz. The raw data were loaded into CASA with importevla, where zero and shadowed data were flagged (note that this is done by the NRAO Data Archive). Then the data were split, so we could average from the native 1-second integrations to 10 seconds, select only antennas with Ka-band receivers, and select only spectral windows (called spw in CASA) with Ka-band data. This produces a significantly smaller dataset for processing.
IRC+10216 (CW Leo) is the brightest star in the sky at 5 microns but only 16th magnitude visually. It was discovered during the first survey of the infrared sky, carried out by Bob Leighton and Gerry Neugebauer in 1965. An Asymptotic Giant Branch star, it is a Mira-type variable going through prodigious episodic mass loss. The dust condensed from the atmosphere during the mass loss is responsible for the millimeter emission; the continuum emission seen at radio wavelengths probes the actual stellar photosphere. Molecules form along with the dust, and a steady state chemistry occurs in the dense inner regions (Tsuji 1973 A&A 23, 411). As the density of material drops, the chemistry freezes. But the molecules continue their long coast outward into the Galaxy, and as the shell thins ultraviolet light from the ambient galactic radiation field penetrates and initiates a new chemistry in the gas.
SiS, a simple molecule created in the dense inner envelope chemistry is photodissociated as it coasts out into the shell. The result is the centrally condensed emission we see in the present observations. Interestingly, in the next lower transition at 18 GHz, the line shape is much different from what we see here. At the extreme velocities in the profile, very bright narrow emission is seen which has been interpreted as maser emission. The interested student can find EVLA observations of this line in the archive.
HC3N, a much more complex species, is created by the photochemistry which becomes active as atoms and pieces of molecules destroyed by ultraviolet radiation undergo the next phase of chemistry in the shell. HC3N has many vibrational modes which may be excited in addition to its rotational modes. Owing to this, it can re-radiate energy absorbed from ultraviolet radiation more effectively than some molecules with a single bond. Eventually it too is destroyed however, but during its brief existence its rise to abundance in the envelope results in a ring of emission, which is what is observed in this image made with the EVLA. A recent model by Cordiner & Millar (2009, ApJ, 697, 68) describes a new chemical model for the shell, which also takes into account the variation of mass loss by the star. They show that in addition to purely chemical effects, local gas and dust density peaks play a role in shaping the observed emission.
Getting the data
The post-split averaged data can be downloaded from http://casa.nrao.edu/Data/EVLA/IRC10216/day2_TDEM0003_10s_norx.tar.gz (data size: 1.1GB)
Once the download is complete, unzip and unpack the file:
# in a terminal, outside of CASA: tar -xzvf day2_TDEM0003_10s_norx.tar.gz
How to Use This casaguide
There are a number of possible ways to run CASA. Many aspects are described in Getting_Started_in_CASA. You should review this page if you are new to CASA. In brief you can run CASA interactively by looking at the inputs to tasks with inp taskname (example: inp clean), setting the parameters one by one (example: selectdata=T) as you desire and then run go. After setting parameters one by one in a task and then looking at the inputs again, you will notice that the parameters that have been set to something other than their defaults are blue. If you have mistyped any parameters, they will be red and must be fixed for the task to run correctly. You can get more detailed help on any task by typing help taskname (example: help clean). Once a task is run you can get the same parameters back by running tget taskname (example: tget clean); subsequent runs will overwrite the previous tget file.
The second way to run CASA is to provide task function calls. This tutorial is made up of such calls, which were developed by looking at the inputs for each task and deciding what needed to be changed from default values. For task function calls, only parameters that you want to be different from their defaults need to be set. A series of task function calls can be combined together into a script, and run with execfile('scriptname.py'). It is possible to extract a script containing all the CASA task function calls in this and other casaguides using the method described at the Extracting_scripts_from_these_tutorials page.
If you are a relative novice or just new to CASA it is strongly recommended to work through this tutorial by cutting and pasting the task function calls provided below after you have read all the associated explanations. Work at your own pace, look at the inputs to the tasks to see what other options exist, and read the help files. Later, when you are more comfortable, you might try to extract the script, modify it for your purposes, and begin to reduce other data.
The Observing Log, Antenna Position Corrections, Opacities, Gaincurves, and other Calibration "Priors"
For all EVLA observations, the operators keep an observing log. You can look at the observing logs at the observing log website by entering a date range covering only 26 Apr 2010, then selecting "TDEM0003". Pertinent information from this observation is repeated below:
INFORMATION FROM OBSERVING LOG: Date of the observation: 26-April-2010 There are no Ka-band receivers on ea11, ea13, ea14, ea16, ea17, ea18, ea26 Antenna ea06 is out of the array Antenna 15 has a focus timeout Antennas ea10, ea12, ea22 do not have good baseline positions
As mentioned in the log, antennas ea10, ea12, and ea22 do not have good baseline positions. Antenna ea10 was not in the array, but for the other two antennas we need to check for any improved baseline positions that were derived after the observations were taken.
In CASA, these corrections can derived by gencal for the JVLA. The parameter caltype='antpos' will us antenna position offset values and create a calibration table for the delay changes that this correction involves. If the 'antenna' parameter is not specified, gencal will query the EVLA webpages for the offsets (EVLA/VLA Baseline Corrections page).
# In CASA gencal(vis='day2_TDEM0003_10s_norx',caltable='antpos.cal', caltype='antpos', antenna='')
which return the antenna offsets in the logger:
Determine antenna position offests from the baseline correction database offsets for antenna ea01 : 0.00000 0.00300 0.00000 offsets for antenna ea02 : -0.00080 0.00000 0.00000 offsets for antenna ea03 : -0.00280 0.00000 0.00000 offsets for antenna ea05 : 0.00000 0.00280 0.00000 offsets for antenna ea12 : -0.01000 0.00450 -0.00170 offsets for antenna ea22 : -0.02570 0.00270 -0.01900 offsets for antenna ea23 : -0.00140 0.00000 0.00000 offsets for antenna ea24 : -0.00150 0.00000 0.00000 offsets for antenna ea27 : 0.00000 0.00190 -0.00160
Please note: if you are reducing VLA data taken before March 1, 2010, you need to set caltype='antposvla'. For the VLA before the conversion to the EVLA or JVLA, the automated lookup will not work and one has to specify the the antenna positions explicitly.
There are two other types of "priors" that we will need to apply at each calibration stage described below:
- Gaincurve -- the gaincurve describes how each antenna behaves as a function of elevation, for each receiver band. Currently only gaincurves for the VLA/EVLA are available (see  for the incorporated models). This option should not be used with any other telescopes.
- Opacity correction -- the opacity of the observation can be computed from a seasonal model and/or weather station information. We are planning to have a task available for this information. At the moment, the CASA_EVLA_Scripts page hosts a script and a contributed task to display the weather information and to calculate the zenith opacities for each spectral window. After the zenith opacities are derived, they will be recomputed for the correct elevation of the data automatically using in gaincal, applycal, bandpass etc.
Now we want to plot the opacity of the atmosphere at 36 GHz at the time this observation was taken so it can be corrected for in subsequent calibrations. The task plotweather takes care of that. The task will plot the weather conditions at the time of the observations and calculate the atmospheric opacities basedon these data in combination with a seasonal model that contains long-term statistics at the VLA site. 'seasonal_weight=0.5' gives both equal weights:
We will be running plotweather in a way that will assign the opacity list (one entry for each spw in ascending order) to the variable myTau:
# In CASA myTau = plotweather(vis='day2_TDEM0003_10s_norx', doPlot=T)
With this input, the task
- returns the opacity values to the logger:
SPW : Frequency (GHz) : Zenith opacity (nepers) plotweather 0 : 36.387 : 0.039 plotweather 1 : 36.305 : 0.039 plotweather wrote weather figure: day2_TDEM0003_10s_norx.plotweather.png
- creates a file "day2_TDEM0003_10s_norx.plotweather.png" with the elevation of the sun, the wind speed and direction, the temperature, precipitable water vapor (PWV) as functions of time over the observation (view this file with your preferred image viewer like gthumb, xv or Preview),
- and assigns the myTau variable to the list of opacities per spectral window:
# In CASA myTau [0.039229977911770589, 0.038993893363737826]
The calibration task opacity parameters can now take either the variable directly opacity = myTau or the list of opacities like opacity = [0.0392, 0.0390] or, if a single value is good enough to correct for all spectral windows, one can simple use opacity = 0.039.
Initial Inspection and Flagging
listobs provides almost all relevant observational parameters such as correlator setup (frequencies, bandwidths, channel number and widths, polarization products), sources, scans, scan intents, and antenna locations:
# In CASA listobs(vis='day2_TDEM0003_10s_norx')
Below we have cut and pasted the most relevant output from the logger.
Fields: 4 ID Code Name RA Decl Epoch SrcId nVis 2 D J0954+1743 09:54:56.8236 +188.8.131.524 J2000 2 65326 3 NONE IRC+10216 09:47:57.3820 +184.108.40.20600 J2000 3 208242 5 F J1229+0203 12:29:06.6997 +02.03.08.5982 J2000 5 10836 7 E J1331+3030 13:31:08.2880 +220.127.116.1189 J2000 7 5814 (nVis = Total number of time/baseline visibilities per field) Spectral Windows: (2 unique spectral windows and 1 unique polarization setups) SpwID #Chans Frame Ch1(MHz) ChanWid(kHz)TotBW(kHz) Ref(MHz) Corrs 0 64 TOPO 36387.2295 125 8000 36387.2295 RR RL LR LL 1 64 TOPO 36304.542 125 8000 36304.542 RR RL LR LL Sources: 10 ID Name SpwId RestFreq(MHz) SysVel(km/s) 0 J1008+0730 0 0.03639232 -0.026 0 J1008+0730 1 0.03639232 -0.026 2 J0954+1743 0 0.03639232 -0.026 2 J0954+1743 1 0.03639232 -0.026 3 IRC+10216 0 0.03639232 -0.026 3 IRC+10216 1 0.03639232 -0.026 5 J1229+0203 0 0.03639232 -0.026 5 J1229+0203 1 0.03639232 -0.026 7 J1331+3030 0 0.03639232 -0.026 7 J1331+3030 1 0.03639232 -0.026 Antennas: 19: ID Name Station Diam. Long. Lat. 0 ea01 W09 25.0 m -18.104.22.168 +22.214.171.124 1 ea02 E02 25.0 m -107.37.04.4 +33.54.01.1 2 ea03 E09 25.0 m -126.96.36.199 +188.8.131.52 3 ea04 W01 25.0 m -107.37.05.9 +33.54.00.5 4 ea05 W08 25.0 m -184.108.40.206 +220.127.116.11 5 ea07 N06 25.0 m -107.37.06.9 +18.104.22.168 6 ea08 N01 25.0 m -107.37.06.0 +33.54.01.8 7 ea09 E06 25.0 m -22.214.171.124 +126.96.36.199 8 ea12 E08 25.0 m -188.8.131.52 +184.108.40.206 9 ea15 W06 25.0 m -220.127.116.11 +18.104.22.168 10 ea19 W04 25.0 m -22.214.171.124 +126.96.36.199 11 ea20 N05 25.0 m -107.37.06.7 +33.54.08.0 12 ea21 E01 25.0 m -107.37.05.7 +188.8.131.52 13 ea22 N04 25.0 m -107.37.06.5 +33.54.06.1 14 ea23 E07 25.0 m -184.108.40.206 +220.127.116.11 15 ea24 W05 25.0 m -18.104.22.168 +22.214.171.124 16 ea25 N02 25.0 m -107.37.06.2 +33.54.03.5 17 ea27 E03 25.0 m -107.37.02.8 +33.54.00.5 18 ea28 N08 25.0 m -107.37.07.5 +126.96.36.199
In addition to source names, antenna names, etc. CASA indexes each of these parameters starting from 0. To specify, e.g., a source field, one can use the index or the name -- both are accepted.
Note that the Rest Frequency and Systemic Velocity are wrong in the listobs log by a factor 10^6 and 1000, respectively, given the quoted units (MHz) and (km/s). This was due to a temporary error in the EVLA Observing Tool that has subsequently been fixed. Because the sky frequencies are correct, and we set the rest frequency explicitly later in the deconvolution stage, this does not present a problem for the data reduction.
Summary of Observing Strategy Gain Calibrator: J0954+1743, field id=2 Bandpass Calibrator: J1229+0203, field id=5 Flux Calibrator: J1331+3030 (3C286), field id=7 Target: IRC+10216, field id=3 Ka-band spws = 0,1
Look at a graphical plot of the antenna locations and save hardcopy in case you want it later. This will be useful for picking a reference antenna -- typically a good choice is an antenna close to the center of the array. Unless it shows problems after inspection of the data, we provisionally choose ea02.
# In CASA plotants(vis='day2_TDEM0003_10s_norx',figfile='ant_locations.png')
Next, let's look at the elevation as a function of time for all sources. It's not the case for these data, but if the elevation is very low (usually at start or end of track) you may want to flag. Also, how near in elevation your flux calibrator is to your target will impact your ultimate absolute flux calibration accuracy.
Unfortunately, the target and flux calibrator are not particularly well-matched for this observation, as you can show by plotting the elevation for each source (each sources has a different colors). We will be using data in spectral window 0, channels 4 to 60 for this plot as given by the spw = '0:4~60' parameter; the general CASA selection syntax is described in the Measurement Selection Syntax Document:
# In CASA plotms(vis='day2_TDEM0003_10s_norx', xaxis='time',yaxis='elevation',correlation='RR,LL', avgchannel='64',spw='0:4~60', coloraxis='field')
Thus we are strongly dependent on the opacity and gaincurve corrections to get the flux scale right for these data. (This is something to keep in mind when planning observations!)
Next, let's look at all the source amplitudes as a function of time:
# In CASA plotms(vis='day2_TDEM0003_10s_norx', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60', coloraxis='field')
Now zoom in on the region very near zero amplitude for sources J0954+1743 and IRC+10216. To zoom, select the Zoom tool in lower left corner of the plotms GUI, then you can left click to draw a box. Look for the low values (you may want to zoom a few times to really see the suspect points clearly). Now use the Mark Region and Locate buttons (located along the bottom of the GUI) to see which antenna is causing problems. The output is be shown in the logger. Since all the "located" baselines include ea12, this is the responsible antenna.
IMPORTANT NOTES ON PLOTMS: * When using the locate button it is important to have only selected a modest number of points with the mark region tool (see example of marked region in the thumbnail), otherwise the response will be very slow and possibly hang the tool (all of the information will be output to your terminal window, not the logger). * Throughout the tutorial, when you are done marking/locate use the Clear Regions tool to get rid of the marked box before plotting other things. * After flagdata command flagging, you have to force a complete reload of the cache to look at the same plot again with the new flags applied. To do this, either check the "force reload" box in the lower left, or do Shift+Plot. * If the plotms tool hangs during a plot, try clicking the cancel button on the load progress GUI, and/or if you see a "table locked" message try typing "clearstat" on the CASA command line. * Occasionally, plotms will get into a strange state that you cannot clear from inside. We are working on these issues, but for now, when all else fails, exit from casapy and restart.
Check the other sideband by changing spw to 1:4~60. You will have to re-zoom. If you have trouble, click on the Mark icon and then back to zoom. In spw 1, ea07 is bad from the beginning until after next pointing run: 03:21:40~04:10:00. To see this, compare the amplitudes when antenna is set to 'ea07' and when it is set to one of the other antennas, such as 'ea08'.
If you set antenna to 'ea12' and zoom in on this initial timerange, you can also see that ea12 is bad during the same time range as for spw 0. You can also see this by entering '!ea07' for antenna, which removes ea07 from the plot (in CASA selection, "!" deselects).
We can set up a flagging command to get both bad antennas for the appropriate time and spw:
# In CASA flagdata(vis='day2_TDEM0003_10s_norx', field=['2,3','2,3'], spw=['','1'], antenna=['ea12','ea07'], timerange=['03:41:00~04:10:00','03:21:40~04:10:00'])
Note: In CASA 3.4 flagdata spits out a warning and suggests to use 'tflagdata'. For the purpose of this tutorial, we continue to use the current 'flagdata'. In CASA 3.5, 'tflagdata' will replace 'flagdata' and the tutorial will switch to the new syntax.
flagdata works by spanning up a matrix. The first entries in each list must be taken as one flagging command, as well as the second entries etc. Lists within lists are fine. In the above example, the first flagging command is issued for fields 2 and 3 for all spws and within the 03:41:00~04:10:00 timerange. A second command is again for the fields 2 and 3 but for spw 1 only and for the second timerange in the list '03:21:40~04:10:00'.
Note that because the chosen timerange is limited to fields 2 and 3, the field parameter is not really needed; however, flagdata will run fastest if you put as many constraints as possible (since it will traverse less data).
Now remove the !ea07 from antenna and replot both spw, zooming in to be sure that all obviously low points are gone. Also zoom in and check 3C286 (J1229+0203 is already obvious because it is so bright!).
Let's look more closely at IRC+10216:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='3', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0~1:4~60', coloraxis='spw')
You can see a that there are some noisy high points. But now try:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='3', xaxis='uvdist',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0~1:4~60', coloraxis='spw')
Now, we can see that most of the high points on IRC+10216 are due to large scale emission on short baselines, but there is still some noisy stuff -- for a target like this with extended emission it's best to wait until later to decide what to do about it. We will not be able to get adequate calibration for antennas that are truly bad (even if they don't stand out here) so these will be obvious later.
Set Up the Model for the Flux Calibrator
Next, we set the model for the flux calibrator. Depending on your observing frequency and angular resolution you can do this several ways. In the past, one typically used a point source (constant flux) model for the flux calibrator, possibly with a uvrange cutoff if necessary. More recently for the VLA/EVLA, model images for the most common flux calibrators have been made available for use in cases where the sources are somewhat resolved. This is most likely to be true at higher frequencies and at higher resolutions (more extended arrays).
Currently, CASA contains models for the most common calibrators at the most common frequencies but not yet the new EVLA frequency bands S and Ka. These will be added as soon as they become available. One may check the availability of calibration models in setjy:
# In CASA setjy(vis='day2_TDEM0003_10s_norx',listmodels=T)
The terminal will now show the available models, e.g. 3C286_C.im, 3C48_K.im etc. (setjy will search in the working directory for images that may contain models, as well as in a CASA directory where known calibrator models are stored.)
As mentioned above, the Ka band does not yet have a model incorporated at this time as the full configuration cycle has not ended yet. However, models for Ka-band that are good for D, C, and B configurations are available from the EVLA calibrator model webpage for the new frequency bands.
From this page, please download the file 36.68 GHz CASA image of 3C286, this is the direct link. To do so, you may run
# in a terminal, outside of CASA: wget http://www.aoc.nrao.edu/evlacal/3C286_36.68G.model.tgz tar xzvf 3C286_36.68G.model.tgz
to create the model 3C286_36.68G.model on disk.
setjy scales the total flux in the model image to that appropriate for your individual spectral window frequencies according to the calibrator's flux and reports this number to the logger -- it is a good idea to save this information for your records.
# In CASA setjy(vis='day2_TDEM0003_10s_norx',field='7',spw='0~1', modimage='3C286_36.68G.model')
The logger output for each spw is: J1331+3030 (fld ind 7) spw 0 [I=1.7762, Q=0, U=0, V=0] Jy, (Perley-Butler 2010) Using model image /Users/jott/casa/data/3.4-scripts/irc/3C286_36.68G.model Scaling spw 0's model image to I = 1.77625 Jy for visibility prediction. J1331+3030 (fld ind 7) spw 1 [I=1.7794, Q=0, U=0, V=0] Jy, (Perley-Butler 2010) Scaling spw 1's model image to I = 1.7794 Jy for visibility prediction.
The absolute fluxes for the frequencies have now been determined, and one can proceed to the bandpass and complex gain calibrations.
Bandpass and Delay
Before determining the bandpass solution, we need to inspect phase and amplitude variations with time and frequency on the bandpass calibrator to decide how best to proceed. We limit the number of antennas to make the plot easier to see. We chose ea02 as it seems like a good candidate for the reference antenna.
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='5', xaxis='channel',yaxis='phase',correlation='RR', avgtime='1e8',spw='0:4~60',antenna='ea02&ea23')
If you find the points difficult to see, go to the "Display" tab, set "Unflagged Points Symbol" to "Custom", and increase the "Style" number, then replot. The phase variation is modest ~10 degrees. Now expand to all baselines that include ea02, and add an extra dimension of color:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='5', xaxis='channel',yaxis='phase',correlation='RR', avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
From this you can see that the phase variation across the bandpass is modest. Next check LL, and spw=1, both correlations. Also check other antennas if you like.
Now look at the phase as a function of time:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='5', xaxis='time',yaxis='phase',correlation='RR', avgchannel='64',spw='0:4~60',antenna='ea02&ea23')
Expand to all antennas with ea02:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='5', xaxis='time',yaxis='phase',correlation='RR', avgchannel='64',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
You can see that the phase variations are smooth, but do vary significantly over the 5 minutes of observation -- in most cases by a few 10s of degrees. Zoom in to see this better if you want.
The conclusion from this investigation is that you need to correct the phase variations with time before solving for the bandpass to prevent decorrelation of the vector averaged bandpass solution. Since the phase variation as a function of channel is modest, you can average over several channels to increase the signal to noise of the phase vs. time solution. If the phase variation as a function of channel is larger you may need to use only a few channels to prevent introducing delay-based closure errors as can happen from averaging over non-bandpass corrected channels with large phase variations.
Let's take care of the delay calibration first. The delay is the slope of phase across frequency. Form the plot we saw that it was rather modest and the bandpass calibration will certainly take care of it. Nevertheless, it might be good to derive a delay calibration first and the calculate the bandpass. The delay calibration is an antenna based calibration solution and it can be derived in gaincal with caltype='K' .
# In CASA gaincal(vis='day2_TDEM0003_10s_norx', caltable='delays.cal', field='5', refant='ea02', gaintype='K', gaintable=['antpos.cal'])
We use the strong bandpass calibrator for the solution. It will be extrapolated in time to all observations.
This call also exemplifies how the CASA calibration table system works. Calibration tables are incremental. So we supply the previous 'antpos.cal' table as input to gaincal. The following calibration step will require gaintable=['antpos.cal','delay.cal'] and so on.
Now we proceed to the actual bandpass calibration: Since the bandpass calibrator is quite strong we do the phase-only solution on the integration time of 10 seconds (solint='int').
Remember that previously we determined the myTau variable as a list of opacities. It's repeated here but may still be in your CASA session.
# In CASA myTau=[0.039229977911770589, 0.038993893363737826] gaincal(vis='day2_TDEM0003_10s_norx',caltable='bpphase.gcal', field='5',spw='0~1:20~40', refant='ea02',calmode='p',solint='int',minsnr=2.0, gaintable=['antpos.cal','delays.cal'], opacity=myTau,gaincurve=T)
Plot the solutions:
# In CASA plotcal(caltable='bpphase.gcal',xaxis='time',yaxis='phase', iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
These solutions will appear in the CASA plotter gui. If you closed it after plotting the antennas above, it should reopen. If it is still open from before, the new plots should just appear. After you are done looking at the first set of plots, push the "Next" button on the GUI to see the next set of antennas.
Now we can apply this phase solution on the fly while determining the bandpass solutions on the timescale of the bandpass calibrator scan (solint='inf').
We also use the opacity list now instead of myTau -- but both options will work.
# In CASA bandpass(vis='day2_TDEM0003_10s_norx',caltable='bandpass.bcal',field='5', refant='ea02',solint='inf',solnorm=T, gaintable=['antpos.cal','delays.cal','bpphase.gcal'], opacity=0.039,gaincurve=T)
A few words about solint and combine:
The use of solint='inf' in bandpass will derive one bandpass solution for the whole J1229+0203 scan. Note that if there had been two observations of the bandpass calibrator (for example), this command would have combined the data from both scans to form one bandpass solution, because the default of the combine parameter for bandpass is combine='scan'. To solve for one bandpass for each bandpass calibrator scan you would also need to include combine=' ' in the bandpass call. In all calibration tasks, regardless of solint, scan boundaries are only crossed when combine='scan'. Likewise, field (spw) boundaries are only crossed if combine='field' (combine='spw'), the latter two are not generally good ideas for bandpass solutions.
Plot the solutions, amplitude and phase:
# In CASA plotcal(caltable='bandpass.bcal',xaxis='chan',yaxis='amp', iteration='antenna',subplot=331)
# In CASA plotcal(caltable='bandpass.bcal',xaxis='chan',yaxis='phase', iteration='antenna',subplot=331)
Note the scale for the phasesm, they are very small since the delays were already taken out in a previous setp.
This step isn't necessary from a calibration perspective, but if you want to go ahead and check the bandpass calibration on the bandpass calibrator you can run applycal here. In the future we hope to plot corrected data on-the-fly without this applycal step. Subsequent applycals will overwrite this one, so no need to worry.
applycal(vis='day2_TDEM0003_10s_norx',field='5', gaintable=['antpos.cal','delays.cal','bandpass.bcal'], gainfield=['','5','5'], opacity=0.039,gaincurve=T,calwt=F)
Similar to flagdata, applycal works like a matrix. The first entries in the lists are to be used together, so are the second entries etc. (except for the opacity list, which is referring to spws). All will be applied to the 'field' selection. In the above example, 'antpos.cal' from any field is applied to source '5', and the 'bandpass.cal' that was obtained for field '5' (the bandpass observation) is also applied to field '5'. Again, lists within the lists are fine.
plotms(vis='day2_TDEM0003_10s_norx',field='5', xaxis='channel',yaxis='phase',ydatacolumn='corrected', correlation='RR', avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
plotms(vis='day2_TDEM0003_10s_norx',field='5', xaxis='channel',yaxis='amp',ydatacolumn='corrected', correlation='RR', avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
Note that the phase and amplitude as a function of channel are very flat now.
Now that we have a bandpass solution to apply, we can solve for the antenna-based phase and amplitude gain calibration. Since the phase changes on a much shorter timescale than the amplitude, we will solve for them separately. In particular, if the phase changes significantly over a scan time, the amplitude would be decorrelated, if the un-corrected phase were averaged over this timescale. Note that we re-solve for the gain solutions of the bandpass calibrator, so we can derive new solutions that are corrected for the bandpass shape. Since the bandpass calibrator will not be used again, this is not strictly necessary, but is useful to check its calibrated flux density for example. We use a minimum signal-to-noise of 2 here as it seems to be a good compromise for using good data without rejecting too many solutions (minsnr = 2).
# In CASA gaincal(vis='day2_TDEM0003_10s_norx',caltable='intphase.gcal', field='2,5,7',spw='0~1:4~60', refant='ea02',calmode='p',solint='int',minsnr=2.0, gaintable=['antpos.cal','delays.cal','bandpass.bcal'], opacity=0.039,gaincurve=T)
Here solint='int' coupled with calmode='p' will derive a single phase solution for each 10 second integration. Note that the bandpass table is applied on-the-fly before solving for the phase solutions, however the bandpass is NOT applied to the data permanently until applycal is run later on.
Note that quite a few solutions are rejected due to SNR < 2 (printed to terminal). For the most part it is only one or two solutions out of >30 so this isn't too worrying. Take note if you see large numbers of rejected solutions per integration. This is likely an indication that solint is too short for the S/N of the data.
Now look at the phase solution, and note the obvious scatter within a scan time:
# In CASA plotcal(caltable='intphase.gcal',xaxis='time',yaxis='phase', iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
Although solint='int' (i.e. the integration time of 10 seconds) is the best choice to apply before for solving for the amplitude solutions, it is not a good idea to use this to apply to the target. This is because the phase-scatter within a scan can dominate the interpolation between calibrator scans. Instead, we also solve for the phase on the scan time, solint='inf' (but combine=' ', since we want one solution per scan) for application to the target later on. Unlike the bandpass task, for gaincal, the default of the combine parameter is combine=' '.
# In CASA gaincal(vis='day2_TDEM0003_10s_norx',caltable='scanphase.gcal', field='2,5,7',spw='0~1:4~60', refant='ea02',calmode='p',solint='inf',minsnr=2.0, gaintable=['antpos.cal','delays.cal','bandpass.bcal'], opacity=0.039,gaincurve=T)
# In CASA plotcal(caltable='scanphase.gcal',xaxis='time',yaxis='phase', iteration='antenna',subplot=331,plotrange=[0,0,-180,180])
Note that there are no failed solutions here because of the added S/N afforded by the longer solint. Alternatively, instead of making a separate phase solution for application to the target, one can also run smoothcal to smooth the solutions derived on the integration time.
Next we apply the bandpass and solint='int' phase-only calibration solutions on-the-fly to derive amplitude solutions. Here the use of solint='inf', but combine=' ' will result in one solution per scan interval.
# In CASA gaincal(vis='day2_TDEM0003_10s_norx',caltable='amp.gcal', field='2,5,7',spw='0~1:4~60', refant='ea02',calmode='ap',solint='inf',minsnr=2.0, gaintable=['antpos.cal','delays.cal','bandpass.bcal','intphase.gcal'], opacity=0.039,gaincurve=T)
Now let's look at the resulting phase solutions. Since we have taken out the phase as best we can by applying the solint='int' phase-only solution, this plot will give a good idea of the residual phase error. If you see scatter of more than a few degrees here, you should consider going back and looking for more data to flag, particularly bad timeranges etc.
# In CASA plotcal(caltable='amp.gcal',xaxis='time',yaxis='phase', iteration='antenna',subplot=331,plotrange=[-1,-1,-20,20])
Note that we have restricted the plot range so that all antennas are plotted with the same scale (instead of autoscaled; the "-1,-1" values for the x-axis allow autoscaling to take place still). We can see that both antenna ea12 (all times) and ea23 (first 1/3 of observation) show particularly large residual phase noise.
# In CASA plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp', iteration='antenna',subplot=331,plotrange=[-1,-1,0,1])
Note that the amplitude solutions for ea12 are very low; this is another indication that this antenna is dubious.
Next we use the flux calibrator (whose flux density was set in setjy above) to derive the flux of the other calibrators. Note that the flux table REPLACES the amp.gcal in terms of future application of the calibration to the data, i.e. the flux table contains both the amp.gcal and flux scaling. Unlike the gain calibration steps, this is not an incremental table.
# In CASA fluxscale(vis='day2_TDEM0003_10s_norx',caltable='amp.gcal', fluxtable='flux.cal',reference='7')
It is a good idea to note down for your records the derived flux densities:
Flux density for J0954+1743 in SpW=0 is: 0.237109 +/- 0.0312564 (SNR = 7.58592, N = 38) Flux density for J0954+1743 in SpW=1 is: 0.240604 +/- 0.0302218 (SNR = 7.96126, N = 38) Flux density for J1229+0203 in SpW=0 is: 28.2701 +/- 0.245906 (SNR = 114.963, N = 38) Flux density for J1229+0203 in SpW=1 is: 28.3918 +/- 0.257522 (SNR = 110.25, N = 38)
Next, check that the flux.cal table looks as expected:
# In CASA plotcal(caltable='flux.cal',xaxis='time',yaxis='amp', iteration='antenna',subplot=331)
Note the odd behavior for antennas ea07, ea12, and ea23.
Applycal and Inspect
Now we apply the calibration to each source, according to which tables are appropriate, and which source should be used to do that particular calibration. For the calibrators, all bandpass solutions come from the bandpass calibrator (id=5), and the phase and amplitude calibration comes from their own solutions.
Note: In applycal we set calwt=F. It is very important to turn off this parameter which determines if the weights are calibrated along with the data. Data from antennas with better receiver performance and/or longer integration times should have higher weights, and it can be advantageous to factor this information into the calibration. During the VLA era, meaningful weights were available for each visibility. However, at the time of this observation, the EVLA was not yet recording the information necessary to calculate meaningful weights. Since these data weights are used at the imaging stage you can get strange results from having calwt=T when the input weights are themselves not meaningful, especially for self-calibration on resolved sources (your flux calibrator and target, for example).
For more recent data, the switched power information is recorded, but we currently do not recommend using this information to calculate data weights without exercising considerable caution.
# In CASA # for the gain/phase calibrator (field '2'): applycal(vis='day2_TDEM0003_10s_norx',field='2', gaintable=['antpos.cal','delays.cal','bandpass.bcal','intphase.gcal','flux.cal'], gainfield=['','5','5','2','2'], opacity=0.039,gaincurve=T,calwt=F)
# In CASA # for the bandpass calibrator (field '5'): applycal(vis='day2_TDEM0003_10s_norx',field='5', gaintable=['antpos.cal','delays.cal','bandpass.bcal','intphase.gcal','flux.cal'], gainfield=['','5','5','5','5'], opacity=0.039,gaincurve=T,calwt=F)
# In CASA # for the flux calibrator (field '7'): applycal(vis='day2_TDEM0003_10s_norx',field='7', gaintable=['antpos.cal','delays.cal','bandpass.bcal','intphase.gcal','flux.cal'], gainfield=['','5','5','7','7'], opacity=0.039,gaincurve=T,calwt=F)
For the target we apply the bandpass from id=5, and the calibration from the gain calibrator (id=2):
# In CASA # for the target source IRC10216 (field '3'): applycal(vis='day2_TDEM0003_10s_norx',field='3', gaintable=['antpos.cal','delays.cal','bandpass.bcal','scanphase.gcal','flux.cal'], gainfield=['','5','5','2','2'], opacity=0.039,gaincurve=T,calwt=F)
Now inspect the corrected data:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='5',ydatacolumn='corrected', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna1')
This plot shows some data deviating from the average amplitudes. Use methods described above to mark a region for a small number of deviant data points, and click "Locate". You will find that ea12 is responsible.
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='2',ydatacolumn='corrected', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna2')
Here we see some problems, with high points. Mark some regions and locate in plotms to find out which antennas and in which spws. Pay special attention to antennas that have been called out already as showing some dubious behavior.
What you find is that ea07 which we flagged spw=1 above, is also bad for the same timerange in spw=0. This was not obvious in the raw data, because spw=0 was adjusted in the on-line system by a gain attenuator, while spw=1 wasn't. So a lack of power on this antenna can look like very low (and obvious) amplitudes in spw=1 but not for spw=0. Looking carefully you'll see that ea07 is actually pretty noisy throughout.
From the locate we also find that ea12 and ea23 show some high points; to see this, replot baselines using each of them alone:
plotms(vis='day2_TDEM0003_10s_norx',field='2',ydatacolumn='corrected', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60',antenna='ea12', coloraxis='antenna2')
plotms(vis='day2_TDEM0003_10s_norx',field='2',ydatacolumn='corrected', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60',antenna='ea23', coloraxis='antenna2')
It may be a a good idea to flag ea12 completely - it's just a bit noisy all around and ea23 is pretty noisy during the first scans between initial and second pointing. Recall that these are antennas we became suspicious of while inspecting the calibration solutions.
Now let's see how the target looks. Because the target has resolved structure, its best to look at it as a function of uvdistance. We'll go ahead and exclude the three antennas we already know have problems.
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected', xaxis='uvdist',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60',antenna='!ea07;!ea12;!ea23', coloraxis='antenna2')
The color indicates that the spikes are caused by a single antenna. Use, zoom, mark, and locate to see which one. Also look at spw=1.
The bad antenna turns out to be ea28; to confirm, replot with antenna=!ea28:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected', xaxis='uvdist',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60',antenna='!ea07;!ea12;!ea23;!ea28', coloraxis='antenna2')
Looks much better! To see if it's restricted to a certain time, do:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60',antenna='ea28', coloraxis='antenna1')
Baselines with ea28 clearly show issues until about two-thirds of the way through the observation. Plot another distant antenna to compare. We will go ahead and flag it all, since it's hanging far out on the north arm by itself.
The additional data we've identified as bad need to be flagged, and then all the calibration steps will need to be run again.
# In CASA flagdata(vis='day2_TDEM0003_10s_norx', field=['',''], spw=['',''], antenna=['ea07,ea12,ea28','ea07,ea23'], timerange=['','03:21:40~04:10:00'])
Redo Calibration after more Flagging
After flagging, you'll need to repeat the calibration steps above. Here, we append _redo to the table names to distinguish them from the first round, in case we want to compare with previous versions.
# In CASA gaincal(vis='day2_TDEM0003_10s_norx',caltable='bpphase_redo.gcal', field='5',spw='0~1:20~40', refant='ea02',calmode='p',solint='int',minsnr=2.0, gaintable=['antpos.cal','delays.cal'], opacity=0.039,gaincurve=T) # bandpass(vis='day2_TDEM0003_10s_norx',caltable='bandpass_redo.bcal', field='5', refant='ea02',solint='inf',solnorm=T, gaintable=['antpos.cal','delays.cal','bpphase_redo.gcal'], opacity=0.039,gaincurve=T) # gaincal(vis='day2_TDEM0003_10s_norx',caltable='intphase_redo.gcal', field='2,5,7',spw='0~1:4~60', refant='ea02',calmode='p',solint='int',minsnr=2.0, gaintable=['antpos.cal','delays.cal','bandpass_redo.bcal'], opacity=0.039,gaincurve=T) # gaincal(vis='day2_TDEM0003_10s_norx',caltable='scanphase_redo.gcal', field='2,5,7',spw='0~1:4~60', refant='ea02',calmode='p',solint='inf',minsnr=2.0, gaintable=['antpos.cal','delays.cal','bandpass_redo.bcal'], opacity=0.039,gaincurve=T) # gaincal(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal', field='2,5,7',spw='0~1:4~60', refant='ea02',calmode='ap',solint='inf',minsnr=2.0, gaintable=['antpos.cal','delays.cal','bandpass_redo.bcal','intphase_redo.gcal'], opacity=0.039,gaincurve=T) # fluxscale(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal', fluxtable='flux_redo.cal',reference='7')
Flux density for J0954+1743 in SpW=0 is: 0.265024 +/- 0.0286248 (SNR = 9.25854, N = 32) Flux density for J0954+1743 in SpW=1 is: 0.264985 +/- 0.0292237 (SNR = 9.06747, N = 32) Flux density for J1229+0203 in SpW=0 is: 28.448 +/- 0.286479 (SNR = 99.3023, N = 32) Flux density for J1229+0203 in SpW=1 is: 28.3255 +/- 0.294061 (SNR = 96.3252, N = 32)
Feel free to pause here and remake the calibration solution plots from above, just be sure to put in the revised table names.
Redo Applycal and Inspect
Now, apply all the new calibrations, which will overwrite the old ones. These commands are identical to those above, with the exception of the _redo part of each calibration filename.
# In CASA applycal(vis='day2_TDEM0003_10s_norx',field='2', gaintable=['antpos.cal','delays.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'], gainfield=['','5','5','2','2'], opacity=0.039,gaincurve=T,calwt=F) # applycal(vis='day2_TDEM0003_10s_norx',field='5', gaintable=['antpos.cal','delays.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'], gainfield=['','5','5','5','5'], opacity=0.039,gaincurve=T,calwt=F) # applycal(vis='day2_TDEM0003_10s_norx',field='7', gaintable=['antpos.cal','delays.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'], gainfield=['','5','5','7','7'], opacity=0.039,gaincurve=T,calwt=F) # applycal(vis='day2_TDEM0003_10s_norx',field='3', gaintable=['antpos.cal','delays.cal','bandpass_redo.bcal','scanphase_redo.gcal','flux_redo.cal'], gainfield=['','5','5','2','2'], opacity=0.039,gaincurve=T,calwt=F)
Now you can inspect the calibrated data again. Except for random scatter things look pretty good.
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='2',ydatacolumn='corrected', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna2')
You can use the Mark and Locate buttons to assess that the remaining scatter seems random, i.e. no particular antenna or time range appears to be responsible.
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected', xaxis='uvdist',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='0~1:4~60',antenna='', coloraxis='spw')
Now we split the data into individual files. This is not strictly necessary, as you can select the appropriate fields in later clean stages, but it is safer in case for example you get confused with later processing and want to fall back to this point (this is especially a good idea if you plan to do continuum subtraction or self calibration later on). It also makes smaller individual files in case you want to copy to another machine or colleague.
Here, we split off the data for the phase calibrator and the target:
# In CASA split(vis='day2_TDEM0003_10s_norx',outputvis='J0954', field='2') # split(vis='day2_TDEM0003_10s_norx',outputvis='IRC10216', field='3')
To reinitialize the scratch columns for use by later tasks, we need to run clearcal for both new datasets
# In CASA clearcal(vis='J0954') # clearcal(vis='IRC10216')
This concludes the calibration phase of the data reductions. The tutorial continues with continuum subtraction, imaging, and image analysis in
UV Continuum Subtraction and Setting Up for Self-Calibration
Now we can make a vector averaged uv-plot of the calibrated target spectral line data. It is important to note that you will only see signal in such a plot if (1) the data are well calibrated, and (2) there is significant signal near the phase center of the observations, or if the line emission (or absorption) is weak but extended. If this isn't true for your data, you won't be able to see the line signal in such a plot and will need to make an initial (dirty or lightly cleaned) line+continuum cube to determine the line-free channels. Generally, this is the recommended course for finding the line-free channels more precisely than is being done here due to time constraints, as weak line signal would not be obvious in this plot.
plotms(vis='IRC10216',field='',ydatacolumn='corrected', xaxis='channel',yaxis='amp',correlation='RR', avgtime='1e8',avgscan=T,spw='0~1:4~60',antenna='', coloraxis='spw')
In the Display tab, change the Unflagged Points Symbol to Custom and Style of 3.
You should see the "horned profile" typical of a rotation shell. From this plot, you can guess that strong line emission is restricted to channels 18 to 47 (zoom in if necessary to see exactly what the channel numbers are).
In the Data tab, under Averaging, you can also click on "All Baselines" to average all baselines, but this is a little harder to see.
Now we want to use the line free channels to create a model of the continuum emission that can be subtracted to form a line-only dataset. We want to refrain from going to close to the edges of the band -- these channels are typically noisy, and we don't want to get too close to the line channels because we could only see strong line emission in the vector averaged uv-plot.
The "want_cont=T" will produce two new datasets, IRC10216.contsub is the continuum subtracted line data, and IRC10216.cont is the continuum estimate (note however, that it is still a multi-channel dataset).
Velocity Systems and Doppler corrections
The current incarnation of the EVLA does not support Doppler tracking. Doppler setting is possible which will calculate the sky frequency based on a velocity of the source at the start of an observation. The sky frequency is then fixed throughout that track. Typically, a fixed frequency is better for the calibration of interferometric data. The downside, however, is that ta spectral line may shift over one or more channels during an observation. clean takes care of such a shift when regridding the visibilities in velocity space (default is LSRK) to form an image. Sometimes, in particular when adding together different observing tracks, it may be advisable to do the regrid all data sets to the same velocity grid, combine all data to a single file, then Fourier transform and deconvolve. The tasks cvel, concat, and clean serve this purpose respectively. The following run of cvel shows an example on how the parameters of cvel may be set.
The IRC10216.contsub visibility spans the following channel range (see also the listobs output in the first part of the tutorial):
# In CASA vishead(vis='IRC10216.contsub', mode='summary')
SpwID #Chans Frame Ch1(MHz) ChanWid(kHz) TotBW(kHz) Corrs 0 64 TOPO 36387.2295 125 8000 RR RL LR LL 1 64 TOPO 36304.542 125 8000 RR RL LR LL
For spw 0, this corresponds to about 1 km/s channel width. If we want to image the HC3N spectral line with a rest frequency of 36.39232 GHz over a velocity range of -50km/s to 0km/s and a channel width of 5 km/s, we may decide to regrid the visibilities in cvel as
Note that this step is not necessary for the processing further down in this tutorial. You may skip it if you wish.
# In CASA cvel(vis='IRC10216.contsub', outputvis='IRC10216.contsub-cveled', mode='velocity', interpolation='linear', nchan=10, start='-50km/s', width='5km/s', restfreq='36.39232GHz',outframe='LSRK', veltype='optical')
This will create a new dataset where the data is binned into the new grid. Since all data in measurement sets are stored in frequency space, an inspection with vishead now gives:
# In CASA vishead(vis='IRC10216.contsub-cveled', mode='summary')
SpwID #Chans Frame Ch1(MHz) ChanWid(kHz) TotBW(kHz) Corrs 0 10 LSRK 36392.927 606.97375 6070.64874 RR RL LR LL
Image the Spectral Line Data
Here we make images from the continuum-subtracted, calibrated spectral line data. Because the spectral line emission from IRC+10216 has significant extended emission, it is very important to run clean interactively, and make a clean mask. To make the cube a bit smaller and stay away from noisy edge channels we restrict the channel range using the spw parameter.
Note that interrupting clean by Ctrl+C may corrupt your visibilities -- you may be better off choosing to let clean finish. We are currently implementing a command that will nicely exit to prevent this from happening, but for the moment try to avoid Ctrl+C.
# In CASA clean(vis='IRC10216.contsub',imagename='IRC10216_HC3N.cube_r0.5', imagermode='csclean', imsize=300,cell=['0.4arcsec'],spw='0:5~58', mode='velocity',interpolation='linear', restfreq='36.39232GHz',outframe='LSRK', weighting='briggs',robust=0.5, interactive=T, threshold='3.0mJy',niter=100000)
- imagermode = csclean will invoke the Cotton-Schwab cleaning algorithm and the data will be regridded into a new output velocity frame, correcting for Doppler shifts of the line during the run (EVLA data for each track is always topocentric at a fixed sky frequency). The iterations are chosen as a high value to allow many clean cycles when needed. Typically, however, the threshold will kick in earlier and stop the cleaning process.
It will take a little while to grid the data, but the viewer will open when it's ready to start an interactive clean. Use the "tape deck" at the bottom of the Viewer display GUI to step through to the channel with the most extended (in angular size) emission, select "all channels" for the clean mask, select the polygon tool (the 'R' with the wobbly line around it) and make a single mask that applies to all channels (see example in thumbnail). Once you make the polygon region, you need to double click inside it to save the mask region -- if you see the polygon turn white you will know you succeeded. Note, that if you had the time and patience you could make a clean mask for each channel, and this would create a slightly better result.
After making the mask you should check that the emission in all the other channels fits within the mask you made using the "tape deck" to move back and forth. If you need to include more area in the mask, you can chose the "erase" toggle at the top, and then encircle your existing mask with a polygon and double click inside. Then go back to "add" toggle at top and make a new mask. Alternatively, you can erase a part of the mask, or you can add to the existing mask by drawing new polygons. Feel free to experiment with this a bit.
Note: If you start an interactive clean, and then do not make a mask, clean will stop when you tell it to go on because it has nothing to clean. There is no default mask.
To continue with clean use the "Next action" buttons in the green area on the Viewer Display GUI: The red X will stop clean where you are, the blue arrow will stop the interactive part of clean, but continue to clean non-interactively until reaching the stopping niter (note that this is "iterations" x "cycles") or threshold (whichever comes first), and the green arrow will clean until it reaches the "iterations" parameter on the left side of the green area. When the interactive viewer comes back use the tape deck to recheck that your mask encompasses what you think is real emission. The middle mouse button by default controls the image stretch.
Note that for this example, threshold has been set to threshold = '3mJy' to protect you from cleaning too deeply. With a careful clean mask you can clean to close to the thermal noise limit (note here I mean the actual observed rms noise limit and not the theoretical one you calculated for the proposal, as flagging, weather, etc. can affect what you actually get). It is ALWAYS best to clean each channel in a cube to a specific threshold than to stop by simply using the niter parameter, which can leave each channel cleaned to different levels. There are many ways to determine a suitable threshold. One way is to make a dirty image (niter = 0), open the cube using the viewer, go to a line free channel, select the box region tool, make a box near the field center about the size of your source, and double click inside. The rms noise of that channel will appear in the terminal window from which the viewer was launched. Try a few different boxes, average the results and this is a good estimate of the rms per channel assuming your data are not dynamic range limited (i.e. noise can be higher in channels with strong signal). This is the absolute minimum for threshold. With no mask you probably shouldn't clean deeper than 3x this rms.
Keep cleaning, by using the green Next Action arrow until the residual displayed in the viewer looks "noise like". To speed things up, you might change the iteration parameter in the viewer to something like 300. This parameter can also be set in the task command. You will notice that in this particular case, there are residuals that cannot be cleaned -- these are due to the extended resolved out structure on size scales larger than the array is sensitive to (the "Largest Angular Scale" or LAS that the array is sensitive to can be calculated from the shortest baseline length), and potential residual phase and amplitude calibration errors. We will explore this in a few sections with self-calibration.
Repeat the process for the SiS line using the call below, note that the emission for this line is less extended than the HC3N -- this has to do with the different excitation requirements of the two different lines. The SiS is excited closer to the central star than the HC3N.
# In CASA clean(vis='IRC10216.contsub',imagename='IRC10216_SiS.cube_r0.5', imagermode='csclean', imsize=300,cell=['0.4arcsec'],spw='1:5~58', mode='velocity',interpolation='linear', restfreq='36.30963GHz',outframe='LSRK', weighting='briggs',robust=0.5, interactive=T, usescratch=True, threshold='3.0mJy',niter=100000)
'usescratch=True' prepares the data for selfcal below.
You can look at both cubes using the viewer, and the tape deck to play the cube as a "movie".
# In CASA viewer
Image the Continuum data
Below the use of mode='mfs' will make a single multi-frequency synthesis image out of the specified spw/channels. Again you should make an interactive clean mask. Since no threshold is set, you will need to stop cleaning when the residual looks noise like using the red x "Next Action" button (it will be done when the viewer comes back the second time). The continuum for IRC10216 is very weak but interesting -- it is essentially tracing the photosphere of the AGB star.
The continuum data set produced with wantcont=True in uvcontsub2 is the model fit. To image the continuum itself, use the line-free channels.
# In CASA clean(vis='IRC10216',imagename='IRC10216.36GHzcont', mode='mfs',imagermode='csclean', imsize=300,cell=['0.4arcsec'],spw='0~1:5~14,0~1:48~59', weighting='briggs',robust=0.5, interactive=T)
Now look at the result in the viewer, if you like:
# In CASA viewer
Image Analysis and Viewing
Next make integrated intensity maps (moment 0) and intensity-weighted velocity maps (moment 1). For HC3N, we also produce a velocity dispersion, peak flux, and median map. All are derived with immoments. To do this, we'll want to know what channels the line emission starts and ends on, and also the rms noise in a single channel. So first lets open the viewer:
# In CASA viewer
Then use the Viewer tape deck to see which channels have significant line emission. For HC3N, the line channel range in the cube is 16 to 45, and it is the same for SiS.
Then use the tape deck to go to a line free channel, select the box region tool and make a box. When you double click in the box, the image statistics for the channel you are on will print to the terminal. Move the box around a bit to see what the variation in rms noise is. You should get something like 2 mJy. Note that the rms is much worse in channels with strong emission because of the low dynamic range of these data. If you want the box tool to go away (i.e. if you want to make a new one), hit the escape key.
Now let's make the moment 0 and moment 1 maps. For moment zero, it's best to limit the calculation to image channels with significant signal in them, but not to apply a flux cutoff, as this will bias the derived integrated intensities upward.
# In CASA immoments(imagename='IRC10216_HC3N.cube_r0.5.image',moments=, axis='spectral', chans='16~45', outfile='IRC10216_HC3N.cube_r0.5.image.mom0')
# In CASA immoments(imagename='IRC10216_SiS.cube_r0.5.image',moments=, axis='spectral', chans='16~45', outfile='IRC10216_SiS.cube_r0.5.image.mom0')
To have a look at these, use the viewer:
# In CASA viewer('IRC10216_HC3N.cube_r0.5.image.mom0') # viewer('IRC10216_SiS.cube_r0.5.image.mom0')
For moment 1, it is essential to apply a conservative flux cutoff to limit the calculation to high signal-to-noise areas. Here we use about 5σ:
# In CASA immoments(imagename='IRC10216_HC3N.cube_r0.5.image',moments=, axis='spectral', chans='16~46',excludepix=[-100,0.01], outfile='IRC10216_HC3N.cube_r0.5.image.mom1')
# In CASA immoments(imagename='IRC10216_SiS.cube_r0.5.image',moments=, axis='spectral', chans='16~45',excludepix=[-100,0.01], outfile='IRC10216_SiS.cube_r0.5.image.mom1')
Finally, we will do velocity dispersion, peak flux, and median map in a single step for HC3N. immoments can perform these steps even though the latter two are not 'moments' in a mathematical sense. Check the help file to find out the options. Peak flux and median are produced with the moment parameter set to 8 and 3:
# In CASA immoments(imagename='IRC10216_HC3N.cube_r0.5.image',moments=[2,8,3], axis='spectral', chans='16~46',excludepix=[-100,0.01], outfile='IRC10216_HC3N.cube_r0.5.image.extramoms')
will create the files IRC10216_HC3N.cube_r0.5.image.extramoms.weighted_dispersion_coord, IRC10216_HC3N.cube_r0.5.image.extramoms.median, and IRC10216_HC3N.cube_r0.5.image.extramoms.maximum.
Now use the viewer to further explore the images you've made.
For fun you can download the VLT V-band image at http://casa.nrao.edu/Data/EVLA/IRC10216/irc_fors1_dec_header.fits kindly provided by Izan Leão and overlay the moment images and 36 GHz continuum. More information about the dust properties can be found in the Leão et al. (2006) paper http://adsabs.harvard.edu/abs/2006A%26A...455..187L.
The creation of position velocity cuts from the viewer is currently being developed and hopefully available soon in CASA. If you are interested in a work-around, you may have a look at the pV casaguide. Masking the data cube to extract the emission is described here.
Frequently, one would like to fit Gaussians or polynomials to the spectral line in the data cube. This can be done with CASA's specfit task. specfit can fit those functions to an average spectrum define by some bounding box, or, alternatively, for each pixel. In the following, we will do both.
Fitting an average spectrum
First, we want to inspect the spectrum. Load the image into the viewer (here: the HC3N image cube), select "spectral profile" from the Tools menu and open a region with the mouse button that is assigned to the rectangular "R" region in the tool bar. Best to do this at a plane that shows the entire extent of the source. The average spectrum will be displayed in a separate panel.
To fit this profile in specfit, we need a region file outlining the 2-D region that is averaged (the green box in the viewer screenshot). In the following we use the new CASA region format (CASA 3.3 and higher) that is described here. Following the guidelines on that page, we create a file named specfit.crtf that describes a box with its [[x1,y1],[x2,y2]] corners in J2000 RA DEC coordinates.
#CRTFv0 box[[09:47:59.2, 13.16.24], [09:47:55.8, 13.17.09]]
We will fit 2 Gaussians to the two peaks of the spectrum (the scientific merit is debatable). A file with initial values for the fit can be provided via the estimates parameter - see specfit for details. Here we will let CASA figure out the start values by itself:
# In CASA myfit = specfit(imagename='IRC10216_HC3N.cube_r0.5.image', region='specfit.crtf', multifit=F, estimates='', ngauss=2)
Note that the output is stored in a Python dictionary called "myfit", as well as printed to the CASA logger. You should get something similar to this (depending on the details of flagging etc.):
Fit : RA : 09:47:57.49 Dec : 188.8.131.52 Stokes : I Pixel : [146.002, 164.499, 0.000, *] Attempted : YES Converged : YES Iterations : 27 Results for component 0: Type : GAUSSIAN Peak : 6.13 +/- 0.49 mJy/beam Center : -16.34 +/- 0.37 km/s 40.41 +/- 0.36 pixel FWHM : 8.81 +/- 0.89 km/s 8.56 +/- 0.87 pixel Integral : 57.5 +/- 7.4 mJy/beam.km/s Results for component 1: Type : GAUSSIAN Peak : 5.40 +/- 0.40 mJy/beam Center : -34.35 +/- 0.51 km/s 22.92 +/- 0.49 pixel FWHM : 13.4 +/- 1.3 km/s 13.0 +/- 1.3 pixel Integral : 77.2 +/- 9.4 mJy/beam.km/s
which seems to have caught the two peaks pretty well.
Spectral Fitting pixel by pixel
specfit can also fit Gaussians to every spectrum in each single spatial pixel. The following command will do this within the spectfit.crtf region defined above:
# In CASA pixfit = specfit(imagename='IRC10216_HC3N.cube_r0.5.image', region='specfit.crtf', ngauss=2, multifit=T, amp='fit.amp.image', center='fitcenter.image', fwhm='fitfwhm.image')
In this example, specfit will produce three images per Gaussian, images that map the best fit values of the Gauss peaks (amplitudes), velocity centers, and full widths at half maximum. For the first Gaussian, the image to the right displays the amplitude image fit.amp.image_0.
The many different aspects of self-calibration could fill several casaguides. Here we describe a simple process for this particular relatively low S/N data (low S/N per channel, at least).
While running clean above, the model column for each channel will have been filled with the clean model (if you made a Fourier transform of this model, you would see an image of the clean components).
We choose to do the self cal on the spw=1 SiS line data because it has the strongest emission in a single channel and is a bit more compact than the HC3N data. We will run gaincal specifying the channel in the uv-data that has the brightest peak in the image (use the viewer to figure out which channel this is for spw=1), note down what the peak flux is. Since we started the image with a channel range we need to account for the fact that the image channel numbers do not map exactly to the uv-data channel numbers (they are off by 5 so that channel 13 in the image is roughly channel 19 in the uv-data).
The next thing we need to understand is the S/N of the data. In particular, to self-cal, you need enough signal on a single baseline over the course of your chosen solint to get a S/N of about 3. Above we calculated an average rms noise of about 2 mJy/beam/channel for the whole timerange (about 95 minutes on source time) and all antennas (16). We can use our knowledge of the radiometer equation (see EVLA Sensitivity) where rms scales as 1/sqrt(time * #baselines), and the number of baselines= N(N-1)/2 and N=# of antennas. So the rms noise on one baseline, for one 10 second integration in this observation is given by:
The 95 minutes of on-source time can be estimated from a plot like this where you can sum up the amount of time on a source:
# In CASA plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected', xaxis='time',yaxis='amp',correlation='RR,LL', avgchannel='64',spw='1:4~60',antenna='')
This analysis suggests that the rms noise on one baseline, for one 10 second integration is only about 500 mJy. In contrast, the peak flux density in the strongest SiS channel is only about 200 mJy (you can check using the viewer). Since the emission is fairly compact, most baselines will see about this peak flux; this is why we choose the more compact of the two possible lines. Thus, a 10 second solution interval is not enough to get a SNR of at least 3 on a 200 mJy peak. We need to use a solint large enough so that the rms noise is not worse than about 1/3 of 200 mJy. Thus, a solint of 10 minutes is about the shortest we can use and be reasonably confident of the solutions.
Now we run gaincal with the solint we have determined. Note that because our desired solint is more than the scan time, we need to include combine='scan'.
# In CASA gaincal(vis='IRC10216.contsub',caltable='pcal_ch19one_10min', spw='1:19',calmode='p',solint='10min',combine='scan', refant='ea02',minsnr=3.0)
Now let's look at the solutions:
# In CASA plotcal(caltable='pcal_ch19one_10min',xaxis='time',yaxis='phase', iteration='antenna',subplot=331,plotrange=[0,0,-50,50])
For some antennas you can see clear global trends away from zero: ea08, ea21, and ea24 are examples, and you can also see some smaller variations with time.
Now let's explore whether applying this solution actually improves matters. To do this we need to run applycal to apply the solutions to the line dataset, both spw. We need to use spwmap to tell it that the solutions derived for spw=1 should be applied to both spw=0 and spw=1. Again it's important to set calwt=F here.
# In CASA applycal(vis='IRC10216.contsub',field='',spw='0,1', gaintable=['pcal_ch19one_10min'],spwmap=[[1,1]],calwt=F)
Note: in this example we ran the self-cal steps on the full uv continuum subtracted spectral line data set. For a more complex iterative self-calibration proceedure, you may find it easier to split off the channel/spw you want to experiment on with split, and then do all the imaging (clean) and gaincal steps with it. The gaincal tables created on the single channel can still be applied with applycal to the multi-channel/spw dataset. If you do this though, keep in mind that once split, the single-channel data will have its spw id reset to 0 (you can check with listobs), no matter what spw it came from. Thus in order to applycal with it you would need spwmap=[[0,0]].
To save time we can use the clean mask we made before and run in a non-interactive mode. You can use a mask over again as long as the number of channels in the clean call haven't changed. You can change cell or imsize and it will still do the right thing.
# In CASA clean(vis='IRC10216.contsub',imagename='IRC10216_HC3N.cube_r0.5.pselfcal', imagermode='csclean', imsize=300,cell=['0.4arcsec'],spw='0:5~58', mode='velocity',interpolation='linear', restfreq='36.39232GHz',outframe='LSRK', weighting='briggs',robust=0.5, mask='IRC10216_HC3N.cube_r0.5.mask', interactive=F,threshold='3.0mJy',niter=100000)
# In CASA clean(vis='IRC10216.contsub',imagename='IRC10216_SiS.cube_r0.5.pselfcal', imagermode='csclean',usescratch=T, imsize=300,cell=['0.4arcsec'],spw='1:5~58', mode='velocity',interpolation='linear', restfreq='36.30963GHz',outframe='LSRK', weighting='briggs',robust=0.5, mask='IRC10216_SiS.cube_r0.5.mask', interactive=F,threshold='3.0mJy',niter=100000)
Now investigate the original and self-cal'ed images in the viewer. You will find that even this single self-cal step significantly improves the images. Try opening both versions of the SiS image cubes. Then select a bright channel from the tape deck like channel 37, then use the "wrench" and "pwrench" guis to make a plot like below setting the same image range for both cubes, and two panels in x, then to see both images of that channel side-by-side click the blink toggle (see image below for more tips on setup.)
Repeat for HC3N:
Now you can redo the moment images if you like with the improved cubes (be sure to change the output file). names.
Last checked on CASA Version 3.4.0.