# I-TRAIN with the European ARC Network # Improving image fidelity through self-calibration # Richards, Moravec, Perez-Sanchez, Toribio # 2021/05/24 #============================================================================================== #-- This script was prepared for a tutorial on improving image fidelity through self-calibration #-- of ALMA data. #-- The dataset used in this script corresponds to Science Verification data of VY CMA in Band 7 #-- that was preprocessed for the purpose #-- Link to the original dataset VY CMA from SV data: #-- https://almascience.nrao.edu/almadata/sciver/VYCMaBand7/ #-- Preparation steps applied the original dataset can be found commented below. #-- The script is organized in STEPS. Search for the variable "step_title" that holds all the steps #-- executed by the script. To execute the full script, use the CASA command 'execfile'. If you #-- would only like to run certain steps, define the variable 'mysteps' and then execute the script #-- with 'execfile', eg, to execute steps 0 & 1 only, run in the CASA terminal: #-- mysteps=[0,1] #-- execfile('itrain-selfcal.py') #-- In this script self-calibration of the maser is performed in several cycles. #-- The script dismisses antenna DA50 from imaging and gaincal in the first self-calibration cycles #-- to recover it in the end. #-- If AnalysisUils is installed, the duration of scans is computed in step 0. #-- AnalysisUtils can be installed from https://casaguides.nrao.edu/index.php?title=Analysis_Utilities """ # PREPARATION ALREADY DONE ON THE ORIGINAL SV DATASET # For future reference, we include here the steps that have been applied to the orignal # Science Verification dataset of VY CMA to obtain the measurement set used in this tutorial. # Starting with the Science Verification data of VY CMA # 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') #Flag antenna DA50 flagdata(vis=vis, mode='manual', antenna='DA50', flagbackup=True) """ #============================================================================ # Python packages #============================================================================ import os import sys import matplotlib import matplotlib.pyplot as plt from scipy import stats import numpy as np #============================================================================ # PARAMETER DEFINITION #============================================================================ #-- Measurement set visname = 'X1de2_VYCMa_325' vis = visname+'.ms' field = 'Vy Cma' #-- reference antenna refantenna ='DV15' #-- Continuum channels selection #-- This selection was made after calibration on the maser and imaging but #-- you could make a selection from the visibility spectrum before self-cal contchans='0:5~224;365~379;425~479;615~739;1185~1249;1460~1569;1615~1699;1820~1919,1:0~214;440~475;487~494;580~609;750~789;830~1049;1305~1364;1490~1529;1610~1659;1775~1829;1815~1834;1895~1919' #-- Maser peak channel ch_peak='0:981' ch_peakname='0-981' startch_peak=981 #-- Test channel that will be imaged for comparison ch_test='1:636' ch_testname='1-636' startch_test=636 #-- 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=True #-- Cleaning mask #cleaningmask='ellipse[[256pix,277pix],[68pix,57pix],90deg]' cleaningmask='ellipse[[272pix,288pix],[129pix,96pix],0deg]', cleaningmask_cont=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'], cleaningmask_maser=['ellipse[[218pix,255pix],[64pix,52pix],90deg]','ellipse[[281pix,359pix],[140pix,70pix],20deg]'] #=========================================================================== # FUNCTIONS #=========================================================================== # Useful functions for our purposes are defined here def get_im_stats(im_name): # Calculate image statistics noise=imstat(imagename=im_name,box='250,50,500,200')['rms'][0] #region representative of the image RMS, large enough and without source signal peak=imstat(imagename=im_name,box='200,200,330,330')['max'][0] #region including our target print('rms {0:.3f}, peak {1:.3f}, snr {2:.0f}'.format(noise, peak, peak/noise)) def plot_gaincal_table(caltable): # make plots for antenna triplets that will be saved in png files for ant in range(3): plotms(vis=caltable,xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=4, gridcols=2,coloraxis='spw', xaxisfont = 7, yaxisfont = 7, highres = True, antenna=str(ant*8)+'~'+str(ant*8+7), showgui = False, plotfile=caltable+'_ant'+str(ant*8)+'-'+str(ant*8+7)+'.png') def plot_gaincal_snr_dist(path, visname, selfcal_cycle, solints): # Make plot of SNR of gaintables for a list of solution intervals # path: path to caltables; plot will be saved in that path too # solints: list of solution intervals (list of strings) print("Percentile of score for snr=6, i.e. percentage of solutions with snr <=6:") for ss in solints: tb.open(path+'/'+visname+'.'+selfcal_cycle+'.solint_'+ss+'.tb') snr = tb.getcol( 'SNR' ).ravel() tb.close() if(matplotlib.__version__>="2.0"): plt.hist( snr, bins=50, density=True, histtype='step', label=ss ) else: plt.hist( snr, bins=50, normed=True, histtype='step', label=ss ) print( 'P(<=6) = {0} ({1})'.format( stats.percentileofscore( snr, 6 ), ss ) ) plt.legend( loc='upper right' ) plt.xlabel( 'SNR' ) plt.savefig(path+'/'+selfcal_cycle+'_SNR_hist_solint_all.png') def explore_scan_duration(vis): # Print information related to the time on source using AnalysisUtils task: https://casaguides.nrao.edu/index.php?title=TimeOnSource # This function requires Analysisi Utilities to be installed: https://casaguides.nrao.edu/index.php?title=Analysis_Utilities try: print("Computing duration of scans") print("---------------------------------------------") import analysisUtils as au mydict = au.timeOnSource(vis, verbose=False) print("---------------------------------------------") print("Number of scans:", mydict[0]['num_of_scans']) print("Minutes of science per scan:", mydict['minutes_on_science_per_scan']) except: print("This calculation requires the package AnalysisUtils to be installed (https://casaguides.nrao.edu/index.php?title=Analysis_Utilities) ") print("Alternatively, you can estimate scan duration times from the listobs file") #=========================================================================== # STEPS #=========================================================================== # List of steps executed by this script thesteps=[] step_title = {0: 'List the data set, unflag antenna DA50 and plot antennas and visibility spectrum, calculate duration of scans', 1: 'Make dirty image of continuum', ### INITIAL MODEL 2: 'Make an initial, conservative cleaning', 3: 'Check and save model', ### FIRST ROUND OF SELF-CALIBRATION - PHASE 4: 'Calculate gain solution table - phase-only, long solution interval', 5: 'Explore different solution intervals', 6: '[ADVANCED] Calculate SNR of the different solution intervals', 7: 'Apply calibration table', 8: 'Make second, conservative cleaning and save model', ### SECOND ROUND OF SELF-CALIBRATION - PHASE 9: 'Explore different solution intervals', 10: '[ADVANCED] Calculate SNR of the different solution intervals', 11: 'Calculate gain solution table - phase-only, shorter solution interval applying round 1 table on-the-fly', 12: 'Apply calibration tables', 13: 'Make image of continuum and save model', ### THIRD ROUND OF SELF-CALIBRATION - AMPLITUDE & PHASE 14: 'Calculate gain solution table - amplitude and phase, long solution interval', 15: 'Apply calibration tables', 16: 'Make image of continuum and save model', ### FOURTH ROUND OF SELF-CALIBRATION - AMPLITUDE & PHASE 17: 'Calculate gain solution table - amplitude and phase, short solution interval', 18: 'Apply calibration table', ### FINAL IMAGE 19: 'Make channel image and save model', 20: '[OPTIONAL] Make low-resolution cube to search for line-free channels ', 21: 'Make image of continuum', ### CONTINUUM SUBTRACTION 22: 'Subtract continuum', ## ANALYSIS 23: 'Maser line cube and moment map', 24: 'Image NaCl 24-23 and moment maps' } # The Python variable 'mysteps' will control which steps are executed when you start the script using # execfile('[This script].py') # e.g. to execute only steps 2, 3, and 4 of the script define first 'mysteps' and then execute the script: # mysteps = [2,3,4] # execfile('script.py') # 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 = list(range(0,len(step_title))) print('Executing all steps: ', thesteps) #--------- mystep = 0 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # print listobs os.system('rm '+visname+'_listobs.txt') listobs(vis=vis, listfile=visname+'_listobs.txt', verbose=True) # Unflag antenna DA50 flagdata(vis=vis,mode='unflag', antenna='DA50', flagbackup=True) # plot antenna locations os.system('rm '+visname+'*plotants*png') plotants(vis, figfile=visname+'_plotants.png') plotants(vis, logpos=True, figfile=visname+'_plotants_log.png') # Plots of amplitude vs. frequency # --You will notice the very bright maser line in spw 0:981 # --even for continuum it is worth plotting this to reveal bad data for s in ['0','1']: os.system('rm -rf '+visname+'_spw'+s+'_vis-spectrum.png') plotms(vis = vis, xaxis='frequency', yaxis='amp', selectdata=True, spw=s, avgtime='1e8',avgscan=True,avgbaseline=True, coloraxis='baseline', showgui=False, plotfile=visname+'_spw'+s+'_vis-spectrum.png') # Plot visibilities of channel 0:981 os.system('rm -rf '+visname+'_'+ch_peakname+'_1st_amp.png') plotms(vis=vis, xaxis='time', yaxis='amp', selectdata=True, antenna=refantenna+'&*', spw=ch_peak, coloraxis='baseline', showgui=True, plotfile=visname+'_'+ch_peakname+'_1st_amp.png') os.system('rm -rf '+visname+'_'+ch_peakname+'_1st_phase.png') plotms(vis=vis, xaxis='time', yaxis='phase', selectdata=True, antenna=refantenna+'&*', spw=ch_peak, coloraxis='baseline', showgui=True, plotfile=visname+'_'+ch_peakname+'_1st_phase.png') # Compute duration of scans (requires AnalysisUtils) explore_scan_duration(vis) # * DA50 has fast phase but it turns out it can be corrected # * The antenna DA50 will be excluded in tclean in the first cycles of self-calibration # * and it will be recovered later. #--------- mystep = 1 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) ## Make a first dirty image of the maser peak to get a sense of the structure of the object imagename = visname + '_'+ch_peakname+'.dirty' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, antenna = '!DA50', #DA50 is excluded spw=ch_peak, specmode='cube',start=startch_peak,nchan=1, cell='0.01arcsec', imsize=512, niter=0, interactive=False) # view image imview(imagename+'.image') ### INITIAL MODEL #--------- mystep = 2 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) #Delete existing models and calibration tables to make sure we have a neat start: delmod(vis=vis, scr=True) clearcal(vis=vis) # make an initial, conservative clean imagename = visname +'_'+ch_peakname +'.init.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, antenna = '!DA50', #DA50 is excluded spw=ch_peak, specmode='cube',start=startch_peak,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') # Image a test channel for sanity checks imagename = visname +'_'+ch_testname +'.init.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, antenna = '!DA50', #DA50 is excluded spw=ch_test, specmode='cube',start=startch_test,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') #--------- mystep = 3 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) #model name modelname=visname+'_'+ch_peakname+'.init.clean.model' # check that model has saved plotms(vis=vis, xaxis='UVwave', yaxis='amp', ydatacolumn='model',showgui=False,plotfile=modelname+'.png') # force model to save ft(vis=vis,model=modelname,usescratch=True) # check that model has saved after ft plotms(vis=vis, xaxis='UVwave', yaxis='amp', ydatacolumn='model',showgui=False,plotfile=modelname+'_ft.png') ### FIRST ROUND OF SELF-CALIBRATION - PHASE #--------- mystep = 4 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # calculate gain table for solint='inf' = scan length solint='60s' caltable=visname+'_'+ch_peakname+'.ph1.solint_'+solint+'.tb' os.system('rm -rf '+caltable) gaincal(vis = vis, field= field, refant=refantenna, caltable=caltable, spw=ch_peak, calmode='p', solint=solint, gaintype='G', minsnr=3) # view gain table for all antennas plotms(vis=caltable,xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3,coloraxis='spw') #--------- mystep = 5 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) ### COMMANDS FOR EXPLORING THE SOLUTION INTERVAL # # We recommend investigating the phases with various solution intervals to get a sense of the interval on which the solutions vary # The following loop calculates gaincal solutions for a list of intervals and makes corresponding plots # The output is saved in a separate folder selfcal_cycle = 'ph1_checks' solint_all = ['int', '20s', '40s', '60s', '80s', '160s', 'inf'] for solint in solint_all: print('Solint:', solint) caltable = visname+'_'+ch_peakname+'.'+selfcal_cycle+'.solint_'+solint+'.tb' gaincal(vis=vis,caltable=caltable,solint=solint,refant=refantenna,spw=ch_peak,calmode='p',gaintype='G',minsnr=3) # make plots for antenna triplets that will be saved in png files plot_gaincal_table(caltable) os.system('rm -r '+selfcal_cycle) if not os.path.exists(selfcal_cycle): os.makedirs(selfcal_cycle) os.system('mv *.'+selfcal_cycle+'*tb* '+selfcal_cycle+'/') print("Check output of this step in folder: "+selfcal_cycle) #--------- mystep = 6 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # [ADVANCED] # Calculate the distribution of SNR of the gaincal tables for different solution intervals # using the tables generated in the previous step; plots are moved to the same output folder as in the step above selfcal_cycle = 'ph1_checks' solint_all = ['int', '20s', '40s', '60s', '80s', '160s', 'inf'] try: path = selfcal_cycle + '/' plot_gaincal_snr_dist(path, visname+'_'+ch_peakname, selfcal_cycle, solint_all) print("Check output of this step in folder: "+selfcal_cycle) except: print("You need to run the previous step for the same list of solution intervals") #--------- mystep = 7 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # apply the solutions to the MS solint='60s' caltable=visname+'_'+ch_peakname+'.ph1.solint_'+solint+'.tb' applycal(vis = vis, field= field, spw='0,1', spwmap=[0,0], gaintable=caltable, calwt = False, applymode='calonly', flagbackup = False) #--------- mystep = 8 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # make a second, conservative clean imagename = visname + '_'+ ch_peakname + '.ph1.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, antenna = '!DA50', #DA50 is excluded spw=ch_peak, specmode='cube',start=startch_peak,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') #force model to save modelname=imagename+'.model' ft(vis=vis,model=modelname,usescratch=True) # Image a test channel for sanity checks imagename = visname +'_'+ch_testname +'.ph1.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, antenna = '!DA50', #DA50 is excluded spw=ch_test, specmode='cube',start=startch_test,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') ### SECOND ROUND OF SELF-CALIBRATION - PHASE #--------- mystep = 9 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) ### COMMANDS FOR EXPLORING THE SOLUTION INTERVAL # # We recommend investigating the phases with various solution intervals to get a sense of the interval on which the solutions vary # The following loop calculates gaincal solutions for a list of intervals and makes corresponding plots # The output is saved in a separate folder selfcal_cycle = 'ph2_checks' solint_all = ['int', '20s', '40s', '60s', '80s', '160s', 'inf'] for solint in solint_all: print('Solint:', solint) solint_1='60s' caltable = visname+'_'+ch_peakname+'.'+selfcal_cycle+'.solint_'+solint+'.tb' gaincal(vis=vis,caltable=caltable,solint=solint,refant=refantenna,spw=ch_peak, gaintable = [visname + '_'+ch_peakname + '.ph1.solint_'+solint_1+'.tb'], spwmap=[0,0], calmode='p',gaintype='G',minsnr=3) # make plots for antenna triplets that will be saved in png files plot_gaincal_table(caltable) os.system('rm -r '+selfcal_cycle) if not os.path.exists(selfcal_cycle): os.makedirs(selfcal_cycle) os.system('mv *.'+selfcal_cycle+'*tb* '+selfcal_cycle+'/') print("Check output of this step in folder: "+selfcal_cycle) #--------- mystep = 10 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # [ADVANCED] # Calculate the distribution of SNR of the gaincal tables for different solution intervals # using the tables generated in the previous step; plots are moved to the same output folder as in the step above selfcal_cycle = 'ph2_checks' solint_all = ['int', '20s', '40s', '60s', '80s', '160s', 'inf'] try: path = selfcal_cycle + '/' plot_gaincal_snr_dist(path, visname+'_'+ch_peakname, selfcal_cycle, solint_all) print("Check output of this step in folder: "+selfcal_cycle) except: print("You need to run the previous step for the same list of solution intervals") #--------- mystep = 11 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # calculate gain table for a new solint while applying the table from round 1 'on the fly' solint_1='60s' solint='20s' caltable = visname + '_' +ch_peakname+ '.ph2.solint_'+solint+'.tb' os.system('rm -rf '+caltable) gaincal(vis = vis, field= field, refant=refantenna, caltable=caltable, spw=ch_peak, gaintable = [visname + '_'+ch_peakname+'.ph1.solint_'+solint_1+'.tb'], spwmap=[0,0], calmode='p', solint=solint, gaintype='G', minsnr=3) # view gain table for all antennas plotms(vis=caltable,xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3,coloraxis='spw') #--------- mystep = 12 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # apply the cumulative solutions to the MS solint_1='60s' solint_2='20s' applycal(vis = vis, field= field, spw='0,1', gaintable=[visname+'_'+ch_peakname+ '.ph1.solint_'+solint_1+'.tb', visname+'_'+ ch_peakname+'.ph2.solint_'+solint_2+'.tb'], spwmap=[[0,0],[0,0]], calwt = False, applymode='calonly', flagbackup = False) #--------- mystep = 13 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # If you would like to compare what the second round of phase self-cal accomplished compared to the first # make a third, conservative clean # Antenna DA50 may also be included from now on imagename = visname + '_'+ ch_peakname + '.ph2.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, #antenna = '!DA50', #DA50 is excluded spw=ch_peak, specmode='cube',start=startch_peak,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') #force model to save modelname=imagename+'.model' ft(vis=vis,model=modelname,usescratch=True) # Image a test channel for sanity checks imagename = visname +'_'+ch_testname +'.ph2.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, #antenna = '!DA50', #DA50 is excluded spw=ch_test, specmode='cube',start=startch_test,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') ### THIRD ROUND OF SELF-CALIBRATION - AMPLITUDE & PHASE #--------- mystep = 14 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # long solution interval, applying previous solutions on the fly solint='120s' caltable = visname + '_'+ch_peakname + '.ap1.solint_'+solint+'.tb' # apply the cumulative solutions to the MS solint_1='60s' solint_2='20s' os.system('rm -rf '+caltable) gaincal(vis = vis, field= field, refant=refantenna, caltable=caltable, gaintable = [visname + '_'+ch_peakname + '.ph1.solint_'+solint_1+'.tb', visname + '_'+ch_peakname +'.ph2.solint_'+solint_2+'.tb' ], spwmap=[[0,0],[0,0]], spw=ch_peak, calmode='ap', solint=solint, gaintype='G', minsnr=3) # view gain table for all antennas plotms(vis=caltable,xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3,coloraxis='spw') #--------- mystep = 15 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # apply the cumulative solutions to the MS caltable=visname+'_'+ch_peakname+'.ap1.solint_'+solint+'.tb' solint_1='60s' solint_2='20s' solint_3='120s' applycal(vis = vis, field= field, spw='0,1', gaintable = [visname + '_'+ch_peakname + '.ph1.solint_'+solint_1+'.tb', visname + '_'+ch_peakname +'.ph2.solint_'+solint_2+'.tb', visname + '_'+ch_peakname +'.ap1.solint_'+solint_3+'.tb' ], spwmap=[[0,0],[0,0],[0,0]], calwt = False, applymode='calonly', flagbackup = False) #--------- mystep = 16 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # If you would like to compare what the second round of phase self-cal accomplished compared to the first # make yet another, conservative clean imagename = visname + '_'+ ch_peakname + '.ap1.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, #antenna = '!DA50', #DA50 is excluded spw=ch_peak, specmode='cube',start=startch_peak,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') #force model to save modelname=imagename+'.model' ft(vis=vis,model=modelname,usescratch=True) # Image a test channel for sanity checks imagename = visname +'_'+ch_testname +'.ap1.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, #antenna = '!DA50', #DA50 is excluded spw=ch_test, specmode='cube',start=startch_test,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') ### FOURTH ROUND OF SELF-CALIBRATION - AMPLITUDE & PHASE #--------- mystep = 17 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # calculate gain table for a new solint while applying previous solutions 'on the fly' solint='40s' caltable = visname + '_'+ch_peakname+'.ap2.solint_'+solint+'.tb' solint_1='60s' solint_2='20s' solint_3='120s' os.system('rm -rf '+caltable) gaincal(vis = vis, field= field, refant=refantenna, caltable=caltable, gaintable = [visname + '_'+ch_peakname + '.ph1.solint_'+solint_1+'.tb', visname + '_'+ch_peakname +'.ph2.solint_'+solint_2+'.tb', visname + '_'+ch_peakname +'.ap1.solint_'+solint_3+'.tb' ], spwmap=[[0,0],[0,0],[0,0]], spw=ch_peak, calmode='ap', solint=solint, #gaintype='T', # notice gaintype option minsnr=3) # view gain table for all antennas plotms(vis=caltable,xaxis='time',yaxis='phase',iteraxis='antenna',gridrows=3, gridcols=3,coloraxis='spw') #--------- mystep = 18 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # apply the cumulative solutions to the MS solint_1='60s' solint_2='20s' solint_3='120s' solint_4='40s' applycal(vis = vis, field= field, spw='0,1', gaintable = [visname + '_'+ch_peakname + '.ph1.solint_'+solint_1+'.tb', visname + '_'+ch_peakname +'.ph2.solint_'+solint_2+'.tb', visname + '_'+ch_peakname +'.ap1.solint_'+solint_3+'.tb', visname + '_'+ch_peakname +'.ap2.solint_'+solint_4+'.tb' ], spwmap=[[0,0],[0,0],[0,0],[0,0]], calwt = False, flagbackup = False) ### FINAL IMAGE #--------- mystep = 19 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) imagename = visname + '_'+ ch_peakname + '.ap2.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, #antenna = '!DA50', #DA50 is excluded spw=ch_peak, specmode='cube',start=startch_peak,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') #force model to save modelname=imagename+'.model' ft(vis=vis,model=modelname,usescratch=True) # Image a test channel for sanity checks imagename = visname +'_'+ch_testname +'.ap2.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, #antenna = '!DA50', #DA50 is excluded spw=ch_test, specmode='cube',start=startch_test,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) #Get image statistics for comparison get_im_stats(imagename+'.image') ### [OPTIONAL] IDENTIFY LINE-FREE CHANNELS #--------- mystep = 20 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Image whole cube at lower resolution to identify line-free channels # This step can be omitted if you trust our selection in variable contchans above. for s in ['0','1']: os.system('rm -rf X1de2_VYCMa_325_spw'+s+'_5ap2.clean*') imagename = visname + '_spw'+ s+ '_av5.ap2.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, spw=s, specmode='cube',width=5, cell='0.02arcsec', imsize=512, niter = 1000, mask=cleaningmask, interactive=cleaninteractive) #--------- mystep = 21 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Image continuum imagename = visname + '_cont.ap2.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename=imagename, field=field, spw=contchans, specmode='mfs', cell='0.01arcsec', imsize=512, niter=300, deconvolver = 'multiscale', scales=[0,4,8,12], mask=['ellipse[[250pix,280pix],[70pix,67pix],90deg]','ellipse[[270pix,310pix],[95pix,42pix],20deg]'], interactive=True) exportfits(imagename=imagename+'.image', fitsimage=imagename+'.fits') ### CONTINUUM SUBTRACTION #--------- mystep = 22 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Subtract continuum os.system('rm -rf '+vis+'.contsub') uvcontsub(vis = vis,field='V*',fitorder=1, combine='spw',spw='',fitspw=contchans) # Make image of test channel imagename = visname +'_'+ch_testname +'.contsub.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis, imagename = imagename, field = field, #antenna = '!DA50', #DA50 is excluded spw=ch_test, specmode='cube',start=startch_test,nchan=1, cell='0.01arcsec', imsize=512, niter = 200, mask=cleaningmask, interactive=cleaninteractive) ### ANALYSIS #--------- mystep = 23 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # High resolution data cube for maser # Line frequencies taken from Kaminski et al. 2013 ApJS 209 38 imagename = visname +'.contsub.H2O.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis+'.contsub', imagename= imagename, field = field, spw='0', specmode='cube',start='-38km/s',nchan=134, restfreq='325.152919GHz', cell='0.01arcsec',imsize=512, deconvolver='hogbom', weighting='briggs',robust=0.5, threshold='10mJy', mask=['ellipse[[218pix,255pix],[64pix,52pix],90deg]','ellipse[[281pix,359pix],[140pix,70pix],20deg]'], interactive=cleaninteractive, savemodel = 'modelcolumn', niter=1000) # Make moment map # Quick estimate of RMS accross the cube thisimage=imagename+'.image' box='250,50,500,100' chanstat=imstat(imagename=thisimage,box=box,chans='3') rms1= chanstat['rms'][0] chanstat=imstat(imagename=thisimage,box=box,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, fitsimage=thisimage+'.fits') exportfits(imagename=thisimage+'.mom0', fitsimage=thisimage+'.mom0.fits') #--------- mystep = 24 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) # Image NaCl 24-23 line in spw 1 # Line frequencies taken from Kaminski et al. 2013 ApJS 209 38 imagename = visname +'.contsub.NaCl.clean' os.system('rm -rf '+imagename+'.*') tclean(vis = vis+'.contsub', imagename= imagename, field = field, spw='1', specmode='cube',start='-8km/s',nchan=30, width='2km/s', restfreq='312.1099004GHz', cell='0.01arcsec',imsize=512, deconvolver='hogbom', weighting='briggs',robust=0.5, threshold='4mJy', mask = 'ellipse[[284pix,263pix],[150pix,115pix],20deg]', interactive=cleaninteractive, savemodel = 'modelcolumn', niter=1000) # Make moment map # Quick estimate of RMS accross the cube thisimage=imagename+'.image' box='250,50,500,100' chanstat=imstat(imagename=thisimage,box=box,chans='5') rms1= chanstat['rms'][0] chanstat=imstat(imagename=thisimage,box=box,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, fitsimage=thisimage+'.fits') exportfits(imagename=thisimage+'.mom0', fitsimage=thisimage+'.mom0.fits') exportfits(imagename=thisimage+'.mom0', fitsimage=thisimage+'.mom1.fits') #--------------------------------------------------------------------------------- #----- End of script. #----------------------------------------------------------------------------------