##############################################################################
# Data reduction script for VY CMa Band 9 658-GHz 12-m data: Part 2: Imaging
# Tested using CASA Version 4.2.2
#
# Datasets used:
#
# uid___A002_X6d34fd_Xb3d - observed Aug 18 2013 with 18 antennas
# Flux scale calibrator J0522-364
# Bandpass calibrator J0522-364
# Phase reference source J0648-3044
# Target Vy CMa
#
##############################################################################
"""
See accompanying README file for details of the necessary input files
and comments on the data.
Before running this script you should have run
VYCMa_Band9_Calibration.py or obtained the calibrated
target dataset VYCMa_658.ms. This script assumes that VYCMa_658.ms
has had all instrumental and reference source calibration applied, and
has been adjusted to the LSRK frame, as at the end of
VYCMa_Band9_Calibration.py
"""
#----------------------------------------------------------------------------------
#----- Optional Steps -------------------------------------------------------------
#----------------------------------------------------------------------------------
#----- Set this option to False if you do not want to do interactive clean.
#----- Note that better results will be obtained if you do make masks interatively,
#----- expanding the area as fainter emission emerges
cleaninteractive=T
#----- Set this option to true if you want to make diagnostic plots along the way.
#----- Note that this will slow down the reduction.
#----- Default is no plots (you can make them later anyway).
calplots=T
#----------------------------------------------------------------------------------
#----- Some setup steps -----------------------------------------------------------
#----------------------------------------------------------------------------------
version = casalog.version()
print "You are using " + version
if (version < '4.2.1'):
print "\033[91m YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
print "\033[91m PLEASE UPDATE IT BEFORE PROCEEDING."
else:
print "Your version of CASA is appropriate for this guide."
# Repeat this if you are using copy and paste, and you re-start CASA
refantenna='DV15'
thesteps=[]
step_title = {0: 'List the data set and plot the visibility spectrum',
1: 'Make first image of maser peak',
2: 'Self-calibrate cycle 1, phase-only, re-image',
3: 'Self-calibrate cycle 2, amp & phase, re-image',
4: 'Self-calibrate cycle 3, phase, better model, re-image',
5: 'Self-calibrate cycle 4, amp & phase, re-image',
6: 'Self-calibrate cycle 4, amp & phase, shorter solint, re-image',
7: 'Self-calibrate cycle 5, amp & phase, shorter solint, re-image',
8: 'Make low-res cubes to identify continuum \n(this step can be omitted if you trust our channel selection)',
9: 'Image continuum',
10: 'Subtract continuum',
11: 'Image a line in spw 0',
}
# The Python variable 'mysteps' will control which steps
# are executed when you start the script using
# execfile('scriptForCalibration.py')
# e.g. setting
# mysteps = [2,3,4]# before starting the script will make the script execute
# only steps 2, 3, and 4
# Setting mysteps = [] will make it execute all steps.
try:
print 'List of steps to be executed ...', mysteps
thesteps = mysteps
except:
print 'global variable mysteps not set.'
if (thesteps==[]):
thesteps = range(0,len(step_title))
print 'Executing all steps: ', thesteps
######################################################
# Do listobs and write info to a text file; plot spectrum and brightest channel
mystep = 0
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm VYCMa_658_listobs.txt')
listobs(vis='VYCMa_658.ms', listfile='VYCMa_658_listobs.txt', verbose=True)
print '<< Printed listobs output to VYCMa_658_listobs.txt'
os.system('rm -rf VYCMa_658_spw_1st-vis-spectrum.png')
plotms(vis='VYCMa_658.ms',
xaxis='frequency', yaxis='amp',
selectdata=T,
avgtime='1e8',avgscan=T,
coloraxis='baseline',
showgui=F,
plotfile='VYCMa_658_spw_1st-vis-spectrum.png')
# You will notice the very bright maser line in spw 0:1966
################################
# Last scan is NOISY?
# ANY MISBEHAVING ANTENNAS?
################################
os.system('rm -rf VYCMa_658_0-1966_1st_amp.png')
plotms(vis='VYCMa_658.ms',
xaxis='time', yaxis='amp',
selectdata=T,
antenna='DV15&*',
spw='0:1966',
coloraxis='baseline',
showgui=F,
plotfile='VYCMa_658_0-1966_1st_amp.png')
os.system('rm -rf VYCMa_658_0-1966_1st_phase.png')
plotms(vis='VYCMa_658.ms',
xaxis='time', yaxis='phase',
selectdata=T,
antenna='DV15&*',
spw='0:1966',
coloraxis='baseline',
showgui=F,
plotfile='VYCMa_658_0-1966_1st_phase.png')
######################################################
# Make first image of maser peak, omitting dubious data
# *Scan 23 is noisy, also DA41 especially in 19 and 24
mystep = 2
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf VYCMa_658_0_1966.clean*')
clean(vis = 'VYCMa_658.ms',
imagename='VYCMa_658_0_1966.clean',
field='Vy CMa',
spw='0:1966',
scan='7~22,24~73',
antenna='!DA41',
mode='channel',start=1966,nchan=1,
cell='0.005arcsec',npercycle=10,imsize=1024,
psfmode='hogbom',
mask='ellipse[[523pix,563pix],[62pix,39pix],90deg]',
interactive=cleaninteractive,
niter=50)
# Get image statistics
imin='VYCMa_658_0_1966.clean.image'
noise=imstat(imagename=imin,
box='500,100,900,400')['rms'][0]
peak=imstat(imagename=imin,
box='400,400,699,699')['max'][0]
print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)
#rms 20.029, peak 342.312, snr 17
######################################################
# Self-calibrate, phase-only first, using all data
mystep = 3
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf cal-VYCMa_658_1966.ph1')
gaincal(vis = 'VYCMa_658.ms',
field='Vy CMa',
refant=refantenna,
minsnr=2,
caltable='cal-VYCMa_658_1966.ph1',
antenna='!DA41',
spw='0:1966',
calmode='p',
solint='120s')
if (calplots):
os.system('rm -rf cal-VYCMa_658_1966.ph1*.png')
for ant in range(3):
plotcal(caltable='cal-VYCMa_658_1966.ph1',
xaxis = 'time', yaxis = 'phase',
subplot=421,
fontsize=7.0,
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration='antenna',
figfile='cal-VYCMa_658_1966.ph1_'+str(ant*8)+'-'+str(ant*8+7)+'.png')
applycal(vis = 'VYCMa_658.ms',
field='Vy CMa',
spw='0',
gaintable='cal-VYCMa_658_1966.ph1',
calwt = F,
applymode='calonly',
flagbackup = F)
os.system('rm -rf VYCMa_658_0_1966_ph1.clean*')
clean(vis = 'VYCMa_658.ms',
imagename='VYCMa_658_0_1966_ph1.clean',
field='Vy CMa',
spw='0:1966',
scan='7~22,24~73',
antenna='!DA41',
mode='channel',start=1966,nchan=1,
cell='0.005arcsec',npercycle=40,imsize=1024,
psfmode='hogbom',
mask='ellipse[[520pix,564pix],[89pix,54pix],90deg]',
interactive=cleaninteractive,
niter=200)
# Get image statistics
imin='VYCMa_658_0_1966_ph1.clean.image'
noise=imstat(imagename=imin,
box='500,100,900,400')['rms'][0]
peak=imstat(imagename=imin,
box='400,400,699,699')['max'][0]
print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)
#rms 4.220, peak 429.636, snr 101
######################################################
# Self-calibrate cycle 2, amplitude and phase,
mystep = 4
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf VYCMa_658_1966.ap1')
gaincal(vis = 'VYCMa_658.ms',
field='Vy CMa',
antenna='!DA41',
refant=refantenna,
caltable='cal-VYCMa_658_1966.ap1',
spw='0:1966',
calmode='ap',
solint='300s',
minsnr=2,
solnorm=T,
gaintable='cal-VYCMa_658_1966.ph1')
if (calplots):
os.system('rm -rf cal-VYCMa_658_1966.ap1*.png')
for ant in range(3):
plotcal(caltable='cal-VYCMa_658_1966.ap1',
xaxis = 'time', yaxis = 'phase',
subplot=421,
fontsize=7.0,
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration='antenna',
figfile='cal-VYCMa_1_658_1966.ap1_ph_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')
plotcal(caltable='cal-VYCMa_658_1966.ap1',
xaxis = 'time', yaxis = 'amp',
subplot=421,
fontsize=7.0,
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration='antenna',
figfile='cal-VYCMa_1_658_1966.ap1_amp_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')
applycal(vis = 'VYCMa_658.ms',
field='Vy CMa',
spw='0',
gaintable=['cal-VYCMa_658_1966.ph1','cal-VYCMa_658_1966.ap1'],
calwt = F,
applymode='calonly',
flagbackup = F)
os.system('rm -rf VYCMa_658_0_1966_ap1.clean*')
clean(vis = 'VYCMa_658.ms',
imagename='VYCMa_658_0_1966_ap1.clean',
field='Vy CMa',
spw='0:1966',
scan='7~22,24~73',
antenna='!DA41',
mode='channel',start=1966,nchan=1,
cell='0.005arcsec',npercycle=50,imsize=1024,
psfmode='hogbom',
mask='ellipse[[519pix,567pix],[90pix,57pix],90deg]',
interactive=cleaninteractive,
niter=500)
# Get image statistics
imin='VYCMa_658_0_1966_ap1.clean.image'
noise=imstat(imagename=imin,
box='500,100,900,400')['rms'][0]
peak=imstat(imagename=imin,
box='400,400,699,699')['max'][0]
print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)
#rms 1.296, peak 404.728, snr 312
######################################################
# Restart gaincal with better model, all data
mystep = 5
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf cal-VYCMa_658_1966.ph2')
gaincal(vis = 'VYCMa_658.ms',
field='Vy CMa',
refant=refantenna,
minsnr=2,
caltable='cal-VYCMa_658_1966.ph2',
spw='0:1966',
calmode='p',
solint='60s')
if (calplots):
os.system('rm -rf cal-VYCMa_658_1966.ph2*.png')
for ant in range(3):
plotcal(caltable='cal-VYCMa_658_1966.ph2',
xaxis = 'time', yaxis = 'phase',
subplot=421,
fontsize=7.0,
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration='antenna',
figfile='cal-VYCMa_1_658_1966.ph1_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')
applycal(vis = 'VYCMa_658.ms',
field='Vy CMa',
spw='0',
gaintable='cal-VYCMa_658_1966.ph2',
calwt = F,
applymode='calonly',
flagbackup = F)
#
os.system('rm -rf VYCMa_658_0_1966_ph2.clean*')
clean(vis = 'VYCMa_658.ms',
imagename='VYCMa_658_0_1966_ph2.clean',
field='Vy CMa',
spw='0:1966',
mode='channel',start=1966,nchan=1,
cell='0.005arcsec',npercycle=30,imsize=1024,
psfmode='hogbom',
mask='ellipse[[519pix,567pix],[90pix,57pix],90deg]',
interactive=cleaninteractive,
niter=100)
# Get image statistics
imin='VYCMa_658_0_1966_ph2.clean.image'
noise=imstat(imagename=imin,
box='500,100,900,400')['rms'][0]
peak=imstat(imagename=imin,
box='400,400,699,699')['max'][0]
print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)
# rms 4.918, peak 463.782, snr 94
######################################################
# Self-calibrate cycle 4, ap,
mystep = 6
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf VYCMa_658_1966.ap2')
gaincal(vis = 'VYCMa_658.ms',
field='Vy CMa',
refant=refantenna,
caltable='cal-VYCMa_658_1966.ap2',
spw='0:1966',
gaintable='cal-VYCMa_658_1966.ph2',
calmode='ap',
minsnr=2,
solint='120s')
if (calplots):
os.system('rm -rf cal-VYCMa_658_1966.ap2*.png')
for ant in range(3):
plotcal(caltable='cal-VYCMa_658_1966.ap2',
xaxis = 'time', yaxis = 'amp',
subplot=421,
fontsize=7.0,
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration='antenna',
figfile='cal-VYCMa_1_658_1966.ap2_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')
applycal(vis = 'VYCMa_658.ms',
field='Vy CMa',
spw='0',
gaintable=['cal-VYCMa_658_1966.ph2','cal-VYCMa_658_1966.ap2'],
calwt = F,
applymode='calonly',
flagbackup = F)
#
os.system('rm -rf VYCMa_658_0_1966_ap2.clean*')
clean(vis = 'VYCMa_658.ms',
imagename='VYCMa_658_0_1966_ap2.clean',
field='Vy CMa',
spw='0:1966',
mode='channel',start=1966,nchan=1,
cell='0.005arcsec',npercycle=100,imsize=1024,
psfmode='hogbom',
mask='ellipse[[519pix,567pix],[90pix,57pix],90deg]',
interactive=cleaninteractive,
niter=1000)
# Get image statistics
imin='VYCMa_658_0_1966_ap2.clean.image'
noise=imstat(imagename=imin,
box='500,100,900,400')['rms'][0]
peak=imstat(imagename=imin,
box='400,400,699,699')['max'][0]
print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)
# rms 1.080, peak 456.564, snr 422
######################################################
# Self-calibrate cycle 5, ap, shorter solint
mystep = 7
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf cal-VYCMa_658_1966.ap3')
gaincal(vis = 'VYCMa_658.ms',
field='Vy CMa',
refant=refantenna,
caltable='cal-VYCMa_658_1966.ap3',
spw='0:1966',
gaintable=['cal-VYCMa_658_1966.ph2','cal-VYCMa_658_1966.ap2'],
calmode='ap',
solint='30s')
if (calplots):
os.system('rm -rf cal-VYCMa_658_1966.ap3*.png')
for ant in range(3):
plotcal(caltable='cal-VYCMa_658_1966.ap3',
xaxis = 'time', yaxis = 'amp',
subplot=421,
fontsize=7.0,
antenna=str(ant*8)+'~'+str(ant*8+7),
iteration='antenna',
figfile='cal-VYCMa_1_658_1966.ap3_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png')
applycal(vis = 'VYCMa_658.ms',
field='Vy CMa',
spw='0',
gaintable=['cal-VYCMa_658_1966.ph2','cal-VYCMa_658_1966.ap2','cal-VYCMa_658_1966.ap3'],
calwt = F,
flagbackup = F)
#
os.system('rm -rf VYCMa_658_0_1966_ap3.clean*')
clean(vis = 'VYCMa_658.ms',
imagename='VYCMa_658_0_1966_ap3.clean',
field='Vy CMa',
spw='0:1966',
mode='channel',start=1966,nchan=1,
cell='0.005arcsec',npercycle=100,imsize=1024,
psfmode='hogbom',
mask='ellipse[[519pix,567pix],[90pix,57pix],90deg]',
interactive=cleaninteractive,
niter=1000)
# Get image statistics
imin='VYCMa_658_0_1966_ap3.clean.image'
noise=imstat(imagename=imin,
box='500,100,900,400')['rms'][0]
peak=imstat(imagename=imin,
box='400,400,699,699')['max'][0]
print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise)
#rms 0.888, peak 455.488, snr 512
#####################################
# Image whole cube at lower resolution to identify line-free channels
mystep = 8
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf VYCMa_658_5ap3.clean*')
clean(vis = 'VYCMa_658.ms',
imagename='VYCMa_658_5ap3.clean',
field='Vy CMa',
mode='channel',width=5,
cell='0.01arcsec',imsize=512,usescratch=T,
psfmode='hogbom',
interactive=cleaninteractive,
niter=1000)
# The following is our best guess at a selection
contchans='0:5~94;425~579;1175~1229;1295~1344;2540~2734;3370~3469;3573~3759'
##############################################
# Image continuum
mystep = 9
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf VYCMa_658_contap4.clean*')
clean(vis = 'VYCMa_658.ms',
imagename='VYCMa_658_contap4.clean',
field='Vy CMa',
spw=contchans,
mode='mfs',
multiscale=[0,4,8,12] ,
mask=['ellipse[[436pix,510pix],[128pix,102pix],90deg]','ellipse[[540pix,718pix],[280pix,140pix],20deg]'],
cell='0.005arcsec',imsize=1024,
psfmode='hogbom',interactive=cleaninteractive)
exportfits(imagename='VYCMa_658_contap4.clean.image',
fitsimage='VYCMa_658_contap4.clean.fits')
# PROBLEM #
#####################################################
#### TRY TO SELFCAL ON THE CONTINUUM EMISSION
#####################################################
# PROBLEM #
#####################################################
#### TRY TO SELFCAL ON A WIDER CHANNEL RANGE
#####################################################
# Subtract continuum
#####################################################
#### FIND THE SPECTRAL RANGE FOR CONTINUUM SUBTRACTION
#####################################################
mystep = 10
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
os.system('rm -rf VYCMa_658.ms.contsub')
uvcontsub(vis = 'VYCMa_658.ms',field='V*',fitorder=1,
solint='30s',spw='',fitspw=contchans)
#----------------------------------------------------------------------------------
#----- Image maser and make moment map --------------------------------------------
#----------------------------------------------------------------------------------
#####################################
# Image H2O line spw 0
mystep = 11
if(mystep in thesteps):
casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
print 'Step ', mystep, step_title[mystep]
# High resolution for maser
os.system('rm -rf VYCMa_658.0_H2O.clean*')
clean(vis = 'VYCMa_658.ms.contsub',
imagename='VYCMa_658.0_H2O.clean',
field='Vy CMa',
spw='0',
mode='channel',start=1825,nchan=255,
restfreq='658.0037GHz',
cell='0.005arcsec',imsize=1024,
psfmode='hogbom',
weighting='briggs',robust=0.5,
threshold='100mJy',
mask='ellipse[[519pix,567pix],[90pix,57pix],90deg]',
interactive=cleaninteractive,
usescratch=T,
niter=5000)
# Make moment map
thisimage='VYCMa_658.0_H2O.clean.image'
chanstat=imstat(imagename=thisimage,box='500,100,900,400',chans='100')
rms1= chanstat['rms'][0]
chanstat=imstat(imagename=thisimage,box='500,100,900,400',chans='156')
rms2= chanstat['rms'][0]
rms=0.5*(rms1+rms2)
print 'rms in a bright channel = '+str(rms)
os.system('rm -rf '+thisimage+'.mom0')
immoments(imagename = thisimage,
moments = [0],
axis = 'spectral',
includepix = [rms*3,5000.],
outfile = thisimage+'.mom0')
# Export as fits
exportfits(imagename=thisimage+'.mom0',
fitsimage=thisimage+'.mom0.fits')
exportfits(imagename='VYCMa_658.0_H2O.clean.image',
fitsimage='VYCMa_658.0_H2O.clean.image.fits')
#----------------------------------------------------------------------------------
#----- End of imaging script.
#----------------------------------------------------------------------------------
imview(raster={'file': 'VYCMa_658_0_1966.clean.image','scaling': -1.4,
'range': [-19.6376, 200.725],
'colorwedge':T, 'colormap': 'Rainbow 2', }, out='VYCMa_658_0_1966.clean.png',
zoom=2)
imview(raster={'file': 'VYCMa_658_0_1966_ph1.clean.image','scaling': -1.4,
'range': [-19.6376, 200.725],
'colorwedge':T, 'colormap': 'Rainbow 2', }, out='VYCMa_658_0_1966_ph1.clean.png',
zoom=2)
imview(raster={'file': 'VYCMa_658_0_1966_ph2.clean.image','scaling': -1.4,
'range': [-19.6376, 200.725],
'colorwedge':T, 'colormap': 'Rainbow 2', }, out='VYCMa_658_0_1966_ph2.clean.png',
zoom=2)
imview(raster={'file': 'VYCMa_658_0_1966_ap2.clean.image','scaling': -1.4,
'range': [-19.6376, 200.725],
'colorwedge':T, 'colormap': 'Rainbow 2', }, out='VYCMa_658_0_1966_ap2.clean.png',
zoom=2)
imview(raster={'file': 'VYCMa_658_0_1966_ap3.clean.image','scaling': -1.4,
'range': [-19.6376, 200.725],
'colorwedge':T, 'colormap': 'Rainbow 2', }, out='VYCMa_658_0_1966_ap3.clean.png',
zoom=2)