EVLA high frequency Spectral Line tutorial - IRC+10216 part2: Difference between revisions

From CASA Guides
Jump to navigationJump to search
 
(46 intermediate revisions by 3 users not shown)
Line 1: Line 1:
This tutorial picks up where [[EVLA Spectral Line Calibration IRC+10216]] leaves off.
This is the second part of the EVLA IRC+10216 data reduction tutorial. [http://casaguides.nrao.edu/index.php?title=EVLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216_part1 IRC+10216 Part 1: editing, calibration]
 
== UV Continuum Subtraction and Setting Up for Self-Calibration==
== UV Continuum Subtraction and Setting Up for Self-Calibration==


[[Image:irc10216_uvspec.png|thumb|UV-plot of the spectral line signal in both spw for IRC+10216.]]
[[Image:irc10216_uvspec.png|thumb|UV-plot of the spectral line signal in both spw for IRC+10216.]]
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.  
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.  


Line 9: Line 9:
plotms(vis='IRC10216',field='',ydatacolumn='corrected',
plotms(vis='IRC10216',field='',ydatacolumn='corrected',
       xaxis='channel',yaxis='amp',correlation='RR',
       xaxis='channel',yaxis='amp',correlation='RR',
       avgtime='1e8',avgscan=T,spw='0~1:4~60',antenna='', colorize=True, coloraxis='spw')
       avgtime='1e8',avgscan=T,spw='0~1:4~60',antenna='', coloraxis='spw')
</source>
</source>


in the Display tab, change the Unflagged points symbol to custom and Style of 3.  
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  
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).  
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 you can also click on "all baselines" to average all baselines, but this is a little harder to see.
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.
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.
Line 27: Line 27:


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).
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).
<!--
'''If you want to try self-cal:''' Unfortunately, at the moment, uvcontsub2 doesn't leave the continuum subtracted line dataset in the state you need *if* you think you might want to self-calibrate the data later. This is because {{clean}} always looks at the "corrected" datacolumn, while {{gaincal}} (also used for self-calibration) always looks at the "data" column. You also need to know that unless the imagename supplied to clean already exists, clean always overwrites the model column with the clean model (this will be the model supplied to the self-calibration process).
<pre style="background-color: #E0FFFF;">
The files produced by uvcontsub2 will have the following in their data columns:
                      Data      Model    Corrected
IRC10216.contsub    Line+Cont    Cont        Line
IRC10216.cont          Cont      Cont        Cont
</pre>
Now the imaging task {{clean}} will clean both of these files fine, and correctly overwrite the model data column with the correct clean model. However, if you try to self-cal (i.e. run gaincal) on the continuum subtracted line data (.contsub), it will use the Line+Cont as its input along with the line only clean model. If the line and continuum have significantly different morphology (almost always) the self-cal process will fail. In contrast the continuum dataset (.cont) will work fine because the Data and Corrected columns agree.
To fix things up, we must {{split}} the "corrected" column (which places "Corrected" in the "Data" column of the new dataset. Then you can either run clearcal to reinitialize the "model" and "corrected" columns or let clean do it for you. We put this step here explicitly for clarity about the process.
<source lang="python">
split(vis='IRC10216.contsub',outputvis='IRC10216.contsub.data')
</source>
<source lang="python">
clearcal(vis='IRC10216.contsub.data')
</source>
Now IRC10216.contsub.data will have the Line in the "data" column, 1 in the "model" column for amplitudes and 0 for phase, and the Line in the "corrected" column as desired. You can always check this using {{plotms}} to
look at the "data", "model", and "corrected" columns.
We expect that these extra shenanigans will be unnecessary in the future.
-->


==Velocity Systems and Doppler corrections==
==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}} take care of such a shift when regridding in velocity space (default is LSRK). Sometimes, in particular when adding together different tracks, it may be advisable to do the regrid all data sets to the same velocity grid, combine all data is 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 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):  
The '''IRC10216.contsub''' visibility spans the following channel range (see also the listobs output in the first part of the tutorial):  
Line 71: Line 45:
</pre>
</pre>


For spw 0, this corresponds to about 1km/s channel width. If we want to image the HC3N spectral line with a rest frequency of 36.39232GHz over a velocity range of -20km/s to 60km/s and a channel width of 4km/s, we may decide to regrid the visibilities in {{cvel}} as
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.'''


<source lang="python">
<source lang="python">
# In CASA
# In CASA
cvel(vis='IRC10216.contsub', outputvis='IRC10216.contsub-cveled', mode='velocity',
cvel(vis='IRC10216.contsub', outputvis='IRC10216.contsub-cveled', mode='velocity',
     interpolation='linear', nchan=20, start='-20km/s', width='4km/s',
     interpolation='linear', nchan=10, start='-50km/s', width='5km/s',
     restfreq='36.39232GHz',outframe='LSRK', veltype='optical')  
     restfreq='36.39232GHz',outframe='LSRK', veltype='optical')  
</source>
</source>


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}} results in
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:
 


<source lang="python">
<source lang="python">
Line 90: Line 65:
<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)  TotBW(kHz)  Corrs
SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)  TotBW(kHz)  Corrs
0          20 LSRK  36385.5234  485.382262    9710.10637 RR  RL  LR  LL  
0          10 LSRK  36392.927  606.97375    6070.64874 RR  RL  LR  LL
</pre>
</pre>


After the {{cvel}} step, the data can then be combined with other observations via {{concat}} and imaged in {{clean}} with mode='channel' to conserve that velocity system and grid.
After the {{cvel}} step, the data can then be combined with other observations via {{concat}} and imaged in {{clean}} with mode='channel' to conserve that velocity system and grid.


{{cvel}} can also Hanning smooth the data in the process if needed (as an alternative to the stand-alone {{hanningsmooth}} task).
Note that {{cvel}} can also Hanning smooth the data, if needed (as an alternative to the stand-alone {{hanningsmooth}} task).


==Image the Spectral Line Data==
==Image the Spectral Line Data==
Line 102: Line 76:
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  
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.
channel range using the spw parameter.
[[Image:viewer_interactive.png|thumb|Channel 28 shown for the HC3N cube shown in the interactive viewer with the white contour showing the mask contour drawn with the polygon tool.]]
[[Image:viewer_interactive.png|thumb|Channel 32 shown for the HC3N cube shown in the interactive viewer with the white contour showing the mask contour drawn with the polygon tool.]]
 


'''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.'''


<source lang="python">
<source lang="python">
Line 118: Line 92:
</source>
</source>


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 (second in from the right) 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.  
* 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.
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.   
'''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 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  
To continue with {{clean}} use the "Next action" buttons in the green area on the Viewer Display GUI: The red X [[File:clean-stop.png]] will stop {{clean}} where you are, the blue arrow [[File:clean-continue.png]] 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 [[File:clean-redo.png]] 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.
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 a pop-up window (rms noise for whole cube will go to terminal). 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.  
Note that for this example, threshold has been set to <tt>threshold = '3mJy'</tt> 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 (<tt>niter = 0</tt>), 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.  


[[Image:SiS_interactive.png|thumb|Channel 16 shown for the SiS cube in the interactive viewer with the white contour showing the mask contour drawn with the polygon tool.]]
[[Image:SiS_interactive.png|thumb|Channel 43 shown for the SiS cube in the interactive viewer with the white contour showing the mask contour drawn with the polygon tool.]]


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.  
Keep cleaning, by using the green Next Action arrow [[File:clean-redo.png]] 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.
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.
Line 147: Line 123:
</source>
</source>


You can look at both cubes using the viewer, and the tape deck to play the cube as "movie".
You can look at both cubes using the viewer, and the tape deck to play the cube as a "movie".
<source lang="python">
<source lang="python">
# In CASA
# In CASA
Line 156: Line 132:


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.  
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.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
clean(vis='IRC10216.cont',imagename='IRC10216.36GHzcont',
clean(vis='IRC10216',imagename='IRC10216.36GHzcont',
       mode='mfs',imagermode='csclean',
       mode='mfs',imagermode='csclean',
       imsize=300,cell=['0.4arcsec'],spw='0~1:5~59',
       imsize=300,cell=['0.4arcsec'],spw='0~1:5~14,0~1:48~59',
       weighting='briggs',robust=0.5,
       weighting='briggs',robust=0.5,
       interactive=T)
       interactive=T)
</source>
</source>


Now look at the result in the viewer if you like
Now look at the result in the viewer, if you like:


<source lang="python">
<source lang="python">
Line 175: Line 153:
==Image Analysis and Viewing==
==Image Analysis and Viewing==


Next make integrated intensity maps (moment 0) and integrated velocity maps (moment 1). 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:
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:


<source lang="python">
<source lang="python">
Line 182: Line 160:
</source>
</source>


Then use the Viewer tape deck to see which channels have significant line emission. For HC3N, the line channel range in the cube is 11 to 40, and it is the same for SiS.  
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.  


[[Image:irc10216.jpg|thumb|HC3N moment 0 map with white continuum contours superposed.]]
[[Image:irc10216.jpg|thumb|HC3N moment 0 map with white continuum contours superposed.]]
Line 189: Line 167:
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.  
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 lets 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.
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.


<source lang="python">
<source lang="python">
Line 195: Line 173:
immoments(imagename='IRC10216_HC3N.cube_r0.5.image',moments=[0],
immoments(imagename='IRC10216_HC3N.cube_r0.5.image',moments=[0],
           axis='spectral',
           axis='spectral',
           chans='11~40',
           chans='16~45',
           outfile='IRC10216_HC3N.cube_r0.5.image.mom0')
           outfile='IRC10216_HC3N.cube_r0.5.image.mom0')
</source>
</source>
Line 203: Line 181:
immoments(imagename='IRC10216_SiS.cube_r0.5.image',moments=[0],
immoments(imagename='IRC10216_SiS.cube_r0.5.image',moments=[0],
           axis='spectral',
           axis='spectral',
           chans='11~40',
           chans='16~45',
           outfile='IRC10216_SiS.cube_r0.5.image.mom0')
           outfile='IRC10216_SiS.cube_r0.5.image.mom0')
</source>
</source>


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 5sigma.
To have a look at these, use the viewer:
 
<source lang="python">
# In CASA
viewer('IRC10216_HC3N.cube_r0.5.image.mom0')
#
viewer('IRC10216_SiS.cube_r0.5.image.mom0')
</source>
 
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&sigma;:


[[Image:irc10216_dust_HC3N_cont.png|thumb|The VLT V-band image from Le&atilde;o et al. (2006, A&A, 455,187) showing dust rings, overlaid with white HC3N moment 0 contours and black 36 GHz continuum contours.]]
[[Image:irc10216_dust_HC3N_cont.png|thumb|The VLT V-band image from Le&atilde;o et al. (2006, A&A, 455,187) showing dust rings, overlaid with white HC3N moment 0 contours and black 36 GHz continuum contours.]]
Line 215: Line 202:
immoments(imagename='IRC10216_HC3N.cube_r0.5.image',moments=[1],
immoments(imagename='IRC10216_HC3N.cube_r0.5.image',moments=[1],
           axis='spectral',
           axis='spectral',
           chans='11~40',excludepix=[-100,0.01],
           chans='16~46',excludepix=[-100,0.01],
           outfile='IRC10216_HC3N.cube_r0.5.image.mom1')
           outfile='IRC10216_HC3N.cube_r0.5.image.mom1')
</source>
</source>
Line 223: Line 210:
immoments(imagename='IRC10216_SiS.cube_r0.5.image',moments=[1],
immoments(imagename='IRC10216_SiS.cube_r0.5.image',moments=[1],
           axis='spectral',
           axis='spectral',
           chans='11~40',excludepix=[-100,0.01],
           chans='16~45',excludepix=[-100,0.01],
           outfile='IRC10216_SiS.cube_r0.5.image.mom1')
           outfile='IRC10216_SiS.cube_r0.5.image.mom1')
</source>
</source>


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:
<source lang="python">
# 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')
</source>
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.
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  
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&atilde;o and overlay the moment images and 36 GHz continuum. More information about the dust properties can be found in the Le&atilde;o et al. (2006) paper http://adsabs.harvard.edu/abs/2006A%26A...455..187L
Izan Le&atilde;o and overlay the moment images and 36 GHz continuum. More information about the dust properties can be found in the Le&atilde;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 [[How_to_rotate_and_slice_a_cube_for_pV_diagrams|pV casaguide]]. Masking the data cube to extract the emission is described [[Masking_images_for_analysis|here]].  


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 [[How_to_rotate_and_slice_a_cube_for_pV_diagrams|pV casaguide]]. Masking the data cube to extract the emission is described [[Masking_images_for_analysis|here]].


== Spectrum Fitting ==
== Spectrum Fitting ==


Frequently, one would like to fit Gaussians or polynomials to the spectral line in the data cube. CASA's {{specfit}} allows us to do this within the package. The {{specfit}} task can fit the functions to an average spectrum define by some bounding box, or, alternatively, for each pixel. In the following, we will do both.
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 ===
=== Fitting an average spectrum ===
[[Image:viewer-irc-specfit1.png|thumb|The viewer showing the cube at plane 34. The green box marks the region that is averaged for the spectral profile. ]]
[[Image:viewer-irc-specfit1-profile.png|thumb|Spectral profile of the HC3N line within the green rectangle.]]


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.
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 [[CASA Region Format | 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.
<pre style="background-color: #fffacd;">
#CRTFv0
box[[09:47:59.2, 13.16.24], [09:47:55.8, 13.17.09]]
</pre>
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: 
<source lang="python">
# In CASA
myfit = specfit(imagename='IRC10216_HC3N.cube_r0.5.image', region='specfit.crtf', multifit=F,
        estimates='', ngauss=2)


[[Image:viewer-irc-specfit1.png|thumb|The viewer showing the cube at plane 34. The green box marks the region that is averaged for the spectral profile. ]]
</source>
 
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.):


[[Image:viewer-irc-specfit1-profile.png|thumb|Spectral profile of the HC3N line within the green rectangle.]]
<pre style="background-color: #fffacd;">
Fit  :
RA          :  09:47:57.49
Dec          :  13.16.46.46
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
</pre>


To fit this profile in {{specfit}}, we need two files: 1) a region file outlining the 2-D region that is averaged (the green box) and 2) a file that contains the initial parameters for the fit.  
which seems to have caught the two peaks pretty well.


For the region file, we use the new CASA region format that is described [[CASA Region Format | here]]. Following the guidelines there, we create a file containing the following and we name it '''specfit.rgn''
=== Spectral Fitting pixel by pixel ===


=== Fitting pixel by pixel ===
[[Image:irc-amp0.png|thumb|Output image of {{specfit}} of the amplitude of the first Gaussian component.]]


{{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:


<source lang="python">
# 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')
</source>


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'''.


== Self-Calibration ==
== Self-Calibration ==


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).
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).  
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 chose 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).  
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 [http://evlaguides.nrao.edu/index.php?title=Category:Status#Sensitivity 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, one 10 second integration is given by:
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 [http://evlaguides.nrao.edu/index.php?title=Observational_Status_Summary_-_Current#Sensitivity 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:


<math> {\rm RMS(baseline)} = {\rm 2\ mJy\ beam^{-1}\ channel^{-1}} \sqrt{ \frac{95\times 60\ {\rm sec}}{10\ {\rm sec}}\times\frac{16\times 15}{2\times 1}}</math>
<math> {\rm RMS(baseline)} = {\rm 2\ mJy\ beam^{-1}\ channel^{-1}} \sqrt{ \frac{95\times 60\ {\rm sec}}{10\ {\rm sec}}\times\frac{16\times 15}{2\times 1}}\ {\sim}\ 500\ {\rm mJy\ beam^{-1}\ channel^{-1}}</math>


[[Image:timeonsource.png|thumb|Plot to estimate the time on source.]]
[[Image:timeonsource.png|thumb|Plot to estimate the time on source.]]


The 95 minutes of on-source time can be estimated from a plot like:
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:


<source lang="python">
<source lang="python">
Line 281: Line 330:
</source>
</source>


This analysis suggests that the rms noise on one baseline, one 10 second integration is about 525 mJy. In contrast the peak flux density in the strongest SiS channel is only about 200 mJy (you can check using the {{viewer}}). Because the emission is fairly compact, most baselines will see about this peak flux, this is why we chose the more compact of the two possible lines. This peak flux density tells us that 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.  
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  
Now we run {{gaincal}} with the solint we have determined. Note that because our desired solint is more  
Line 302: Line 351:
</source>
</source>


For some antennas you can see clear global trends away from zero: ea08, ea21,ea24 are examples, and you can also see some smaller variations with time.
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.
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.
Line 351: Line 400:
residuals.]]
residuals.]]


Now you can redo the moment images if you like with the improved cubes, be sure to change the output file names.  
Now you can redo the moment images if you like with the improved cubes (be sure to change the output file). names.  


[[Main Page | &#8629; '''CASAguides''']]
[[Main Page | &#8629; '''CASAguides''']]


--[[User:Cbrogan|Crystal Brogan]]  additions:[[User:Jott|Juergen Ott]]
--[[User:Cbrogan|Crystal Brogan]]   
 
--additions:[[User:Jott|Juergen Ott]], [[User:Mkrauss|Miriam Krauss]]


{{Checked 3.2.0}}
{{Checked 3.3.0}}

Latest revision as of 20:37, 29 February 2012

This is the second part of the EVLA IRC+10216 data reduction tutorial. IRC+10216 Part 1: editing, calibration

UV Continuum Subtraction and Setting Up for Self-Calibration

UV-plot of the spectral line signal in both spw for IRC+10216.

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.

uvcontsub2(vis='IRC10216',fitspw='0~1:4~13;52~60',
       want_cont=T)

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  

After the cvel step, the data can then be combined with other observations via concat and imaged in clean with mode='channel' to conserve that velocity system and grid.

Note that cvel can also Hanning smooth the data, if needed (as an alternative to the stand-alone hanningsmooth task).

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.

Channel 32 shown for the HC3N cube shown in the interactive viewer with the white contour showing the mask contour drawn with the polygon tool.

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 Clean-stop.png will stop clean where you are, the blue arrow Clean-continue.png 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 Clean-redo.png 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.

Channel 43 shown for the SiS cube in the interactive viewer with the white contour showing the mask contour drawn with the polygon tool.

Keep cleaning, by using the green Next Action arrow Clean-redo.png 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,
      threshold='3.0mJy',niter=100000)

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.

HC3N moment 0 map with white continuum contours superposed.
SiS moment 0 map with white continuum contours superposed.

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=[0],
          axis='spectral',
          chans='16~45',
          outfile='IRC10216_HC3N.cube_r0.5.image.mom0')
# In CASA
immoments(imagename='IRC10216_SiS.cube_r0.5.image',moments=[0],
          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σ:

The VLT V-band image from Leão et al. (2006, A&A, 455,187) showing dust rings, overlaid with white HC3N moment 0 contours and black 36 GHz continuum contours.
# In CASA
immoments(imagename='IRC10216_HC3N.cube_r0.5.image',moments=[1],
          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=[1],
          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.

Spectrum Fitting

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

The viewer showing the cube at plane 34. The green box marks the region that is averaged for the spectral profile.
Spectral profile of the HC3N line within the green rectangle.

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          :   13.16.46.46
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

Output image of specfit of the amplitude of the first Gaussian component.

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.

Self-Calibration

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:

[math]\displaystyle{ {\rm RMS(baseline)} = {\rm 2\ mJy\ beam^{-1}\ channel^{-1}} \sqrt{ \frac{95\times 60\ {\rm sec}}{10\ {\rm sec}}\times\frac{16\times 15}{2\times 1}}\ {\sim}\ 500\ {\rm mJy\ beam^{-1}\ channel^{-1}} }[/math]

Plot to estimate the time on source.

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~19',calmode='p',solint='10min',combine='scan',
        refant='ea02',minsnr=3.0)
Phase-only self-calibration solutions with 10 minute solint.

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',calready=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.)

Original and self-cal SiS images for channel 37, notice the decrease in residuals.

Repeat for HC3N:

Original and self-cal HC3N images for channel 13, notice the decrease in residuals.

Now you can redo the moment images if you like with the improved cubes (be sure to change the output file). names.

CASAguides

--Crystal Brogan --additions:Juergen Ott, Miriam Krauss

Last checked on CASA Version 3.3.0.