#-- Self-calibration on water maser peak, applied to all channels/spw. #-- Imaging of continuum, continuum subtraction, make cubes for #-- water maser and another line in the other spw. #-- Test channel 1:636 containing extended emission is also imaged at some #-- stages to check the calibration, compare these images by hand. #-- 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=False #-- * When making an image for self-cal, make sure that whatever clean #-- * version you use puts the required model properly into the MS, #-- * and this is not overwritten by making test images. If in doubt, #-- * use 'ft'. #-- 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=True refantenna='DV15' thesteps=[] step_title = {0: 'List the data set', 1: 'Plot the visibility spectrum and peak', 2: 'Make first image of maser peak and test channel', 3: 'Self-calibrate cycle 1, phase-only, re-image and test channel', 4: 'Self-calibrate cycle 2, phase-only, re-image', 5: 'Self-calibrate cycle 3, amp & phase, re-image', 6: 'Self-calibrate cycle 4, amp & phase, re-image and test channel', 7: 'Make low-res cubes to identify continuum \n(this step can be omitted if you trust our channel selection)', 8: 'Image continuum', 9: 'Subtract continuum, image test channel', 10: 'Image a line in spw 0', 11: 'Image a line in spw 1', } # 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 """ # PREPARATION # Select just one of the 3 EBs, with Tsys, WVR, bandpass and phase-ref # solutions applied. Shift to constant LSR velocity and average target basename=['uid___A002_X6d0f96_X1de2'] refantenna='DV15' vislist=[] for name in basename: for s in range(2): os.system('rm -rf '+name+'spw'+str(s)+'_VYCMa_325.ms') mstransform(vis=name+'.ms.split', spw=str(s), outputvis=name+'spw'+str(s)+'_VYCMa_325.ms',datacolumn='corrected', field='V*', regridms=True, width=2, outframe='lsrk', timeaverage=True, timebin='12.1s') vislist.append(name+'spw'+str(s)+'_VYCMa_325.ms') concat(vis=vislist, concatvis='X1de2_VYCMa_325.ms') """ mystep = 0 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm X1de2_VYCMa_325_listobs.txt') listobs(vis='X1de2_VYCMa_325.ms', listfile='X1de2_VYCMa_325_listobs.txt', verbose=True) print '<< Printed listobs output to X1de2_VYCMa_325_listobs.txt' mystep = 1 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] for s in ['0','1']: os.system('rm -rf X1de2_VYCMa_325_spw'+s+'_1st-vis-spectrum.png') plotms(vis='X1de2_VYCMa_325.ms', xaxis='frequency', yaxis='amp', selectdata=True, spw=s, avgtime='1e8',avgscan=True, coloraxis='baseline', showgui=False, plotfile='X1de2_VYCMa_325_spw'+s+'_1st-vis-spectrum.png') # You will notice the very bright maser line in spw 0:981 os.system('rm -rf X1de2_VYCMa_325_0-981_1st_amp.png') plotms(vis='X1de2_VYCMa_325.ms', xaxis='time', yaxis='amp', selectdata=True, antenna='DV15&*', spw='0:981', coloraxis='baseline', showgui=False, plotfile='X1de2_VYCMa_325_0-981_1st_amp.png') os.system('rm -rf X1de2_VYCMa_325_0-981_1st_phase.png') plotms(vis='X1de2_VYCMa_325.ms', xaxis='time', yaxis='phase', selectdata=True, antenna='DV15&*', spw='0:981', coloraxis='baseline', showgui=False, plotfile='X1de2_VYCMa_325_0-981_1st_phase.png') # * DA50 has fast phase but it turns out it can be corrected ###################################################### # Make first image of maser peak, omitting dubious data 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 X1de2_VYCMa_325_0_981.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_0_981.clean', field='Vy Cma', spw='0:981', antenna='!DA50', mode='channel',start=981,nchan=1, mask='ellipse[[258pix,281pix],[63pix,59pix],90deg]', cell='0.01arcsec',npercycle=10,imsize=512, interactive=cleaninteractive, niter=100) # Get image statistics imin='X1de2_VYCMa_325_0_981.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) # rms 4.072, peak 184.486, snr 45 ### image resolved channel in other spw as test: chan 1:636 os.system('rm -rf X1de2_VYCMa_325_1_636.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_1_636.clean', field='Vy Cma', spw='1:636', antenna='!DA50', mode='channel',start=636,nchan=1, cell='0.01arcsec',npercycle=10,imsize=512, interactive=cleaninteractive, niter=1000) ###################################################### # 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 X1de2_VYCMa_325_981.ph1') gaincal(vis = 'X1de2_VYCMa_325.ms', field='Vy Cma', refant=refantenna, caltable='X1de2_VYCMa_325_981.ph1', spw='0:981', calmode='p', solint='60s') if (calplots): os.system('rm -rf X1de2_VYCMa_325_981.ph1*.png') for ant in range(3): plotcal(caltable='X1de2_VYCMa_325_981.ph1', xaxis = 'time', yaxis = 'phase', antenna=str(ant*8)+'~'+str(ant*8+7), subplot=421, fontsize=7.0, iteration='antenna', figfile='X1de2_VYCMa_325_981.ph1_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') applycal(vis = 'X1de2_VYCMa_325.ms', field='Vy Cma', spw='0,1', spwmap=[0,0], gaintable='X1de2_VYCMa_325_981.ph1', calwt = False, applymode='calonly', flagbackup = False) os.system('rm -rf X1de2_VYCMa_325_0_981_ph1.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_0_981_ph1.clean', field='Vy Cma', spw='0:981', mode='channel',start=981,nchan=1, cell='0.01arcsec',npercycle=50,imsize=512, psfmode='hogbom', mask='ellipse[[258pix,281pix],[63pix,59pix],90deg]', interactive=cleaninteractive, niter=250) # Get image statistics imin='X1de2_VYCMa_325_0_981_ph1.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) #rms 1.319, peak 220.227, snr 166 os.system('rm -rf X1de2_VYCMa_325_1_636_ph1.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_1_636_ph1.clean', field='Vy Cma', spw='1:636', antenna='!DA50', mode='channel',start=636,nchan=1, cell='0.01arcsec',npercycle=10,imsize=512, interactive=cleaninteractive, niter=1000) ###################################################### # Self-calibrate cycle 2, phase-only, shorter solint, using all data 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 X1de2_VYCMa_325_981.ph2') gaincal(vis = 'X1de2_VYCMa_325.ms', field='Vy Cma', refant=refantenna, caltable='X1de2_VYCMa_325_981.ph2', spw='0:981', calmode='p', solint='36s') if (calplots): os.system('rm -rf X1de2_VYCMa_325_981.ph2*.png') for ant in range(3): plotcal(caltable='X1de2_VYCMa_325_981.ph2', xaxis = 'time', yaxis = 'phase', antenna=str(ant*8)+'~'+str(ant*8+7), subplot=421, fontsize=7.0, iteration='antenna', figfile='X1de2_VYCMa_325_981.ph2_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') applycal(vis = 'X1de2_VYCMa_325.ms', field='Vy Cma', spw='0,1', spwmap=[0,0], gaintable='X1de2_VYCMa_325_981.ph2', calwt = False, applymode='calonly', flagbackup = False) os.system('rm -rf X1de2_VYCMa_325_0_981_ph2.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_0_981_ph2.clean', field='Vy Cma', spw='0:981', mode='channel',start=981,nchan=1, cell='0.01arcsec',npercycle=50,imsize=512, psfmode='hogbom', mask='ellipse[[256pix,277pix],[68pix,57pix],90deg]', interactive=cleaninteractive, niter=700) # Get image statistics imin='X1de2_VYCMa_325_0_981_ph2.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) # rms 0.963, peak 220.139, snr 228 ###################################################### # Self-calibrate cycle 3, a&p, using 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 X1de2_VYCMa_325_981.ap1') gaincal(vis = 'X1de2_VYCMa_325.ms', field='Vy Cma', refant=refantenna, caltable='X1de2_VYCMa_325_981.ap1', spw='0:981', gaintable='X1de2_VYCMa_325_981.ph2', spwmap=[0,0], calmode='ap', solint='120s') if (calplots): os.system('rm -rf X1de2_VYCMa_325_981.ap1*.png') for ant in range(3): plotcal(caltable='X1de2_VYCMa_325_981.ap1', xaxis = 'time', yaxis = 'amp', antenna=str(ant*8)+'~'+str(ant*8+7), subplot=421, fontsize=7.0, iteration='antenna', figfile='X1de2_VYCMa_325_981.ap1_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') applycal(vis = 'X1de2_VYCMa_325.ms', field='Vy Cma', spw='0,1', spwmap=[[0,0],[0,0]], applymode='calonly', gaintable=['X1de2_VYCMa_325_981.ph2','X1de2_VYCMa_325_981.ap1'], calwt = False, flagbackup = False) # os.system('rm -rf X1de2_VYCMa_325_0_981_ap1.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_0_981_ap1.clean', field='Vy Cma', spw='0:981', mode='channel',start=981,nchan=1, cell='0.01arcsec',npercycle=50,imsize=512, psfmode='hogbom', mask='ellipse[[272pix,288pix],[129pix,96pix],0deg]', interactive=cleaninteractive, niter=1000) # Get image statistics imin='X1de2_VYCMa_325_0_981_ap1.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) # rms 0.418, peak 224.552, snr 537 ###################################################### # Self-calibrate cycle 4, ap, short solint, using all data 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 X1de2_VYCMa_325_981.ap2') gaincal(vis = 'X1de2_VYCMa_325.ms', field='Vy Cma', refant=refantenna, caltable='X1de2_VYCMa_325_981.ap2', spw='0:981', gaintable='X1de2_VYCMa_325_981.ph2', spwmap=[0,0], calmode='ap', solint='36s') if (calplots): os.system('rm -rf X1de2_VYCMa_325_981.ap2*.png') for ant in range(3): plotcal(caltable='X1de2_VYCMa_325_981.ap2', xaxis = 'time', yaxis = 'amp', antenna=str(ant*8)+'~'+str(ant*8+7), subplot=421, fontsize=7.0, iteration='antenna', figfile='X1de2_VYCMa_325_981.ap2_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') applycal(vis = 'X1de2_VYCMa_325.ms', field='Vy Cma', spw='0,1', spwmap=[[0,0],[0,0]], gaintable=['X1de2_VYCMa_325_981.ph2','X1de2_VYCMa_325_981.ap2'], calwt = False, # try true? NB flagging failed flagbackup = False) # os.system('rm -rf X1de2_VYCMa_325_0_981_ap2.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_0_981_ap2.clean', field='Vy Cma', spw='0:981', mode='channel',start=981,nchan=1, cell='0.01arcsec',npercycle=100,imsize=512, psfmode='hogbom', mask='ellipse[[272pix,288pix],[129pix,96pix],0deg]', interactive=cleaninteractive, niter=2000) # Get image statistics imin='X1de2_VYCMa_325_0_981_ap2.clean.image' noise=imstat(imagename=imin, box='250,50,500,200')['rms'][0] peak=imstat(imagename=imin, box='200,200,330,330')['max'][0] print 'rms %8.3f, peak %8.3f, snr %5i' % (noise, peak, peak/noise) # rms 0.273, peak 223.915, snr 821 os.system('rm -rf X1de2_VYCMa_325_1_636_ap2.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_1_636_ap2.clean', field='Vy Cma', spw='1:636', antenna='!DA50', mode='channel',start=636,nchan=1, cell='0.01arcsec',npercycle=10,imsize=512, interactive=cleaninteractive, niter=1000) ##################################### # Image whole cube at lower resolution to identify line-free channels # This step can be omitted if you trust our selection below. mystep = 7 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] for s in ['0','1']: os.system('rm -rf X1de2_VYCMa_325_spw'+s+'_5ap2.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_spw'+s+'_5ap2.clean', field='Vy Cma', spw=s, mode='channel',width=5, cell='0.02arcsec',imsize=512,usescratch=True, psfmode='hogbom', interactive=cleaninteractive, niter=1000) # The following is our best guess at a selection contchans='0:5~224;365~379;425~479;615~739;1185~1249;1460~1569;1615~1699;1820~1919,1:0~214;440~494;580~609;750~789;830~1049;1305~1364;1490~1529;1610~1659;1775~1829;1815~1834;1895~1919' ############################################## # Image continuum 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 X1de2_VYCMa_325_contap2.clean*') clean(vis = 'X1de2_VYCMa_325.ms', imagename='X1de2_VYCMa_325_contap2.clean', field='Vy Cma', spw=contchans, mode='mfs', multiscale=[0,4,8,12] , mask=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'], cell='0.01arcsec',imsize=512, psfmode='hogbom',interactive=cleaninteractive) exportfits(imagename='X1de2_VYCMa_325_contap2.clean.image', fitsimage='X1de2_VYCMa_325_contap2.clean.fits') ##################################### # Subtract 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 X1de2_VYCMa_325.ms.contsub') uvcontsub(vis = 'X1de2_VYCMa_325.ms',field='V*',fitorder=1, combine='spw',spw='',fitspw=contchans) os.system('rm -rf X1de2_VYCMa_325_1_636_ap2_sub.clean*') clean(vis = 'X1de2_VYCMa_325.ms.contsub', imagename='X1de2_VYCMa_325_1_636_ap2_sub.clean', field='Vy Cma', spw='1:636', antenna='!DA50', mode='channel',start=636,nchan=1, cell='0.01arcsec',npercycle=10,imsize=512, interactive=cleaninteractive, niter=1000) #---------------------------------------------------------------------------------- #----- Image maser and line emission and make moment maps ------------------------- #---------------------------------------------------------------------------------- # Some lines from Kaminski et al. 2013 ApJS 209 38 are imaged below. # However, note that it is possible that the named line may be blended with other lines, # especially for the thermal lines. Conversely, it is possible (especially # for the maser) that faint wings exceed the range imaged. # ##################################### # Image H2O 5(1,5)-4(2,2) line spw 0 mystep = 10 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] # Line frequencies taken from Kaminski et al. 2013 ApJS 209 38 # High resolution for maser os.system('rm -rf X1de2_VYCMa_325.15_H2O.clean*') clean(vis = 'X1de2_VYCMa_325.ms.contsub', imagename='X1de2_VYCMa_325.15_H2O.clean', field='Vy Cma', spw='0', mode='velocity',start='-38km/s',nchan=134, restfreq='325.152919GHz', cell='0.01arcsec',imsize=512, psfmode='hogbom', weighting='briggs',robust=0.5, threshold='10mJy', mask=['ellipse[[218pix,255pix],[64pix,52pix],90deg]','ellipse[[281pix,359pix],[140pix,70pix],20deg]'], interactive=cleaninteractive, usescratch=True, niter=1000) # Make moment map thisimage='X1de2_VYCMa_325.15_H2O.clean.image' chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='3') rms1= chanstat['rms'][0] chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='130') rms2= chanstat['rms'][0] rms=0.5*(rms1+rms2) print 'rms in a channel = '+str(rms) os.system('rm -rf '+thisimage+'.mom0') immoments(imagename = thisimage, moments = [0], axis = 'spectral', includepix = [rms*2,1000.], outfile = thisimage+'.mom0') # Export as fits exportfits(imagename=thisimage+'.mom0', fitsimage=thisimage+'.mom0.fits') exportfits(imagename='X1de2_VYCMa_325.15_H2O.clean.image', fitsimage='X1de2_VYCMa_325.15_H2O.clean.fits') # Note that since the dynamic range of the maser cube is 500~1000 around the peak, # implying a noise level more than an order of magnitude # higher than that for the fainter line wings, we don't make a first moment map. # However, it may be possible to obtain a meaningful first moment map with # additional data selection and/or blanking ##################################### # Image NaCl 24-23 line spw 1 mystep = 11 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] # Line frequencies taken from Kaminski et al. 2013 ApJS 209 38 # Lower resolution for weaker lines os.system('rm -rf X1de2_VYCMa_312.11_NaCl.clean*') clean(vis = 'X1de2_VYCMa_325.ms.contsub', imagename='X1de2_VYCMa_312.11_NaCl.clean', field='Vy Cma', spw='1', mode='velocity', start='-8km/s',nchan=30,width='2km/s', restfreq='312.1099004GHz', cell='0.01arcsec',imsize=512, psfmode='hogbom', threshold='4mJy', mask='ellipse[[284pix,263pix],[150pix,115pix],20deg]', interactive=cleaninteractive, usescratch=True, niter=1000) # Make moment maps thisimage='X1de2_VYCMa_312.11_NaCl.clean.image' chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='5') rms1= chanstat['rms'][0] chanstat=imstat(imagename=thisimage,box='250,50,500,100',chans='24') rms2= chanstat['rms'][0] rms=0.5*(rms1+rms2) print 'rms in a channel = '+str(rms) os.system('rm -rf '+thisimage+'.mom0') immoments(imagename = thisimage, moments = [0], axis = 'spectral', includepix = [rms*3,100.], outfile = thisimage+'.mom0') os.system('rm -rf '+thisimage+'.mom1') immoments(imagename = thisimage, moments = [1], axis = 'spectral', includepix = [rms*4,100.], outfile = thisimage+'.mom1') # Export as fits exportfits(imagename=thisimage+'.mom0', fitsimage=thisimage+'.mom0.fits') exportfits(imagename=thisimage+'.mom1', fitsimage=thisimage+'.mom1.fits') exportfits(imagename='X1de2_VYCMa_312.11_NaCl.clean.image', fitsimage='X1de2_VYCMa_312.11_NaCl.clean.image.fits') #---------------------------------------------------------------------------------- #----- End of imaging script. #----------------------------------------------------------------------------------