############################################################################################## #################################### General Information ##################################### ############################################################################################## ##### The following script has been designed to image the Band 4 observations of SDP81 obtained during the long baseline campaign: https://casaguides.nrao.edu/index.php/ALMA2014_LBC_SVDATA#SDP.81_-_Lensed_Galaxy # The Observed line is CO(5-4) at z=3.042, redshifted to freq=142.5700GHz ##### This script has been validated for execution in CASA 6.4. It makes use of the analysisutils environment (aU), which can be downloaded separately from CASA. ##### This script is meant to be run interactively ##### At any time, if you think you modified SDP.81_Band4_continuum.ms, SDP.81_Band4.ms or bandpass.ms incorrectly, you can erase them from your work area and copy them 'fresh' from the common repository ############################################################################################## ################################## Startup and Orientation ################################### ############################################################################################## ##### start CASA in your terminal: type casa ##### importing packages and aliases import re import os import analysisUtils as aU es = aU.stuffForScienceDataReduction() ##### listobs: prints the metadata associated to the measurement set to the logger listobs('SDP.81_Band4.ms') ##### we recommend recording the listobs output to a text file, with the following command: os.system('rm -fr SDP.81_Band4.ms.listobs') listobs(vis = 'SDP.81_Band4.ms', listfile = 'SDP.81_Band4.ms.listobs') ##### Check the Fourier plane coverage plotms(vis='SDP.81_Band4.ms', xaxis='u', yaxis='v',avgchannel='10000', avgspw=False, avgtime='1e9', avgscan=False,coloraxis="observation") ############################################################################################## ############################## Imaging the Bandpass Calibrator ############################### ############################################################################################## ##### let's start with imaging a bright and point-like source ##### May need multiple rounds of cleaning needed before see a clear change in residuals os.system('rm -rf bandpass_natural.*') tclean(vis='bandpass.ms', imagename='bandpass_natural', field='0', spw='', specmode='mfs',deconvolver='hogbom',gridder='standard', imsize=[512,512], cell=['0.005arcsec'], weighting='natural', threshold='0mJy', niter=10000,interactive=True) imview('bandpass_natural.image') os.system('rm -rf bandpass_robust.*') tclean(vis='bandpass.ms', imagename='bandpass_robust', field='0', spw='', specmode='mfs', deconvolver='hogbom',gridder='standard', imsize=[512,512], cell=['0.005arcsec'], weighting='briggs',robust=0.0, threshold='0mJy', niter=10000, interactive=True) imview('bandpass_robust.image') ##### Different weighting schemes result in different synthesized beam size and often sensitivity ##### You can play with the weigthing scheme - uniform, natural or briggs (robust), and for briggs weighting, the robust parameter ##### What happens if one chooses a very large pixel size? os.system('rm -rf bandpass_bigpix.*') tclean(vis='bandpass.ms', imagename='bandpass_bigpix', field='0', spw='', specmode='mfs', deconvolver='hogbom',gridder='standard', imsize=[128,128], cell=['0.05arcsec'], weighting='briggs',robust=-1, threshold='0mJy', niter=10000,interactive=True) imview('bandpass_bigpix.image') ############################################################################################## ################################# Imaging SDP.81 Continuum ################################### ############################################################################################## ##### SDP81 has a complex structure, and will need tuning of several parameters to achieve a statisfactory imaging quality ##### We image the continuum data, using briggs (robust) weighting (with robust=1), and multiscale cleaning. The mask is selected during interactive mode. ##### based on the LBC campaign SV data CASA guide, here is a summary description of the use of multiscale #####'''"Extended emission at high angular resolution can best be imaged with multiscale: Cleaning diffuse emission on size scales significantly larger than the beam with the standard delta-function deconvolution method results in the "clean instability": giving the diffuse emission a "pointilated" or "cotton-candy" like morphology. This can be ameliorated by using a range of scales to do the cleaning. In CASA, this technique is accessed through the tclean multiscale parameter. The scales are given as multipliers of the cell (pixel) size. It is essential to always use 0 for the first scale as this corresponds to the normal clean beam, multiscale=[0,5,15] is usually a good starting point for experimentation. For more details see: Cornwell 2008 (IEEE Journal of Sig Proc., 2, 793).''' os.system('rm -rf SDP.81.continuum_multiscale.*') tclean(vis='SDP.81_Band4_continuum.ms', imagename='SDP.81.continuum_multiscale', spw='',field='SDP*', specmode='mfs',gridder='standard',deconvolver='multiscale', imsize=1500,cell='0.01arcsec', scales=[0,5,15,45], interactive=True, mask = '', weighting='briggs',robust=1.0, niter=10000,threshold='0.02mJy') # sigma = 0.008mJy/beam ##### view the resulting clean image imview('SDP.81.continuum_multiscale.image') ##### It is difficult to define a mask. The clean image only presents emission on very small scales, and there may be diffuse emission on different scales. We can try to play with the robust parameter to try to 'pick up' more emission, as well as with the multiscale parameter. ##### Since some emission is still resolved out at this angular resolution, we image the target tapering the uv-data at long baselines to emphasize and recover more of the extended emission. os.system('rm -rf SDP.81.continuum_smooth.*') tclean(vis='SDP.81_Band4_continuum.ms', imagename='SDP.81.continuum_smooth', spw='',field='SDP*', specmode='mfs',gridder='standard',deconvolver='multiscale', imsize=1500,cell='0.01arcsec', scales=[0,5,15,45], interactive=True, mask='', weighting='briggs',robust=1.0, uvtaper=['1000klambda'], niter=10000,threshold='0.025mJy') ##### view the resulting clean image imview('SDP.81.continuum_smooth.image') ##### view the resulting clean image as a countour plot imview({'file':'SDP.81.continuum_smooth.image','levels':[0.2,0.4,0.6,0.8],'unit':0.0002'}) ##### modify the contour levels in Data -> Adjust data display, basic settings ##### Export to fits file: exportfits(imagename='SDP.81.continuum_smooth.image', fitsimage='SDP.81_B4_cont_smooth.fits') ############################################################################################## ################################## Imaging SDP.81 CO Line #################################### ############################################################################################## ##### The Spectral line is CO(5-4) at z=3.042, redshifted to freq=142.5700GHz ##### Next, we split out the continuum and the line spectral windows using the CASA task 'split'. ##### We then subtract the continuum in the line spws after identifying the channels with line emission. ##### The first two steps in this section are commented out, and the resulting intermediate file has been placed ##### in your imaging directory. The 'split' and 'uvcontsub' tasks can take some time, depending on the size of ##### your ms and the specs of your workstation. ##### Define spectral windows for the continuum and line emission ##### TDM spectral windows (Time Domain Mode) have wide channels and are dedicated to continuum emission ##### FDM spectral windows (Frequency Domain Mode) have narrow channels and are targeting the line emission ##### one can find this information from the listobs file # spw_cont = '0~2,4~6,8~10,12~14,16~18,20~22,24~26,28~30,32~34,36~38,40~42,44~46' # spw_line = '3,7,11,15,19,23,27,31,35,39,43,47' ##### Split continuum and spectral line observations # os.system('rm -rf SDP.81_Band4_COline.ms') # split(vis='SDP.81_Band4.ms',outputvis='SDP.81_Band4_COline.ms', # spw=spw_line,datacolumn='data') ##### Subtracting the continuum ##### Plot the data to determine uvcontsub parameters # Note point shape/size default is smaller than in presentation, can adjust through 'Display' tab # plotms('SDP.81_Band4_COline.ms',yaxis='amp',xaxis='channel',avgtime='1e8', # coloraxis='spw',restfreq='142.5700GHz',freqframe='LSRK',transform=True,avgantenna=True,avgscan=True) # os.system('rm -rf SDP.81_Band4_COline.ms.contsub') # uvcontsub(vis='SDP.81_Band4_COline.ms', # fitorder=1,fitspw='0~11:5~45;170~187') ##### Take a look at the resulting continuum-subtracted data # Note point shape/size default is smaller than in presentation, can adjust through 'Display' tab plotms('SDP.81_Band4_COline.ms.contsub',yaxis='amp',xaxis='channel',avgtime='1e8', coloraxis='spw',restfreq='142.5700GHz',freqframe='LSRK',transform=True,avgantenna=True,avgscan=True) ##### Imaging target - robust=1 weighting and uvtapering of longest baselines: ##### ----------------------------------------------------------------------- ##### CO emission detected between +200 to -400 km/s, hence only these channels ##### have a cleaning box, defined during interactive cleaning os.system('rm -rf SDP.81.Band4.CO_smooth.*') tclean(vis='SDP.81_Band4_COline.ms.contsub', imagename='SDP.81.Band4.CO_smooth', mask='', specmode='cube',gridder='standard',deconvolver='multiscale', imsize=672, cell='0.02arcsec', start='-520km/s',width='21km/s',nchan=45, outframe='LSRK',restfreq='142.5700GHz', scales=[0,5,15,45],#negcomponent=10, interactive=True, restoringbeam='common', weighting='briggsbwtaper',robust=1.0, uvtaper=['1000klambda'], perchanweightdensity = True, niter=10000,threshold='0.52mJy') # imview('SDP.81.Band4.CO_smooth.image') immoments(imagename='SDP.81.Band4.CO_smooth.image', chans='10~36',moments=[0],includepix=[2*0.20e-3,100], # -300 to +240 km/s outfile='SDP.81.Band4.CO_smooth.image.mom0_2sigma') immoments(imagename='SDP.81.Band4.CO_smooth.image', chans='10~36',moments=[1],includepix=[4*0.20e-3,100], # -300 to +240 km/s outfile='SDP.81.Band4.CO_smooth.image.mom1_4sigma') imview('SDP.81.Band4.CO_smooth.image.mom0_2sigma') imview('SDP.81.Band4.CO_smooth.image.mom1_4sigma') exportfits(imagename='SDP.81.Band4.CO_smooth.image', fitsimage='SDP.81.Band4.CO_smooth_z3.042.fits') exportfits(imagename='SDP.81.Band4.CO_smooth.image.mom0_2sigma', fitsimage='SDP.81.Band4.CO_smooth.mom0_z3.042.fits') exportfits(imagename='SDP.81.Band4.CO_smooth.image.mom1_4sigma', fitsimage='SDP.81.Band4.CO_smooth.mom1_z3.042.fits')