TWHydraBand7 Imaging for CASA 3.3: Difference between revisions

From CASA Guides
Jump to navigationJump to search
← Older edit
Cbrogan (talk | contribs)
operators
2,345 edits
Thunter (talk | contribs)
1,106 edits
 
(126 intermediate revisions by 6 users not shown)
Line 1: Line 1:
This tutorial picks up where [[TWHydraBand7_Calibration]] leaves off with the fully calibrated split science target ms. If you are starting directly from the TWHYA_Band7_CalibratedData.tar.gz file, untar and unzip it and if necessary rename it TWHydra_corrected.ms.  
[[Category:ALMA]][[Category:Imaging]][[Category:Spectral Line]]
* This tutorial picks up where '''[[TWHydraBand7_Calibration]]''' leaves off with the fully calibrated split science target ms. If you are starting directly from the TWHYA_Band7_CalibratedData.tar.gz file, tar -zxvf it and cd to that directory to work in.  
 
* Details of the ALMA data are provided at [[TWHydraBand7]]
 


== Prepare Continuum Data for Further Processing==
== Prepare Continuum Data for Further Processing==
Line 5: Line 9:
In order to speed up imaging/self-calibration for the continuum data we will average channels together, and then flag the few spectral features. This method works well in the case of only a few spectral lines, while retaining enough channelization to aid in multi-frequency synthesis and being able to remove spectral features in the averaged data. The amount of channel averaging should be adjusted to your particular data. '''Note for wide bandwidth data it is never a good idea to use the continuum estimate generated by uv-continuum subtraction as your continuum data'''. This worked OK in the past with the narrow band systems like the VLA and low dynamic range instruments like the SMA, but is not advised for ALMA wide bandwidth, high dynamic range data.   
In order to speed up imaging/self-calibration for the continuum data we will average channels together, and then flag the few spectral features. This method works well in the case of only a few spectral lines, while retaining enough channelization to aid in multi-frequency synthesis and being able to remove spectral features in the averaged data. The amount of channel averaging should be adjusted to your particular data. '''Note for wide bandwidth data it is never a good idea to use the continuum estimate generated by uv-continuum subtraction as your continuum data'''. This worked OK in the past with the narrow band systems like the VLA and low dynamic range instruments like the SMA, but is not advised for ALMA wide bandwidth, high dynamic range data.   


[[Image:cont_spw2.png|thumb|CO(3-2) in spw=2 from the "continuum" channel averaged data.]]
[[Image:Cont_spw2.png|thumb|CO(3-2) in spw=2 from the "continuum" channel averaged data.]]
[[Image:cont_spw3.png|thumb|Notice upswing at high channel numbers in spw=3.]]
[[Image:Cont_spw3.png|thumb|Notice upswing at high channel numbers in spw=3.]]


One can also flag the spectral features first but in that case make a backup of the file that is analogous to TWHydra_corrected.ms as you will need it later (unflagged) for spectral line imaging.
One can also flag the spectral features first but in that case make a backup of the file that is analogous to TWHydra_corrected.ms as you will need it later (unflagged) for spectral line imaging. We begin by excluding ~20 channels on each edge, and averaging over 100 channel intervals:


<source lang="python">
<source lang="python">
# In CASA
# In CASA
os.system('rm -rf TWHydra_cont.ms*')
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_cont.ms',
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_cont.ms',
       spw='0~3:22~3820',width=100,datacolumn='data')
       spw='0~3:22~3820',width=100,datacolumn='data')
Line 26: Line 31:
The CO(3-2) and HCO+(4-3) lines are obvious. The 4th spw; spw=3 also shows an upswing on the highest channels that is problematic, it was seen on the calibrators too and is likely due to a weak atmospheric feature, we will flag the worst of it.
The CO(3-2) and HCO+(4-3) lines are obvious. The 4th spw; spw=3 also shows an upswing on the highest channels that is problematic, it was seen on the calibrators too and is likely due to a weak atmospheric feature, we will flag the worst of it.


[[Image:cont_uvplot.png|thumb|UV-plot of the continuum emission, colors show spws.]]
[[Image:Cont_uvplot.png|thumb|UV-plot of the continuum emission, colors show spws.]]


<source lang="python">
<source lang="python">
Line 34: Line 39:
</source>
</source>


Have a look at the continuum as a function of uv-distance (in the plot below, the 4 spw are shown in different colors on top of each other). If the data are flat as a function of uv-distance the source is completely unresolved. If you see structure (in this case decreasing amplitude with uv-distance) the source is resolved.  
Have a look at the continuum as a function of uv-distance (in the plot to the right, the 4 spw are shown in different colors on top of each other). If the data are flat as a function of uv-distance, the source is completely unresolved. If instead you see structure (in this case decreasing amplitude with uv-distance), the source is resolved.  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotms(vis='TWHydra_cont.ms',spw='',xaxis='uvdist',yaxis='amp',field='',avgchannel='38',
plotms(vis='TWHydra_cont.ms',spw='',xaxis='uvdist',yaxis='amp',field='',avgchannel='38',
       coloraxis='spw')
       coloraxis='spw',plotfile='cont_uvplot.png')
</source>
</source>


Notice that the spw do not have the same amplitude. Unfortunately, the red spw (highest frequency) should be higher than the green one (lowest frequency) because optically thin dust goes as nu^4, and optically thick as nu^2. This inconsistency is due to imperfect amplitude/absolute flux calibration. The self-calibration below will bring them into agreement with each other, but this is the sort of thing that results in considerable absolute flux uncertainty at submillimeter wavelengths. We are working on how to do this better.
Notice that the spw do not have the same amplitude. Unfortunately, the red spw (highest frequency) should be higher than the green one (lowest frequency) because optically thin dust goes as nu^4, and optically thick as nu^2. This inconsistency is due to imperfect amplitude/absolute flux calibration. The self-calibration below will bring them into agreement with each other, but this is the sort of thing that results in considerable absolute flux uncertainty at submillimeter wavelengths. We are working on how to do this better.
==A priori Noise Estimates==
To estimate what the expected noise levels are, we first need to know how much time we had on source. Looking at
a plot of amplitude vs time for these data (see below) there are 9x8min scans + 6x6min scans + 3x1min + 2x2min scans for about 1.9 hours on source. There are 2 polarizations and about 1.75 GHz of total continuum bandwidth taking into account flagging.
Using this information and https://almascience.nrao.edu/call-for-proposals/sensitivity-calculator, we find that the '''expected rms noise for the continuum is about 0.3 mJy/beam'''.
For the spectral line data which we will image below with 0.12 km/s channels (needs to be entered in the field ''bandwidth per polarization'', the '''expected rms noise for the spectral lines is about 32 mJy/beam'''.
<gallery widths="300px" heights="180">
File:Corrected_time.png|''Plot for estimating time.
File:sensitivity_cont.png|''Sensitivity calculator setup with observed continuum parameters for these observations.
</gallery>


==Image and Self-Calibrate the Continuum==
==Image and Self-Calibrate the Continuum==


Make an initial image using {{clean}}. We use '''mode='mfs'''' to do multifrequency synthesis, in other words grid each uv-channel independently. Since the uv-spacing is a function of frequency, this will actually achieve greater uv-coverage than if all the channels and spws had been averaged in advance. The '''imsize=100''' and '''cellsize=0.3arcsec''' were chosen to be appropriate for the Band and antenna configuration. At ALMA Band 7, the Full Width Half Power primary beam is about 20", so we want to make our image at least this big. We've made it a bit larger out to about the 25% power point. For the configuration used to take these data, the resolution is expected to be about 1.5" and it is a good idea to oversample the beam by a factor of 4 or 5 (we used 5 here).   
<pre style="background-color: #E0FFFF;">
If you are unfamiliar with the CASA viewer it is a good idea to pause here and watch the demo at
http://casa.nrao.edu/CasaViewerDemo/casaViewerDemo.html
</pre>
 
====Create Initial Clean Image====
 
First we will make an initial image using {{clean}}, this task will both deconvolve and clean.
 
<pre style="background-color: #E0FFFF;">
If you are unfamiliar with the basic concepts of deconvolution and clean, pause here and
review for example http://www.aoc.nrao.edu/events/synthesis/2010/lectures/wilner_synthesis10.pdf
</pre>
 
In {{clean}}, we use '''mode='mfs'''' to do multifrequency synthesis, in other words grid each uv-channel independently. Since the uv-spacing is a function of frequency, this will actually achieve greater uv-coverage than if all the channels and spws had been averaged in advance. The '''imsize=100''' and '''cellsize=0.3arcsec''' were chosen to be appropriate for the Band and antenna configuration. At ALMA Band 7, the Full Width Half Power primary beam is about 20", so we want to make our image at least this big. We've made it a bit larger out to about the 25% power point. For the configuration used to take these data, the resolution is expected to be about 1.5" and it is a good idea to oversample the beam by setting the cell size to be a factor of 4 or 5 smaller than the beam (we have used 5 here).   


For doing self-calibration it is essential that only "real" signal be included in the image model, thus it is important to make a clean mask, and clean the first iteration rather conservatively (shallowly) so that you are sure that you do not have clean components in your model that do not represent real source emission. In the example below we set '''interactive=T''' so you can interactively make the clean mask and stop clean when the signal from the source starts to approach the surrounding "noise" areas in the residual map.
For doing self-calibration it is essential that only "real" signal be included in the image model, thus it is important to make a clean mask, and clean the first iteration rather conservatively (shallowly) so that you are sure that you do not have clean components in your model that do not represent real source emission. In the example below we set '''interactive=T''' so you can interactively make the clean mask and stop clean when the signal from the source starts to approach the surrounding "noise" areas in the residual map.


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:''' Real emission will not go away if you do not put a clean box on it. Conversely if you put a clean box on something that is NOT real and clean deeply you can amplify that feature. Thus, its best to be conservative and build up your clean box as you go, you can generally see fainter "real" features in the residuals.
 
'''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.


[[Image:viewer_clean1.png|thumb|Interactive viewer.]]
[[Image:Viewer_clean1.png|thumb|Interactive viewer.]]
'''If you are unfamiliar with the CASA viewer it is a good idea to pause here and watch the demo at http://casa.nrao.edu/CasaViewerDemo/casaViewerDemo.html'''


<source lang="python">
<source lang="python">
# In CASA
# In CASA
os.system('rm -rf TWHydra_contall.*')
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall',
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall',
       mode='mfs',calready=T,imagermode='csclean',
       mode='mfs',calready=T,imagermode='csclean',
Line 65: Line 100:
</source>
</source>


Once the Clean Viewer opens, click on the polygon tool with your left mouse button, and draw a polygon around the continuum eimssion with the left mouse button. Double click inside the region with the left mouse button when you are happy with it. It will turn from green to white when the clean mask is accepted. 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 circle arrow will clean until it reaches the "iterations" parameter on the left side of the green area. '''You should chose the green circle arrow'''. It will clean 100 iterations (number under iteration in the green area on the Viewer gui) and then show the '''residual image''' again. The middle mouse button is automatically assigned to the "plus" icon, this adjust the colorscale. '''If you cannot easily access multiple mouse buttons, you can always assign your one mouse button the icon you want by clicking on it'''. Adjust the color so you can see it better. Also note that the logger is telling the amount of clean flux recovered so far, the min/max residuals, and the synthesize beam:
'''NOTE:'''{{clean}} will generate several files automatically with the '''imagename''' prepended:
* <imagename>.image # the final restored image
* <imagename>.flux # the effective primary beam response (e.g. for pbcor)
* <imagename>.flux.pbcoverage # the primary beam coverage (ftmachine=’mosaic’ only)
* <imagename>.model # the model image
* <imagename>.residual # the residual image
* <imagename>.psf # the synthesized (dirty) beam
 
Once the Clean Viewer opens, click on the polygon tool with your left mouse button, and draw a polygon around the continuum emission with the left mouse button. Double click inside the region with the left mouse button when you are happy with it. It will turn from green to white when the clean mask is accepted. 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 number of iterations requested (niter) or the flux density threshold (whichever comes first), and the green circle arrow will clean until it reaches the "iterations" parameter on the left side of the green area. '''You should chose the green circle arrow'''. It will clean 100 iterations (number under iteration in the green area on the Viewer gui) and then show the '''residual image''' again. The middle mouse button is automatically assigned to the "plus" icon, which adjusts the colorscale. '''If you cannot easily access multiple mouse buttons, you can always assign your one mouse button the icon you want by clicking on it'''. Adjust the color scale so you can see the structure in the noise better. Also, note that the logger is telling you the amount of clean flux recovered so far, the min/max residuals, and the synthesized beam:


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
Line 75: Line 118:
Already with only 100 iterations, the residuals outside the clean mask are larger than that inside. So we will stop here by hitting the Red X symbol.
Already with only 100 iterations, the residuals outside the clean mask are larger than that inside. So we will stop here by hitting the Red X symbol.


Next we can use the clean model to self-calibrate the uv-data. The process of self-calibration tries to adjust the data to better match the model you give it. That is why it is important to make the model as good as you can. The parameter '''calready=T''' ensured that a uv-model of the resulting image was placed in the model data column. This model can then be used to self-calibrate the data using {{gaincal}}. We use '''gaintype='T'''' to tell it to average the polarizations before deriving solutions which gains us a sqrt(2) in sensitivity. A critical parameter to play with is '''solint''', we use '''solint=30s''' to gain a bit of sensitivity. From {{listobs}} the integration time of these data are about 10 seconds, using '''solint=30s''' then gives us a sqrt(3) in sensitivity. It is often best to start with a larger '''solint''' and then work your way down to the integration time if the S/N of the permits. For the same reasons as described in the [[TWHydraBand7_Calibration]] tutorial we want to fix up the phases before tackling amplitude by setting '''calcode='p''''.
====Self-calibrate the Continuum Data====
 
Next we can use the clean model to self-calibrate the uv-data. The process of self-calibration tries to adjust the data to better match the model you give it. That is why it is important to make the model as good as you can. The parameter '''calready=T''' ensured that a uv-model of the resulting image was placed in the model data column. This model can then be used to self-calibrate the data using {{gaincal}}. We use '''gaintype='T'''' to tell it to average the polarizations before deriving solutions which gains us a sqrt(2) in sensitivity. A critical parameter to play with is '''solint''', we use '''solint=30s''' to gain a bit of sensitivity. From {{listobs}} the integration time of these data are about 10 seconds, so using '''solint=30s''' then gives us a sqrt(3) in sensitivity. It is often best to start with a larger '''solint''' and then work your way down to the integration time if the S/N of the data permits. For the same reasons as described in the [[TWHydraBand7_Calibration]] tutorial we want to fix up the phases before tackling amplitude by setting '''calmode='p''''.


<source lang="python">
<source lang="python">
Line 85: Line 130:
</source>
</source>


[[Image:self_1_phase.png|thumb|Phase solutions from self_1.pcal with solint=30s.]]
[[Image:Self_1_phase.png|thumb|Phase solutions from self_1.pcal with solint=30s.]]


Because '''combine=''''', we are getting independent solutions for each spw. If you are desperate for S/N, you can try combining them with combine='spw'. We have plenty of S/N in these data so we keep them separate. Look at the solutions.
Because '''combine=' '''', we are getting independent solutions for each spw. If you are desperate for S/N, you can try combining them with combine='spw'. We have plenty of S/N in these data so we keep them separate. Look at the solutions.


<source lang="python">
<source lang="python">
Line 102: Line 147:
<source lang="python">
<source lang="python">
# In CASA
# In CASA
applycal(vis='TWHydra_cont.ms',field='',
applycal(vis='TWHydra_cont.ms',field='',gaintable=['self_1.pcal'],calwt=F)
        gaintable=['self_1.pcal'],calwt=F)
</source>
</source>


Next clean the source again, to save time we can start with the clean mask from the last iteration. If you make adjust the clean mask they will be written to the '''new''' mask TWHydra_contall_1pcal.mask, and not the old one.  
[[Image:Viewer_pcal1.png|thumb|Clean residuals after 100 iterations and the 1st phase only self-cal. Residuals inside mask are higher than the "noise" outside.]]
 
Next clean the source again, to save time we can start with the clean mask from the last iteration. If you make interactive adjustments to the clean mask, they will be written to the '''new''' mask TWHydra_contall_1pcal.mask, and not the old one.  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
os.system('rm -rf TWHydra_contall_1pcal.*')
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_1pcal',
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_1pcal',
       mode='mfs',calready=T,imagermode='csclean',
       mode='mfs',calready=T,imagermode='csclean',
Line 118: Line 165:
</source>
</source>


Now after the first 100 iterations you see that the residuals inside the clean mask are higher than those outside. This is because flux that had been spread out in the map due to poorly correlated phases has been increased where it belongs on the source. You may also notice that your original clean mask is a little too small. If so you can make it a bit larger, just by drawing again around where you want the new mask to be. When you double click it will incorporate the new mask. You might want to use the magnifying glass to zoom in. Again you assign mouse buttons but clicking with it on the icon you want. Clean another 100 iterations and stop.
Now after the first 100 iterations you see that the residuals inside the clean mask are higher than those outside. This is because flux that had been spread out in the map due to poorly correlated phases has been placed where it belongs on the source, thereby increasing its flux. You may also notice that your original clean mask is a little too small. If so, you can make it a bit larger, just by drawing a new polygon where you want the mask to extend to. When you double click it will incorporate the new mask into the current mask (if your box goes awry you can use the erase toggle [radio button] to erase and redo the mask), they need not be overlapping (though for this source that wouldn't make sense). You might also want to use the magnifying glass to zoom in. Again you assign mouse buttons by clicking on the icon you want. Clean another 100 iterations and stop (for total of 200).


Now open  
Now open initial and 1st self-cal images to compare them using {{imview}}, this is an alternative version of the {{viewer}} that is somewhat scriptable.


Now that we have a better model, we can try making the '''solint''' smaller. We'll try doing it on the
<source lang="python">
integration time by setting '''solint='int''''
# In CASA
imview(raster=[{'file':'TWHydra_contall.image','range':[-0.01,0.05]},
      {'file':'TWHydra_contall_1pcal.image','range':[-0.01,0.05]}])
</source>


[[Image:self_2_phase.png|thumb|Phase solutions from self_2.pcal with solint='int'.]]
Once the Viewer Display panel opens, click on the "blink" toggle and then use the tapedeck buttons to alternate between the two images. Because the range was explicitly set, they are on the same color scale. It is obvious that the self-calibrated image is better. To quantify, 1st toggle back to "normal" and then draw a region that excludes the central source with the polygon tool and then double click inside (if you stay in the "blink" view, only the statistics for the last image opened will appear on the terminal instead of both). The basic statistics of the polygon region will print to the '''terminal''':
<pre style="background-color: #fffacd;">
(TWHydra_contall.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean          Rms            Std dev        Minimum        Maximum
2175      -3.477803e-01  -1.022026e-02  -1.598990e-04  9.450922e-03  9.451743e-03  -2.887061e-02  3.247919e-02
 
 
(TWHydra_contall_1pcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean          Rms            Std dev        Minimum        Maximum
2175      3.493550e-01  1.026653e-02  1.606230e-04  3.792825e-03  3.790294e-03  -1.135493e-02  1.252072e-02
</pre>
 
The self-calibrated image has almost a factor of 3 lower rms noise!
 
Now that we have a better model, we solve for solutions again. This time we will try making the solint smaller -- the integration time is about 10 seconds, so we set '''solint='int'''' to have it use the integration time.
 
[[Image:Self_2_phase.png|thumb|Phase solutions from self_2.pcal with solint='int'.]]


<source lang="python">
<source lang="python">
Line 144: Line 212:
</source>
</source>


This shows a bit more structure in the phase, but there does appear to be enough S/N to use the shorter solint. It may help to zoom in on one of the plotcal panels with the magnify icon so you can see it better. Alternatively you can set '''subplot=111''' and page through each plot.
The phases on this shorter solint still seem to be tracking well so we will apply it and make a new image and check to see if there is improvement.  
 
Having decided this new iteration of phase-only self-cal was sucessful, we apply it.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
applycal(vis='TWHydra_cont.ms',field='',
applycal(vis='TWHydra_cont.ms',field='',gaintable=['self_2.pcal'],calwt=F)
        gaintable=['self_2.pcal'],calwt=F)
</source>
</source>


Now clean again, updating to start with the most recent mask file.
[[Image:Viewer_pcal2.png|thumb|Clean residuals after 2nd phase only self-cal and 250 iterations. The residuals inside the clean mask are lower than those outside, so its time to stop.]]
 
Now clean again, updating '''mask''' with the most recent mask file.


<source lang="python">
<source lang="python">
# In CASA
# In CASA
os.system('rm -rf TWHydra_contall_2pcal.*')
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_2pcal',
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_2pcal',
       mode='mfs',calready=T,imagermode='csclean',
       mode='mfs',calready=T,imagermode='csclean',
Line 166: Line 234:
</source>
</source>


[[Image:self_ap_phase.png|thumb|Residual phase solutions from self_ap.cal.]]
After 200 iterations
[[Image:self_ap_amp.png|thumb|Amplitude solutions from self_ap.cal.]]


This time clean a total of about 250 iterations before stopping (adjust niter in the gui).  
<pre style="background-color: #fffacd;">
Model 0: max, min residuals = 0.00503853, -0.00479908 clean flux 1.47338
</pre>
 
It looks like it could use just a little more cleaning, so set '''iterations=50 in the gui'''.


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
  Model 0: max, min residuals = 0.00407018, -0.00397106 clean flux 1.48179
Clean used 50 iterations to approach a threshhold of 0.00166574
0.0115307 Jy <- cleaned in this cycle for model 0
Final maximum residual = 0.00328126
Model 0: max, min residuals = 0.00328126, -0.00321554 clean flux 1.48492
</pre>
 
The small amount of flux cleaned this cycle indicates that clean is converging and it's about time to stop. Also the residuals inside the clean mask are now smaller than the residuals outside (i.e. the noise) also indicating its time to stop by hitting the Red X on the gui.
 
Now repeat the display, blinking, and image statistics steps described above for the 1st and 2nd phase-only self-cal iterations.
 
<source lang="python">
# In CASA
imview(raster=[{'file':'TWHydra_contall_1pcal.image','range':[-0.01,0.05]},
      {'file':'TWHydra_contall_2pcal.image','range':[-0.01,0.05]}])
</source>
 
The difference is small:
<pre style="background-color: #fffacd;">
(TWHydra_contall_1pcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean          Rms            Std dev        Minimum        Maximum
2511      5.658062e-01  1.662741e-02  2.253310e-04  3.828434e-03  3.822558e-03  -1.183686e-02  1.252072e-02
 
 
(TWHydra_contall_2pcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean          Rms            Std dev        Minimum        Maximum
2511      6.660774e-01  1.957409e-02  2.652638e-04  3.735169e-03  3.726479e-03  -1.143428e-02  1.235826e-02
</pre>
</pre>


This is only a small improvement suggesting we have gone as far as we can with phase self-cal. Now we will apply the best phase solution table on-the-fly while solving for the amplitude self-cal solutions. Since amplitude changes more slowly and is less constrained than phase, we use a longer '''solint='inf'''', as long as '''combine=''''' this means to do one solution per scan.  
[[Image:Self_ap_phase.png|thumb|Residual phase solutions from self_ap.cal.]]
[[Image:Self_ap_amp.png|thumb|Amplitude solutions from self_ap.cal.]]
 
This small improvement suggests that we have gone as far as we can with phase self-cal. Now we will apply the best phase solution table on-the-fly while solving for the amplitude self-cal solutions. Since amplitude changes more slowly and is less constrained than phase, we use a longer '''solint='inf'''', as long as '''combine=' '''' this means to do one solution per scan.  


<source lang="python">
<source lang="python">
Line 183: Line 284:
         refant='DV06',spw='',minblperant=4,
         refant='DV06',spw='',minblperant=4,
         gaintable=['self_2.pcal'],
         gaintable=['self_2.pcal'],
         calmode='ap',minsnr=2,)
         calmode='ap',minsnr=2)
</source>
</source>


Again its useful to look at the residual phase in the amplitude calibration table.
It's useful to look at the residual phase in the amplitude calibration table, i.e. since we applied our best phase-only solution while solving for this new amplitude and phase ('''calmode='ap'''') solution, the phase should already be corrected.


<source lang="python">
<source lang="python">
Line 195: Line 296:
</source>
</source>


The residuals are very small as expected. Now look at the amplitude solutions.
The residuals are very small as expected. (If you see large residuals here, it means that the phase-only solution is suspect, and you were mostly just moving noise around.)  Now look at the amplitude solutions.


<source lang="python">
<source lang="python">
Line 204: Line 305:
</source>
</source>


For the most part these look very good except one rather large deviation on DV06. We will need to look at this time in the calibrated data. Apply both the final phase and the amplitude calibration tables.
For the most part these look very good except one rather large deviation on DV06. Apply both the final phase and the amplitude calibration tables.
 
[[Image:Selfcal_time.png|thumb|Post self-calibration amplitude vs. time]]
[[Image:Selfcal_uvdist.png|thumb|Post self-calibration amplitude vs. uv-distance]]


<source lang="python">
<source lang="python">
# In CASA
# In CASA
applycal(vis='TWHydra_cont.ms',field='',
applycal(vis='TWHydra_cont.ms',field='',gaintable=['self_2.pcal','self_ap.cal'],calwt=F)
        gaintable=['self_2.pcal','self_ap.cal'],calwt=F)
</source>
</source>


[[Image:selfcal_uvdist.png|thumb|Post self-calibration uv-plot]]
Plot the corrected data as a function of both time and uv-distance. It is always a good idea to check these after an amplitude self-calbration to be sure that it worked properly. Occasionally, a few spurious bits of data will get blown up by the amplitude self-cal if it is not well-constrained.


Plot the corrected data as a function of uv-distance.
<source lang="python">
# In CASA
plotms(vis='TWHydra_cont.ms',xaxis='time',yaxis='amp',
      avgchannel='38',ydatacolumn='corrected',coloraxis='spw',
      plotfile='selfcal_time.png')
</source>
 
The time where we saw a large amplitude correction on DV06 looks consistent here so the large correction we saw in the solutions corrected a real issue.  


<source lang="python">
<source lang="python">
Line 223: Line 333:
</source>
</source>


You will see that this looks much better than before.  
You will see that this plot looks much better than before.
 
[[Image:Viewer_apcal.png|thumb|Clean residuals after final amplitude and phase self-cal.]]


Now make the final continuum image.
Now make the final continuum image.
Line 229: Line 341:
<source lang="python">
<source lang="python">
# In CASA
# In CASA
os.system('rm -rf TWHydra_contall_apcal.*')
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_apcal',
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_apcal',
       mode='mfs',calready=T,imagermode='csclean',
       mode='mfs',calready=T,imagermode='csclean',
Line 237: Line 350:
</source>
</source>


clean until the region inside the mask is noise like.
We know we need to clean at least 200 iterations, so set this in the gui and hit green arrow. After this notice how much better the residual looks. Clean another 100 iterations. The "circular pattern", and 4 bright spots at each corner you see now in the residual is likely due to non-closing amplitude errors that cannot be removed with an antenna based self-calibration, so this is as good as it gets with pure self-cal. Hit the Red X to stop clean.
 
<pre style="background-color: #fffacd;">
Clean used 300 iterations to approach a threshhold of 0.00105507
0.0061808 Jy <- cleaned in this cycle for  model 0
Final maximum residual = 0.00242926
Model 0: max, min residuals = 0.00238278, -0.00242926 clean flux 1.48201
</pre>
 
Now repeat the display, blinking, and image statistics described above for the 2nd phase only and the new amplitude & phase self-cal'ed images.
 
<source lang="python">
# In CASA
imview(raster=[{'file':'TWHydra_contall_2pcal.image','range':[-0.01,0.05]},
      {'file':'TWHydra_contall_apcal.image','range':[-0.01,0.05]}])
</source>
 
The amplitude and phase self cal'ed image is notably better than phase-only. Quantitatively:
<pre style="background-color: #fffacd;">
(TWHydra_contall_2pcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean          Rms            Std dev        Minimum        Maximum
2621      8.885072e-01  2.611066e-02  3.389955e-04  3.717286e-03  3.702503e-03  -1.089106e-02  1.235826e-02
 
 
(TWHydra_contall_apcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean          Rms            Std dev        Minimum        Maximum
2621      2.142941e-01  6.297484e-03  8.176043e-05  1.660546e-03  1.658849e-03  -5.046569e-03  5.200869e-03
</pre>
 
and indeed, we've gained another factor of 2.2 in sensitivity.
 
'''Overall with this self-calibration we've improved the noise by a factor of 5.7! The peak flux density also increased from 0.96 to 1.0 Jy, for a change in the S/N from 101 initially to just over 600 after. So well worth the trouble.'''
 
It is notable that we are still a factor of more than 5 from the Apriori noise estimate of 0.3 mJy/beam. This is likely due to residual sources of non-closing (not antenna based) errors in the data like baseline errors and baseline-dependent bandpass errors. Additionally, polarization errors can contribute (which we did not calibrate at all), as well as imperfect uv-coverage with only 8 antennas. A dynamic range of 600 is actually quite good in the submillimeter.


==Apply self-calibration to Full Dataset and Split Line Data==
==Apply self-calibration to Full Dataset and Split Line Data==
Line 252: Line 400:
<source lang="python">
<source lang="python">
# In CASA
# In CASA
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_CO3_2.ms',
os.system('rm -rf TWHydra_CO3_2.ms*')
      datacolumn='corrected',spw='2')
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_CO3_2.ms',datacolumn='corrected',spw='2')
</source>
</source>


<source lang="python">
<source lang="python">
# In CASA
# In CASA
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_HCOplus.ms',
os.system('rm -rf TWHydra_HCOplus.ms*')
      datacolumn='corrected',spw='0')
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_HCOplus.ms',datacolumn='corrected',spw='0')
</source>
</source>


==Continuum Subtraction and Spectral Line Imaging==
==Continuum Subtraction==  


Make channel plot to get line-free channel ranges. Since the line emission is simple in this case and there are more than enough line-free channels we do not need to be very precise just stay well away from the line and best to chose ranges on either side.
Next we need to subtract the continuum emission from the spectral line data. It is best to do this in the uv-plane. The first step is to make a channel plot of the spectral data to get the line-free channel ranges for both spectral lines. Since the line emission is simple (one spectral line per spw) and there are more than enough line-free channels we do not need to be very precise just stay well away from the line and whenever possible it is best to chose line-free channel ranges on either side of the line for the best possible baseline subtraction.  


[[Image:CO3_2_channel.png|thumb|CO(3-2) line in channel space.]]
[[Image:HCOp4_3_channel.png|thumb|HCO+(4-3) in channel space.]]


'''Note''': after the split above, the original spw ids will be relabeled as 0 in each file because there is only one spectral window in each file.
<source lang="python">
# In CASA
plotms(vis='TWHydra_CO3_2.ms',spw='0',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T,coloraxis='spw',plotfile='CO3_2_channel.png')
</source>
It is a good idea to stay a little away from the very edges of the baseband which can be noisy. The channel ranges
20~2000 and 2400~3800 look like good choices for CO(3-2).
<source lang="python">
# In CASA
plotms(vis='TWHydra_HCOplus.ms',spw='0',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T,coloraxis='spw',plotfile='HCOp4_3_channel.png')
</source>
The channel ranges 20~1500 and 1900~3800 look like good choices for HCO+(4-3).
Next use the identified line-free channel ranges and {{uvcontsub2}} to subtract the continuum model from the line emission.


<source lang="python">
<source lang="python">
Line 275: Line 445:
<source lang="python">
<source lang="python">
# In CASA
# In CASA
uvcontsub2(vis='TWHydra_HCOplus.ms',fitorder=1,fitspw='0:20~1500,0:1900~3800')
</source>
Note: If you do not have line-free channels on both sides of the line emission, it is safer to use fitorder=0, or the fit may diverge.
==Spectral Line Imaging==
[[Image:CO3_2_vel.png|thumb|CO(3-2) spectrum in LSRK velocity space.]]
[[Image:HCOp4_3_vel.png|thumb|HCO+(4-3) spectrum in LSRK velocity space.]]
Now we are ready to make cubes of the line emission. The imaging parameters are similar to the continuum except for
those dealing with the spectral setup: '''mode, start, width, nchan, restfreq, and outframe''' parameters. When making spectral images you have three choices for the '''mode''' parameter: '''channel''', '''velocity''', and '''frequency'''. Data are taken using constant frequency channels. For spectral line analysis it's often more useful to have constant velocity channels, and this is also the best way to make images of multiple lines with the exact same channelization for later comparison. For '''mode='velocity'''', the desired '''start''' and '''width''' also need to be given in velocity units for the desired output frame.
It is important to note that ALMA does not do on-line Doppler Tracking and the native frame of the data is TOPO. If you do not specify '''outframe''' the output cube will also be in TOPO, which is not very useful for spectral line work. The Doppler Shift is taken out during the regridding to the desired outframe in {{clean}} or alternatively it can be done separately by the {{cvel}} task which would need to be run before {{clean}}.
To see what velocity parameters you want to set in clean it is useful to make a plot in {{plotms}} in the desired output frame. To make it easier to see, we set a '''plotrange''' since we know the LSRK velocity of TW Hya is about 2.88 km/s. Note these plots take a little longer because of the frame shift.
<source lang="python">
# In CASA
plotms(vis='TWHydra_CO3_2.ms.contsub',xaxis='velocity',yaxis='amp',
      avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',
      restfreq='345.79599GHz',plotrange=[-20,23,0,0],plotfile='CO3_2_vel.png')
</source>
<source lang="python">
# In CASA
plotms(vis='TWHydra_HCOplus.ms.contsub',xaxis='velocity',yaxis='amp',
      avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',
      restfreq='356.7342GHz',plotrange=[-20,23,0,0],plotfile='HCOp4_3_vel.png')
</source>
From these plots it looks like imaging the cubes from about -5 to +8 km/s will encompass the line, but still provide  several line-free channels on either side so the noise level can be estimated. The channel width is 122 kHz, which at
345.79599 GHz is 0.106 km/s. Recall that the spectral ''resolution'' is a factor of two poorer. We will use a velocity channel width of 0.12 km/s for a little oversampling, but the spectral resolution remains 0.2 km/s.
[[Image:Viewer_55.png|thumb|Interactive Viewer window with CO(3-2) clean mask, before cleaning starts.]]
[[Image:Viewer_500.png|thumb|Interactive Viewer window for CO(3-2), after 500 iterations, more cleaning is needed.]]
As before, it is very important that you make a clean mask. There are many ways to do this ranging from the complicated to simple. For this example we will make a single clean mask that encompasses the line emission in every channel and apply it to all channels. This is much better than no clean mask, though not quite as good as making an individual mask for each channel. For some sources, one can use the continuum mask to clean the line, however TW Hya is a bit curious in that the continuum emission is not as extended as the line emission (Hughes et al. 2008).
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_CO3_2line.*')
clean(vis='TWHydra_CO3_2.ms.contsub',imagename='TWHydra_CO3_2line',
clean(vis='TWHydra_CO3_2.ms.contsub',imagename='TWHydra_CO3_2line',
       imagermode='csclean',calready=T,spw='',
       imagermode='csclean',calready=T,spw='',
Line 285: Line 497:
</source>
</source>


The viewer can be used to make nice figures. In this example we opened the CO(3-2) cube as a raster image and then the continuum as a contour image. Then we used the "wrench" icon and "P wrench" icons to set up the channel images, contour levels etc.
It will take a little while to grid the cube, when the interactive viewer opens, use the tapedeck to cycle through the channels. You should see a progress of line emission in channels 40 to 70 or so.
*Go to channel 57 (about line center)
*Select the polygon tool with the left mouse button
*Select the "All Channels" toggle in the green area
*Draw a polygon around the channel 57 emission, double click inside to activate.
*Now go to channel 60, this is the channel with the most extent to the SE, draw another polygon that encompasses this emission if necessary,double click inside to activate.
*Now go to channel 55, this is the channel with the most extent to the NW, draw another polygon that encompasses this emission if necessary,double click inside to activate.
*Now double check the mask by cycling through channels with real line emission to be sure that your clean box encompasses the emission, but does not extend far beyond it. Again the "erase" toggle can come in handy if things go awry or if you want to just delete a small portion of an existing mask (just draw a polygon around what you want to erase and double click inside), be sure to switch back to "add" before going on.


[[Image:TWHya_channel_co3_2.png|center|800px]]''Channel maps of the CO(3-2) emission with white continuum contours at 3 and 100 sigma.''
The figure shows what the final clean box should look like. Now clean using the green circle arrow. You can watch the progress in the logger. When 100 iterations are done, the viewer will show the residual  map for each channel. Cycle through the channels and see whether you are still happy with the clean box in each signal channel. The "erase button" can help you fix mistakes. If necessary adjust. It is often useful to adjust the colorscale with the "plus" symbol icon. Keep going until the residual emission looks like the surrounding noise.  If it seems to be going slowly, increase the iterations to 300 for a round or two. I stopped cleaning (Red x) when the logger showed (about 900 iterations).
<pre style="background-color: #fffacd;">
12.9846 Jy <- cleaned in this cycle for  model 0
Final maximum residual = 0.157312
Model 0: max, min residuals = 0.157312, -0.0871295 clean flux 279.304
</pre>
 
NOTE: The residuals are still a bit high in the strongest emission channels after it is well cleaned in the weaker channels. This is an indication that the sensitivity is "dynamic range" limited. In other words even with more
intrinsic sensitivity per integration (i.e. better system temperature), the dynamic range limit would remain the same unless you get more uv-coverage. The fundamental dynamic range limitations of ALMA will be considerably better with 16 antennas in Early Science and MUCH better with the full 50-antenna array.


Now repeat for the HCO+(4-3) line.
----


<source lang="python">
Now image and clean the HCO+(4-3) line using what you learned above.
# In CASA
uvcontsub2(vis='TWHydra_HCOplus.ms',fitorder=1,fitspw='0:20~1500,0:1900~3800')
</source>


<source lang="python">
<source lang="python">
# In CASA
# In CASA
os.system('rm -rf TWHydra_HCOplusline.*')
clean(vis='TWHydra_HCOplus.ms.contsub',imagename='TWHydra_HCOplusline',
clean(vis='TWHydra_HCOplus.ms.contsub',imagename='TWHydra_HCOplusline',
       imagermode='csclean',calready=T,spw='',
       imagermode='csclean',calready=T,spw='',
Line 308: Line 533:
</source>
</source>


[[Image:TWHya_channel_HCOp4_3.png|center|800px]]''Channel maps of the HCO+(4-3) emission with white continuum contours at 3 and 100 sigma. The color intensity scale is the same as that used for CO(3-2).''
I stopped cleaning when:
<pre style="background-color: #fffacd;">
4.79833 Jy <- cleaned in this cycle for  model 0
Final maximum residual = 0.0977348
Model 0: max, min residuals = 0.0977348, -0.060484 clean flux 192.066
</pre>
 
==Image Analysis==
 
====Moment Maps====


==Further Image Analysis==
Next we will make moment maps for the image cubes.


To determine channel ranges it is useful to open cube in viewer as a  
To determine channel ranges it is useful to open the cube in {{viewer}} as a  
raster. Using the statistics tool the rms in a line free channel is  
raster image. Using the statistics tool the rms in a line-free channel is  
about 25 mJy/beam. Now open the cube again but as a contour map and  
about 25 mJy/beam, however in one of the peak line channels it is as high as 50 mJy/beam. This difference is again because of the dynamic range limitations of the data (as described above for the continuum). Now open the cube again but as a contour map and set the Unit Contour Level to 0.025 (this is in Jy/beam) and  
set the Unit Contour Level to 0.025 (this is in Jy/beam) and  
in Relative Contour levels set [3]. i.e. this is showing 3sigma contours.
in Relative Contour levels set [3]. i.e. this is showing 3sigma contours.
Now find out where the CO(3-2) emission drops below the 3sigma level.
Now find out where the CO(3-2) emission drops below the 3sigma level.
Line 330: Line 563:


For higher order moments it is very important to set a conservative flux  
For higher order moments it is very important to set a conservative flux  
threshold. Typically something like 10sigma works well.
threshold. Typically something like 6sigma, using sigma from peak line channel works well.


<source lang="python">
<source lang="python">
Line 336: Line 569:
immoments(imagename='TWHydra_CO3_2line.image',moments=[1,2],
immoments(imagename='TWHydra_CO3_2line.image',moments=[1,2],
           outfile='TWHydra_CO3_2line.image.mom',
           outfile='TWHydra_CO3_2line.image.mom',
           chans='40~76',includepix=[0.25,100])
           chans='40~76',includepix=[0.3,100])
</source>
</source>
Display all three moment maps in the viewer, overlaid with continuum contours.
<source lang="python">
# In CASA
imview( raster=[ {'file':'TWHydra_CO3_2line.image.mom0'},
                {'file':'TWHydra_CO3_2line.image.mom.weighted_coord'},
                {'file':'TWHydra_CO3_2line.image.mom.weighted_dispersion_coord'} ],
        contour={'file':'TWHydra_contall_apcal.image', 'base':0, 'unit':0.0025, 'levels':[3,100]} )
</source>
To see all three raster images simultaneously, open the viewer's Panel Display Options (the "P-wrench" icon in the viewer tool bar), click "Number of panels", and change "Number of panels in x" to 3.  Then select the "Blink" toggle (radio button) from the Animator panel.


[[Image:TWHya_CO3_2_moments.png|center|800px]] ''ALMA CO(3-2) moment maps,  with white continuum contours at 3 and 100 sigma. From left to right: integrated intensity, intensity weighted velocity field, intensity weighted velocity dispersion are shown.''
[[Image:TWHya_CO3_2_moments.png|center|800px]] ''ALMA CO(3-2) moment maps,  with white continuum contours at 3 and 100 sigma. From left to right: integrated intensity, intensity weighted velocity field, intensity weighted velocity dispersion are shown.''
Line 346: Line 591:
<source lang="python">
<source lang="python">
# In CASA
# In CASA
immoments(imagename='TWHydra_CO3_2line.image',moments=[0],
immoments(imagename='TWHydra_HCOplusline.image',moments=[0],
           outfile='TWHydra_CO3_2line.image.mom0',
           outfile='TWHydra_HCOplusline.image.mom0',
           chans='43~74')
           chans='43~74')
</source>
</source>
Line 353: Line 598:
<source lang="python">
<source lang="python">
# In CASA
# In CASA
immoments(imagename='TWHydra_CO3_2line.image',moments=[1,2],
immoments(imagename='TWHydra_HCOplusline.image',moments=[1,2],
           outfile='TWHydra_CO3_2line.image.mom',
           outfile='TWHydra_HCOplusline.image.mom',
           chans='43~74',includepix=[0.25,100])
           chans='43~74',includepix=[0.3,100])
</source>
</source>
Display all three moment maps in the viewer using the same method described above for the CO(3-2) maps.


[[Image:TWHya_HCOp4_3_moments.png|center|800px]] ''ALMA HCO+(4-3) moment maps,  with white continuum contours at 3 and 100 sigma. From left to right: integrated intensity, intensity weighted velocity field, intensity weighted velocity dispersion are shown.''
[[Image:TWHya_HCOp4_3_moments.png|center|800px]] ''ALMA HCO+(4-3) moment maps,  with white continuum contours at 3 and 100 sigma. From left to right: integrated intensity, intensity weighted velocity field, intensity weighted velocity dispersion are shown.''
====Channel Maps====
Using the Viewer you can make channel map figures. Start the Viewer and then open the CO(3-2) cube as a raster image and then the continuum as a contour image. Then we use the "wrench" icon and "P wrench" icons to set up the channel images, contour levels etc.
[[Image:TWHya_channel_co3_2.png|center|800px]]''Channel maps of the CO(3-2) emission with white continuum contours at 3 and 100 sigma.''
[[Image:Viewer_setup.png|center|500px]]''The setup windows looked like this''
Repeat for HCO+(4-3)
[[Image:TWHya_channel_HCOp4_3.png|center|800px]]''Channel maps of the HCO+(4-3) emission with white continuum contours at 3 and 100 sigma. The color intensity scale is the same as that used for CO(3-2).''
====Exporting Fits Images====
If you want to analyze the data using another software package it is easy to convert from CASA format to FITS.
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_CO3_2line.image.fits')
exportfits(imagename='TWHydra_CO3_2line.image',fitsimage='TWHydra_CO3_2line.image.fits')
</source>
Although "fits format" is supposed to be a standard, in fact most packages expect slightly different things from a fits image. If you are having difficulty, try setting '''velocity=T''' and/or '''dropstokes=T'''.
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_HCOplusline.image.fits')
exportfits(imagename='TWHydra_HCOplusline.image',fitsimage='TWHydra_HCOplusline.image.fits')
</source>
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_CO3_2line.image.mom0.fits')
exportfits(imagename='TWHydra_CO3_2line.image.mom0',
        fitsimage='TWHydra_CO3_2line.image.mom0.fits')
</source>
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_CO3_2line.image.mom.weighted_coord.fits')
exportfits(imagename='TWHydra_CO3_2line.image.mom.weighted_coord',
        fitsimage='TWHydra_CO3_2line.image.mom.weighted_coord.fits')
</source>
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_CO3_2line.image.mom.weighted_dispersion_coord.fits')
exportfits(imagename='TWHydra_CO3_2line.image.mom.weighted_dispersion_coord',
        fitsimage='TWHydra_CO3_2line.image.mom.weighted_dispersion_coord.fits')
</source>
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_HCOplusline.image.mom0.fits')
exportfits(imagename='TWHydra_HCOplusline.image.mom0',
        fitsimage='TWHydra_HCOplusline.image.mom0.fits')
</source>
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_HCOplusline.image.mom.weighted_coord.fits')
exportfits(imagename='TWHydra_HCOplusline.image.mom.weighted_coord',
        fitsimage='TWHydra_HCOplusline.image.mom.weighted_coord.fits')
</source>
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_HCOplusline.image.mom.weighted_dispersion_coord.fits')
exportfits(imagename='TWHydra_HCOplusline.image.mom.weighted_dispersion_coord',
        fitsimage='TWHydra_HCOplusline.image.mom.weighted_dispersion_coord.fits')
</source>
<source lang="python">
# In CASA
os.system('rm -rf TWHydra_contall_apcal.image.fits')
exportfits(imagename='TWHydra_contall_apcal.image',fitsimage='TWHydra_contall_apcal.image.fits')
</source>
{{Checked 3.3.0}}

Latest revision as of 01:20, 24 April 2012

  • This tutorial picks up where TWHydraBand7_Calibration leaves off with the fully calibrated split science target ms. If you are starting directly from the TWHYA_Band7_CalibratedData.tar.gz file, tar -zxvf it and cd to that directory to work in.


Prepare Continuum Data for Further Processing

In order to speed up imaging/self-calibration for the continuum data we will average channels together, and then flag the few spectral features. This method works well in the case of only a few spectral lines, while retaining enough channelization to aid in multi-frequency synthesis and being able to remove spectral features in the averaged data. The amount of channel averaging should be adjusted to your particular data. Note for wide bandwidth data it is never a good idea to use the continuum estimate generated by uv-continuum subtraction as your continuum data. This worked OK in the past with the narrow band systems like the VLA and low dynamic range instruments like the SMA, but is not advised for ALMA wide bandwidth, high dynamic range data.

CO(3-2) in spw=2 from the "continuum" channel averaged data.
Notice upswing at high channel numbers in spw=3.

One can also flag the spectral features first but in that case make a backup of the file that is analogous to TWHydra_corrected.ms as you will need it later (unflagged) for spectral line imaging. We begin by excluding ~20 channels on each edge, and averaging over 100 channel intervals: 

# In CASA
os.system('rm -rf TWHydra_cont.ms*')
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_cont.ms',
       spw='0~3:22~3820',width=100,datacolumn='data')

Now make a plot of amplitude vs. channel to see what needs to be flagged 

# In CASA
plotms(vis='TWHydra_cont.ms',spw='0~3',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T,coloraxis='spw',iteraxis='spw',xselfscale=T)

The CO(3-2) and HCO+(4-3) lines are obvious. The 4th spw; spw=3 also shows an upswing on the highest channels that is problematic, it was seen on the calibrators too and is likely due to a weak atmospheric feature, we will flag the worst of it.

UV-plot of the continuum emission, colors show spws.
# In CASA
flagdata(vis='TWHydra_cont.ms',
         spw=['0:16~16','2:21~21','3:33~37'])

Have a look at the continuum as a function of uv-distance (in the plot to the right, the 4 spw are shown in different colors on top of each other). If the data are flat as a function of uv-distance, the source is completely unresolved. If instead you see structure (in this case decreasing amplitude with uv-distance), the source is resolved. 

# In CASA
plotms(vis='TWHydra_cont.ms',spw='',xaxis='uvdist',yaxis='amp',field='',avgchannel='38',
       coloraxis='spw',plotfile='cont_uvplot.png')

Notice that the spw do not have the same amplitude. Unfortunately, the red spw (highest frequency) should be higher than the green one (lowest frequency) because optically thin dust goes as nu^4, and optically thick as nu^2. This inconsistency is due to imperfect amplitude/absolute flux calibration. The self-calibration below will bring them into agreement with each other, but this is the sort of thing that results in considerable absolute flux uncertainty at submillimeter wavelengths. We are working on how to do this better.

A priori Noise Estimates

To estimate what the expected noise levels are, we first need to know how much time we had on source. Looking at a plot of amplitude vs time for these data (see below) there are 9x8min scans + 6x6min scans + 3x1min + 2x2min scans for about 1.9 hours on source. There are 2 polarizations and about 1.75 GHz of total continuum bandwidth taking into account flagging.

Using this information and https://almascience.nrao.edu/call-for-proposals/sensitivity-calculator, we find that the expected rms noise for the continuum is about 0.3 mJy/beam.

For the spectral line data which we will image below with 0.12 km/s channels (needs to be entered in the field bandwidth per polarization, the expected rms noise for the spectral lines is about 32 mJy/beam.

  • Plot for estimating time.
    Plot for estimating time.
  • Sensitivity calculator setup with observed continuum parameters for these observations.
    Sensitivity calculator setup with observed continuum parameters for these observations.

Image and Self-Calibrate the Continuum

If you are unfamiliar with the CASA viewer it is a good idea to pause here and watch the demo at 
http://casa.nrao.edu/CasaViewerDemo/casaViewerDemo.html

Create Initial Clean Image

First we will make an initial image using clean, this task will both deconvolve and clean. 

If you are unfamiliar with the basic concepts of deconvolution and clean, pause here and 
review for example http://www.aoc.nrao.edu/events/synthesis/2010/lectures/wilner_synthesis10.pdf

In clean, we use mode='mfs' to do multifrequency synthesis, in other words grid each uv-channel independently. Since the uv-spacing is a function of frequency, this will actually achieve greater uv-coverage than if all the channels and spws had been averaged in advance. The imsize=100 and cellsize=0.3arcsec were chosen to be appropriate for the Band and antenna configuration. At ALMA Band 7, the Full Width Half Power primary beam is about 20", so we want to make our image at least this big. We've made it a bit larger out to about the 25% power point. For the configuration used to take these data, the resolution is expected to be about 1.5" and it is a good idea to oversample the beam by setting the cell size to be a factor of 4 or 5 smaller than the beam (we have used 5 here).

For doing self-calibration it is essential that only "real" signal be included in the image model, thus it is important to make a clean mask, and clean the first iteration rather conservatively (shallowly) so that you are sure that you do not have clean components in your model that do not represent real source emission. In the example below we set interactive=T so you can interactively make the clean mask and stop clean when the signal from the source starts to approach the surrounding "noise" areas in the residual map.

NOTE: Real emission will not go away if you do not put a clean box on it. Conversely if you put a clean box on something that is NOT real and clean deeply you can amplify that feature. Thus, its best to be conservative and build up your clean box as you go, you can generally see fainter "real" features in the residuals.

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.

Interactive viewer.
# In CASA
os.system('rm -rf TWHydra_contall.*')
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall',
      mode='mfs',calready=T,imagermode='csclean',
      imsize=100,cell=['0.3arcsec'],spw='',
      weighting='briggs',robust=0.5, 
      mask='',
      interactive=T,threshold='1mJy',niter=10000)

NOTE:clean will generate several files automatically with the imagename prepended:

  • <imagename>.image # the final restored image
  • <imagename>.flux # the effective primary beam response (e.g. for pbcor)
  • <imagename>.flux.pbcoverage # the primary beam coverage (ftmachine=’mosaic’ only)
  • <imagename>.model # the model image
  • <imagename>.residual # the residual image
  • <imagename>.psf # the synthesized (dirty) beam

Once the Clean Viewer opens, click on the polygon tool with your left mouse button, and draw a polygon around the continuum emission with the left mouse button. Double click inside the region with the left mouse button when you are happy with it. It will turn from green to white when the clean mask is accepted. 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 number of iterations requested (niter) or the flux density threshold (whichever comes first), and the green circle arrow will clean until it reaches the "iterations" parameter on the left side of the green area. You should chose the green circle arrow. It will clean 100 iterations (number under iteration in the green area on the Viewer gui) and then show the residual image again. The middle mouse button is automatically assigned to the "plus" icon, which adjusts the colorscale. If you cannot easily access multiple mouse buttons, you can always assign your one mouse button the icon you want by clicking on it. Adjust the color scale so you can see the structure in the noise better. Also, note that the logger is telling you the amount of clean flux recovered so far, the min/max residuals, and the synthesized beam: 

Model 0: max, min residuals = 0.0203202, 0.00159066 clean flux 1.31662
Threshhold not reached yet.
Beam used in restoration: 1.69902 by 1.59083 (arcsec) at pa 22.3441 (deg)

Already with only 100 iterations, the residuals outside the clean mask are larger than that inside. So we will stop here by hitting the Red X symbol.

Self-calibrate the Continuum Data

Next we can use the clean model to self-calibrate the uv-data. The process of self-calibration tries to adjust the data to better match the model you give it. That is why it is important to make the model as good as you can. The parameter calready=T ensured that a uv-model of the resulting image was placed in the model data column. This model can then be used to self-calibrate the data using gaincal. We use gaintype='T' to tell it to average the polarizations before deriving solutions which gains us a sqrt(2) in sensitivity. A critical parameter to play with is solint, we use solint=30s to gain a bit of sensitivity. From listobs the integration time of these data are about 10 seconds, so using solint=30s then gives us a sqrt(3) in sensitivity. It is often best to start with a larger solint and then work your way down to the integration time if the S/N of the data permits. For the same reasons as described in the TWHydraBand7_Calibration tutorial we want to fix up the phases before tackling amplitude by setting calmode='p'. 

# In CASA
gaincal(vis='TWHydra_cont.ms',caltable='self_1.pcal',
        solint='30s',combine='',gaintype='T',
        refant='DV06',spw='',minblperant=4,
        calmode='p',minsnr=2)
Phase solutions from self_1.pcal with solint=30s.

Because combine=' ', we are getting independent solutions for each spw. If you are desperate for S/N, you can try combining them with combine='spw'. We have plenty of S/N in these data so we keep them separate. Look at the solutions. 

# In CASA
plotcal(caltable='self_1.pcal',xaxis='time',yaxis='phase',
        spw='',field='',antenna='1~8',iteration='antenna',
        subplot=421,plotrange=[0,0,-80,80],figfile='self_1_phase.png')

The magnitude of corrections gives you a sense of how accurate the phase transfer from the calibrators was. You are checking to see that the phases are being well tracked by the solutions, i.e. smoothly varying functions of time within a given scan. Also the solutions for each spw are lying basically on top of each other which indicates that the spw to spw phase offsets are well calibrated already. Overall, these solutions look quite good. If they didn't you could try increasing the solint and/or combining the spw into a single solution.

After checking that you are happy with the solutions, apply them with applycal 

# In CASA
applycal(vis='TWHydra_cont.ms',field='',gaintable=['self_1.pcal'],calwt=F)
Clean residuals after 100 iterations and the 1st phase only self-cal. Residuals inside mask are higher than the "noise" outside.

Next clean the source again, to save time we can start with the clean mask from the last iteration. If you make interactive adjustments to the clean mask, they will be written to the new mask TWHydra_contall_1pcal.mask, and not the old one. 

# In CASA
os.system('rm -rf TWHydra_contall_1pcal.*')
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_1pcal',
      mode='mfs',calready=T,imagermode='csclean',
      imsize=100,cell=['0.3arcsec'],spw='',
      weighting='briggs',robust=0.5, 
      mask='TWHydra_contall.mask',
      interactive=T,threshold='1mJy',niter=10000)

Now after the first 100 iterations you see that the residuals inside the clean mask are higher than those outside. This is because flux that had been spread out in the map due to poorly correlated phases has been placed where it belongs on the source, thereby increasing its flux. You may also notice that your original clean mask is a little too small. If so, you can make it a bit larger, just by drawing a new polygon where you want the mask to extend to. When you double click it will incorporate the new mask into the current mask (if your box goes awry you can use the erase toggle [radio button] to erase and redo the mask), they need not be overlapping (though for this source that wouldn't make sense). You might also want to use the magnifying glass to zoom in. Again you assign mouse buttons by clicking on the icon you want. Clean another 100 iterations and stop (for total of 200).

Now open initial and 1st self-cal images to compare them using imview, this is an alternative version of the viewer that is somewhat scriptable. 

# In CASA
imview(raster=[{'file':'TWHydra_contall.image','range':[-0.01,0.05]},
       {'file':'TWHydra_contall_1pcal.image','range':[-0.01,0.05]}])

Once the Viewer Display panel opens, click on the "blink" toggle and then use the tapedeck buttons to alternate between the two images. Because the range was explicitly set, they are on the same color scale. It is obvious that the self-calibrated image is better. To quantify, 1st toggle back to "normal" and then draw a region that excludes the central source with the polygon tool and then double click inside (if you stay in the "blink" view, only the statistics for the last image opened will appear on the terminal instead of both). The basic statistics of the polygon region will print to the terminal: 

(TWHydra_contall.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean           Rms            Std dev        Minimum        Maximum
2175      -3.477803e-01  -1.022026e-02  -1.598990e-04  9.450922e-03   9.451743e-03   -2.887061e-02  3.247919e-02


(TWHydra_contall_1pcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean           Rms            Std dev        Minimum        Maximum
2175      3.493550e-01   1.026653e-02   1.606230e-04   3.792825e-03   3.790294e-03   -1.135493e-02  1.252072e-02

The self-calibrated image has almost a factor of 3 lower rms noise!

Now that we have a better model, we solve for solutions again. This time we will try making the solint smaller -- the integration time is about 10 seconds, so we set solint='int' to have it use the integration time.

Phase solutions from self_2.pcal with solint='int'.
# In CASA
gaincal(vis='TWHydra_cont.ms',caltable='self_2.pcal',
        solint='int',combine='',gaintype='T',
        refant='DV06',spw='',minblperant=4,
        calmode='p',minsnr=2)

Now plot the new solutions. 

# In CASA
plotcal(caltable='self_2.pcal',xaxis='time',yaxis='phase',
        spw='',field='',antenna='1~8',iteration='antenna',
        subplot=421,plotrange=[0,0,-80,80],figfile='self_2_phase.png')

The phases on this shorter solint still seem to be tracking well so we will apply it and make a new image and check to see if there is improvement. 

# In CASA
applycal(vis='TWHydra_cont.ms',field='',gaintable=['self_2.pcal'],calwt=F)
Clean residuals after 2nd phase only self-cal and 250 iterations. The residuals inside the clean mask are lower than those outside, so its time to stop.

Now clean again, updating mask with the most recent mask file. 

# In CASA
os.system('rm -rf TWHydra_contall_2pcal.*')
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_2pcal',
      mode='mfs',calready=T,imagermode='csclean',
      imsize=100,cell=['0.3arcsec'],spw='',
      weighting='briggs',robust=0.5, 
      mask='TWHydra_contall_1pcal.mask',
      interactive=T,threshold='1mJy',niter=10000)

After 200 iterations 

Model 0: max, min residuals = 0.00503853, -0.00479908 clean flux 1.47338

It looks like it could use just a little more cleaning, so set iterations=50 in the gui. 

Clean used 50 iterations to approach a threshhold of 0.00166574
0.0115307 Jy <- cleaned in this cycle for  model 0
Final maximum residual = 0.00328126
Model 0: max, min residuals = 0.00328126, -0.00321554 clean flux 1.48492

The small amount of flux cleaned this cycle indicates that clean is converging and it's about time to stop. Also the residuals inside the clean mask are now smaller than the residuals outside (i.e. the noise) also indicating its time to stop by hitting the Red X on the gui.

Now repeat the display, blinking, and image statistics steps described above for the 1st and 2nd phase-only self-cal iterations. 

# In CASA
imview(raster=[{'file':'TWHydra_contall_1pcal.image','range':[-0.01,0.05]},
       {'file':'TWHydra_contall_2pcal.image','range':[-0.01,0.05]}])

The difference is small: 

(TWHydra_contall_1pcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean           Rms            Std dev        Minimum        Maximum
2511      5.658062e-01   1.662741e-02   2.253310e-04   3.828434e-03   3.822558e-03   -1.183686e-02  1.252072e-02


(TWHydra_contall_2pcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean           Rms            Std dev        Minimum        Maximum
2511      6.660774e-01   1.957409e-02   2.652638e-04   3.735169e-03   3.726479e-03   -1.143428e-02  1.235826e-02
Residual phase solutions from self_ap.cal.
Amplitude solutions from self_ap.cal.

This small improvement suggests that we have gone as far as we can with phase self-cal. Now we will apply the best phase solution table on-the-fly while solving for the amplitude self-cal solutions. Since amplitude changes more slowly and is less constrained than phase, we use a longer solint='inf', as long as combine=' ' this means to do one solution per scan. 

# In CASA
gaincal(vis='TWHydra_cont.ms',caltable='self_ap.cal',
        solint='inf',combine='',gaintype='T',
        refant='DV06',spw='',minblperant=4,
        gaintable=['self_2.pcal'],
        calmode='ap',minsnr=2)

It's useful to look at the residual phase in the amplitude calibration table, i.e. since we applied our best phase-only solution while solving for this new amplitude and phase (calmode='ap') solution, the phase should already be corrected. 

# In CASA
plotcal(caltable='self_ap.cal',xaxis='time',yaxis='phase',
        spw='',field='',antenna='1~8',iteration='antenna',
        subplot=421,plotrange=[0,0,-1,1],figfile='self_ap_phase.png')

The residuals are very small as expected. (If you see large residuals here, it means that the phase-only solution is suspect, and you were mostly just moving noise around.) Now look at the amplitude solutions. 

# In CASA
plotcal(caltable='self_ap.cal',xaxis='time',yaxis='amp',
        spw='',field='',antenna='1~8',iteration='antenna',
        subplot=421,plotrange=[0,0,0.6,1.4],figfile='self_ap_amp.png')

For the most part these look very good except one rather large deviation on DV06. Apply both the final phase and the amplitude calibration tables.

Post self-calibration amplitude vs. time
Post self-calibration amplitude vs. uv-distance
# In CASA
applycal(vis='TWHydra_cont.ms',field='',gaintable=['self_2.pcal','self_ap.cal'],calwt=F)

Plot the corrected data as a function of both time and uv-distance. It is always a good idea to check these after an amplitude self-calbration to be sure that it worked properly. Occasionally, a few spurious bits of data will get blown up by the amplitude self-cal if it is not well-constrained. 

# In CASA
plotms(vis='TWHydra_cont.ms',xaxis='time',yaxis='amp',
       avgchannel='38',ydatacolumn='corrected',coloraxis='spw',
       plotfile='selfcal_time.png')

The time where we saw a large amplitude correction on DV06 looks consistent here so the large correction we saw in the solutions corrected a real issue. 

# In CASA
plotms(vis='TWHydra_cont.ms',xaxis='uvdist',yaxis='amp',
       avgchannel='38',ydatacolumn='corrected',coloraxis='spw',
       plotfile='selfcal_uvdist.png')

You will see that this plot looks much better than before.

Clean residuals after final amplitude and phase self-cal.

Now make the final continuum image. 

# In CASA
os.system('rm -rf TWHydra_contall_apcal.*')
clean(vis='TWHydra_cont.ms',imagename='TWHydra_contall_apcal',
      mode='mfs',calready=T,imagermode='csclean',
      imsize=100,cell=['0.3arcsec'],spw='',
      weighting='briggs',robust=0.5, 
      mask='TWHydra_contall_2pcal.mask',
      interactive=T,threshold='1mJy',niter=10000)

We know we need to clean at least 200 iterations, so set this in the gui and hit green arrow. After this notice how much better the residual looks. Clean another 100 iterations. The "circular pattern", and 4 bright spots at each corner you see now in the residual is likely due to non-closing amplitude errors that cannot be removed with an antenna based self-calibration, so this is as good as it gets with pure self-cal. Hit the Red X to stop clean. 

Clean used 300 iterations to approach a threshhold of 0.00105507
0.0061808 Jy <- cleaned in this cycle for  model 0
Final maximum residual = 0.00242926
Model 0: max, min residuals = 0.00238278, -0.00242926 clean flux 1.48201

Now repeat the display, blinking, and image statistics described above for the 2nd phase only and the new amplitude & phase self-cal'ed images. 

# In CASA
imview(raster=[{'file':'TWHydra_contall_2pcal.image','range':[-0.01,0.05]},
       {'file':'TWHydra_contall_apcal.image','range':[-0.01,0.05]}])

The amplitude and phase self cal'ed image is notably better than phase-only. Quantitatively: 

(TWHydra_contall_2pcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean           Rms            Std dev        Minimum        Maximum
2621      8.885072e-01   2.611066e-02   3.389955e-04   3.717286e-03   3.702503e-03   -1.089106e-02  1.235826e-02


(TWHydra_contall_apcal.image) image/raster
Stokes=I  Frequency=3.50845e+11Hz  Velocity=0km/s  Frame=TOPO  Doppler=RADIO  BrightnessUnit=Jy/beam  BeamArea=34.0285
Npts      Sum            Flux (Jy)      Mean           Rms            Std dev        Minimum        Maximum
2621      2.142941e-01   6.297484e-03   8.176043e-05   1.660546e-03   1.658849e-03   -5.046569e-03  5.200869e-03

and indeed, we've gained another factor of 2.2 in sensitivity.

Overall with this self-calibration we've improved the noise by a factor of 5.7! The peak flux density also increased from 0.96 to 1.0 Jy, for a change in the S/N from 101 initially to just over 600 after. So well worth the trouble.

It is notable that we are still a factor of more than 5 from the Apriori noise estimate of 0.3 mJy/beam. This is likely due to residual sources of non-closing (not antenna based) errors in the data like baseline errors and baseline-dependent bandpass errors. Additionally, polarization errors can contribute (which we did not calibrate at all), as well as imperfect uv-coverage with only 8 antennas. A dynamic range of 600 is actually quite good in the submillimeter.

Apply self-calibration to Full Dataset and Split Line Data

Now apply the self-calibration derived from the continuum emission to the full unchannel-averaged data. 

# In CASA
applycal(vis='TWHydra_corrected.ms',gaintable=['self_2.pcal','self_ap.cal'],calwt=F)

Split off the CO(3-2) spectral line data in spw=2 and the HCO+(4-3) data in spw=0 for further processing. 

# In CASA
os.system('rm -rf TWHydra_CO3_2.ms*')
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_CO3_2.ms',datacolumn='corrected',spw='2')

# In CASA
os.system('rm -rf TWHydra_HCOplus.ms*')
split(vis='TWHydra_corrected.ms',outputvis='TWHydra_HCOplus.ms',datacolumn='corrected',spw='0')

Continuum Subtraction

Next we need to subtract the continuum emission from the spectral line data. It is best to do this in the uv-plane. The first step is to make a channel plot of the spectral data to get the line-free channel ranges for both spectral lines. Since the line emission is simple (one spectral line per spw) and there are more than enough line-free channels we do not need to be very precise just stay well away from the line and whenever possible it is best to chose line-free channel ranges on either side of the line for the best possible baseline subtraction.

CO(3-2) line in channel space.
HCO+(4-3) in channel space.

Note: after the split above, the original spw ids will be relabeled as 0 in each file because there is only one spectral window in each file. 

# In CASA
plotms(vis='TWHydra_CO3_2.ms',spw='0',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T,coloraxis='spw',plotfile='CO3_2_channel.png')

It is a good idea to stay a little away from the very edges of the baseband which can be noisy. The channel ranges 20~2000 and 2400~3800 look like good choices for CO(3-2). 

# In CASA
plotms(vis='TWHydra_HCOplus.ms',spw='0',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T,coloraxis='spw',plotfile='HCOp4_3_channel.png')

The channel ranges 20~1500 and 1900~3800 look like good choices for HCO+(4-3).

Next use the identified line-free channel ranges and uvcontsub2 to subtract the continuum model from the line emission. 

# In CASA
uvcontsub2(vis='TWHydra_CO3_2.ms',fitorder=1,fitspw='0:20~2000,0:2400~3800')

# In CASA
uvcontsub2(vis='TWHydra_HCOplus.ms',fitorder=1,fitspw='0:20~1500,0:1900~3800')

Note: If you do not have line-free channels on both sides of the line emission, it is safer to use fitorder=0, or the fit may diverge.

Spectral Line Imaging

CO(3-2) spectrum in LSRK velocity space.
HCO+(4-3) spectrum in LSRK velocity space.

Now we are ready to make cubes of the line emission. The imaging parameters are similar to the continuum except for those dealing with the spectral setup: mode, start, width, nchan, restfreq, and outframe parameters. When making spectral images you have three choices for the mode parameter: channel, velocity, and frequency. Data are taken using constant frequency channels. For spectral line analysis it's often more useful to have constant velocity channels, and this is also the best way to make images of multiple lines with the exact same channelization for later comparison. For mode='velocity', the desired start and width also need to be given in velocity units for the desired output frame.

It is important to note that ALMA does not do on-line Doppler Tracking and the native frame of the data is TOPO. If you do not specify outframe the output cube will also be in TOPO, which is not very useful for spectral line work. The Doppler Shift is taken out during the regridding to the desired outframe in clean or alternatively it can be done separately by the cvel task which would need to be run before clean.

To see what velocity parameters you want to set in clean it is useful to make a plot in plotms in the desired output frame. To make it easier to see, we set a plotrange since we know the LSRK velocity of TW Hya is about 2.88 km/s. Note these plots take a little longer because of the frame shift. 

# In CASA
plotms(vis='TWHydra_CO3_2.ms.contsub',xaxis='velocity',yaxis='amp',
       avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',
       restfreq='345.79599GHz',plotrange=[-20,23,0,0],plotfile='CO3_2_vel.png')

# In CASA
plotms(vis='TWHydra_HCOplus.ms.contsub',xaxis='velocity',yaxis='amp',
      avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',
      restfreq='356.7342GHz',plotrange=[-20,23,0,0],plotfile='HCOp4_3_vel.png')

From these plots it looks like imaging the cubes from about -5 to +8 km/s will encompass the line, but still provide several line-free channels on either side so the noise level can be estimated. The channel width is 122 kHz, which at 345.79599 GHz is 0.106 km/s. Recall that the spectral resolution is a factor of two poorer. We will use a velocity channel width of 0.12 km/s for a little oversampling, but the spectral resolution remains 0.2 km/s.

Interactive Viewer window with CO(3-2) clean mask, before cleaning starts.
Interactive Viewer window for CO(3-2), after 500 iterations, more cleaning is needed.

As before, it is very important that you make a clean mask. There are many ways to do this ranging from the complicated to simple. For this example we will make a single clean mask that encompasses the line emission in every channel and apply it to all channels. This is much better than no clean mask, though not quite as good as making an individual mask for each channel. For some sources, one can use the continuum mask to clean the line, however TW Hya is a bit curious in that the continuum emission is not as extended as the line emission (Hughes et al. 2008). 

# In CASA
os.system('rm -rf TWHydra_CO3_2line.*')
clean(vis='TWHydra_CO3_2.ms.contsub',imagename='TWHydra_CO3_2line',
      imagermode='csclean',calready=T,spw='',
      imsize=100,cell=['0.3arcsec'],
      mode='velocity',start='-4km/s',width='0.12km/s',nchan=118,
      restfreq='345.79599GHz',outframe='LSRK',
      weighting='briggs',robust=0.5,
      mask='',
      interactive=T,threshold='1mJy',niter=100000)

It will take a little while to grid the cube, when the interactive viewer opens, use the tapedeck to cycle through the channels. You should see a progress of line emission in channels 40 to 70 or so.

  • Go to channel 57 (about line center)
  • Select the polygon tool with the left mouse button
  • Select the "All Channels" toggle in the green area
  • Draw a polygon around the channel 57 emission, double click inside to activate.
  • Now go to channel 60, this is the channel with the most extent to the SE, draw another polygon that encompasses this emission if necessary,double click inside to activate.
  • Now go to channel 55, this is the channel with the most extent to the NW, draw another polygon that encompasses this emission if necessary,double click inside to activate.
  • Now double check the mask by cycling through channels with real line emission to be sure that your clean box encompasses the emission, but does not extend far beyond it. Again the "erase" toggle can come in handy if things go awry or if you want to just delete a small portion of an existing mask (just draw a polygon around what you want to erase and double click inside), be sure to switch back to "add" before going on.

The figure shows what the final clean box should look like. Now clean using the green circle arrow. You can watch the progress in the logger. When 100 iterations are done, the viewer will show the residual map for each channel. Cycle through the channels and see whether you are still happy with the clean box in each signal channel. The "erase button" can help you fix mistakes. If necessary adjust. It is often useful to adjust the colorscale with the "plus" symbol icon. Keep going until the residual emission looks like the surrounding noise. If it seems to be going slowly, increase the iterations to 300 for a round or two. I stopped cleaning (Red x) when the logger showed (about 900 iterations). 

12.9846 Jy <- cleaned in this cycle for  model 0
Final maximum residual = 0.157312
Model 0: max, min residuals = 0.157312, -0.0871295 clean flux 279.304

NOTE: The residuals are still a bit high in the strongest emission channels after it is well cleaned in the weaker channels. This is an indication that the sensitivity is "dynamic range" limited. In other words even with more intrinsic sensitivity per integration (i.e. better system temperature), the dynamic range limit would remain the same unless you get more uv-coverage. The fundamental dynamic range limitations of ALMA will be considerably better with 16 antennas in Early Science and MUCH better with the full 50-antenna array.

Now image and clean the HCO+(4-3) line using what you learned above. 

# In CASA
os.system('rm -rf TWHydra_HCOplusline.*')
clean(vis='TWHydra_HCOplus.ms.contsub',imagename='TWHydra_HCOplusline',
      imagermode='csclean',calready=T,spw='',
      imsize=100,cell=['0.3arcsec'],
      mode='velocity',start='-4km/s',nchan=118,width='0.12km/s',
      restfreq='356.7342GHz',outframe='LSRK',
      weighting='briggs',robust=0.5,
      mask='',
      interactive=T,threshold='1mJy',niter=100000)

I stopped cleaning when: 

4.79833 Jy <- cleaned in this cycle for  model 0
Final maximum residual = 0.0977348
Model 0: max, min residuals = 0.0977348, -0.060484 clean flux 192.066

Image Analysis

Moment Maps

Next we will make moment maps for the image cubes.

To determine channel ranges it is useful to open the cube in viewer as a raster image. Using the statistics tool the rms in a line-free channel is about 25 mJy/beam, however in one of the peak line channels it is as high as 50 mJy/beam. This difference is again because of the dynamic range limitations of the data (as described above for the continuum). Now open the cube again but as a contour map and set the Unit Contour Level to 0.025 (this is in Jy/beam) and in Relative Contour levels set [3]. i.e. this is showing 3sigma contours. Now find out where the CO(3-2) emission drops below the 3sigma level.

For moment 0 (integrated intensity) maps you do not typically want to set a flux threshold because this will tend to noise bias your integrated intensity. 

# In CASA
immoments(imagename='TWHydra_CO3_2line.image',moments=[0],
          outfile='TWHydra_CO3_2line.image.mom0',
          chans='40~76')

For higher order moments it is very important to set a conservative flux threshold. Typically something like 6sigma, using sigma from peak line channel works well. 

# In CASA
immoments(imagename='TWHydra_CO3_2line.image',moments=[1,2],
          outfile='TWHydra_CO3_2line.image.mom',
          chans='40~76',includepix=[0.3,100])

Display all three moment maps in the viewer, overlaid with continuum contours. 

# In CASA
imview( raster=[ {'file':'TWHydra_CO3_2line.image.mom0'},
                 {'file':'TWHydra_CO3_2line.image.mom.weighted_coord'},
                 {'file':'TWHydra_CO3_2line.image.mom.weighted_dispersion_coord'} ], 
        contour={'file':'TWHydra_contall_apcal.image', 'base':0, 'unit':0.0025, 'levels':[3,100]} )

To see all three raster images simultaneously, open the viewer's Panel Display Options (the "P-wrench" icon in the viewer tool bar), click "Number of panels", and change "Number of panels in x" to 3. Then select the "Blink" toggle (radio button) from the Animator panel.

ALMA CO(3-2) moment maps, with white continuum contours at 3 and 100 sigma. From left to right: integrated intensity, intensity weighted velocity field, intensity weighted velocity dispersion are shown.

Repeat for HCO+(4-3) which has a bit narrower velocity extent than the CO(3-2), but similar rms noise. 

# In CASA
immoments(imagename='TWHydra_HCOplusline.image',moments=[0],
          outfile='TWHydra_HCOplusline.image.mom0',
          chans='43~74')

# In CASA
immoments(imagename='TWHydra_HCOplusline.image',moments=[1,2],
          outfile='TWHydra_HCOplusline.image.mom',
          chans='43~74',includepix=[0.3,100])

Display all three moment maps in the viewer using the same method described above for the CO(3-2) maps.

ALMA HCO+(4-3) moment maps, with white continuum contours at 3 and 100 sigma. From left to right: integrated intensity, intensity weighted velocity field, intensity weighted velocity dispersion are shown.

Channel Maps

Using the Viewer you can make channel map figures. Start the Viewer and then open the CO(3-2) cube as a raster image and then the continuum as a contour image. Then we use the "wrench" icon and "P wrench" icons to set up the channel images, contour levels etc.

Channel maps of the CO(3-2) emission with white continuum contours at 3 and 100 sigma.

The setup windows looked like this

Repeat for HCO+(4-3)

Channel maps of the HCO+(4-3) emission with white continuum contours at 3 and 100 sigma. The color intensity scale is the same as that used for CO(3-2).

Exporting Fits Images

If you want to analyze the data using another software package it is easy to convert from CASA format to FITS. 

# In CASA
os.system('rm -rf TWHydra_CO3_2line.image.fits')
exportfits(imagename='TWHydra_CO3_2line.image',fitsimage='TWHydra_CO3_2line.image.fits')

Although "fits format" is supposed to be a standard, in fact most packages expect slightly different things from a fits image. If you are having difficulty, try setting velocity=T and/or dropstokes=T. 

# In CASA

os.system('rm -rf TWHydra_HCOplusline.image.fits')
exportfits(imagename='TWHydra_HCOplusline.image',fitsimage='TWHydra_HCOplusline.image.fits')

# In CASA
os.system('rm -rf TWHydra_CO3_2line.image.mom0.fits')
exportfits(imagename='TWHydra_CO3_2line.image.mom0',
         fitsimage='TWHydra_CO3_2line.image.mom0.fits')

# In CASA
os.system('rm -rf TWHydra_CO3_2line.image.mom.weighted_coord.fits')
exportfits(imagename='TWHydra_CO3_2line.image.mom.weighted_coord',
         fitsimage='TWHydra_CO3_2line.image.mom.weighted_coord.fits')

# In CASA
os.system('rm -rf TWHydra_CO3_2line.image.mom.weighted_dispersion_coord.fits')
exportfits(imagename='TWHydra_CO3_2line.image.mom.weighted_dispersion_coord',
         fitsimage='TWHydra_CO3_2line.image.mom.weighted_dispersion_coord.fits')

# In CASA
os.system('rm -rf TWHydra_HCOplusline.image.mom0.fits')
exportfits(imagename='TWHydra_HCOplusline.image.mom0',
         fitsimage='TWHydra_HCOplusline.image.mom0.fits')

# In CASA
os.system('rm -rf TWHydra_HCOplusline.image.mom.weighted_coord.fits')
exportfits(imagename='TWHydra_HCOplusline.image.mom.weighted_coord',
         fitsimage='TWHydra_HCOplusline.image.mom.weighted_coord.fits')

# In CASA
os.system('rm -rf TWHydra_HCOplusline.image.mom.weighted_dispersion_coord.fits')
exportfits(imagename='TWHydra_HCOplusline.image.mom.weighted_dispersion_coord',
         fitsimage='TWHydra_HCOplusline.image.mom.weighted_dispersion_coord.fits')

# In CASA
os.system('rm -rf TWHydra_contall_apcal.image.fits')
exportfits(imagename='TWHydra_contall_apcal.image',fitsimage='TWHydra_contall_apcal.image.fits')

Last checked on CASA Version 3.3.0.

Retrieved from "https://casaguides.nrao.edu/index.php?title=TWHydraBand7_Imaging_for_CASA_3.3&oldid=10596"