""" This script was written for CASA 5.1.1 Datasets calibrated (in order of date observed): SB1: 2015.1.00964.S Observed 14 June 2016 (2 execution blocks) SB2: 2013.1.00964.S Observed 02 July 2016 (1 execution block) SB3: 2015.1.00484.L Observed 14 May 2017, 17 May 2017, 19 May 2017 (3 execution blocks) LB1: 2016.1.00484.L Observed 26 September 2017 and 26 November 2017 (2 execution blocks) reducer: L. Perez """ """ Starting matter """ import os execfile('reduction_utils.py') skip_plots = True # if True, can run script non-interactively """ Input for loading data """ prefix = 'HD143006' SB1_path = '/full_path/to_calibrated/msfile.ms' SB2_path = '/full_path/to_calibrated/msfile.ms' SB3_path = '/full_path/to_calibrated/msfile.ms' LB1_path = '/full_path/to_calibrated/msfile.ms' # Note that if you are downloading data from the archive, your SPW numbering # may differ from this script, depending on how you split your data out! data_params = {'SB1': {'vis' : SB1_path, 'name' : 'SB1', 'field': 'HD143006', 'line_spws': np.array([0,9]), # CO SPWs 'line_freqs': np.array([2.30538e11,2.30538e11]), 'cont_spws': '0:264~715,1,2,9:264~715,10,11', 'width_array': [1,960,3840,1,960,3840], }, 'SB2': {'vis' : SB2_path, 'name' : 'SB2', 'field': 'HD143006', 'line_spws': np.array([0]), # CO SPWs 'line_freqs': np.array([2.30538e11]), 'cont_spws': '0:264~715,1,2', 'width_array': [1,960,3840], }, 'SB3': {'vis' : SB3_path, 'name' : 'SB3', 'field': 'HD_143006', 'line_spws': np.array([0,4,8]), # CO SPWs 'line_freqs': np.array([2.30538e11,2.30538e11,2.30538e11]), 'cont_spws': '0:1828~2053,1,2,3,4:1828~2053,5,6,7,8:1828~2053,9,10,11', 'width_array': [1,128,128,128,1,128,128,128,1,128,128,128], }, 'LB1': {'vis' : LB1_path, 'name' : 'LB1', 'field' : 'HD_143006', 'line_spws': np.array([3,7]), # CO SPWs 'line_freqs': np.array([2.30538e11, 2.30538e11]), 'cont_spws': '0,1,2,3:1818~2043,4,5,6,7:1818~2043', 'width_array': [128,128,128,1, 128,128,128,1], } } """ Additional flagging """ flagmanager(vis=data_params['SB2']['vis'], mode='save', versionname='original_flags', comment='Original flag states') flagdata(vis=data_params['SB2']['vis'], mode='manual', spw='0:0~200', flagbackup=False, field=data_params['SB2']['field']) """ Split out the data, with the continuum SPWs averaged down """ SB1_CO = prefix+'_SB1_CO' SB2_CO = prefix+'_SB2_CO' SB3_CO = prefix+'_SB3_CO' LB1_CO = prefix+'_LB1_CO' os.system('rm -rf '+SB1_CO+'.ms') split(vis=SB1_path, outputvis=SB1_CO+'.ms', field=data_params['SB1']['field'], spw=data_params['SB1']['cont_spws'], width=data_params['SB1']['width_array'], timebin='0s', datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf '+SB2_CO+'.ms') split(vis=SB2_path, outputvis=SB2_CO+'.ms', field=data_params['SB2']['field'], spw=data_params['SB2']['cont_spws'], width=data_params['SB2']['width_array'], timebin='0s', datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf '+SB3_CO+'.ms') split(vis=SB3_path, outputvis=SB3_CO+'.ms', field=data_params['SB3']['field'], spw=data_params['SB3']['cont_spws'], width=data_params['SB3']['width_array'], timebin='0s', datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) os.system('rm -rf '+LB1_CO+'.ms') split(vis=LB1_path, outputvis=LB1_CO+'.ms', field=data_params['LB1']['field'], spw=data_params['LB1']['cont_spws'], width=data_params['LB1']['width_array'], timebin='6s', datacolumn='data', intent='OBSERVE_TARGET#ON_SOURCE', keepflags=False) """ Apply same shifts and re-scalings as for continuum """ for i in data_params.keys(): shift = prefix+'_'+i+'_lines_shift.ms' os.system('rm -rf '+shift+'*') split(vis=prefix+'_'+i+'_CO.ms', outputvis=shift, datacolumn='data') fixvis(vis=shift, outputvis=shift, field=data_params[i]['field'], phasecenter='J2000 15h58m36.89967s -22d57m15.602730s') #Split all executions to correct for flux scale differences os.system('rm -rf *_lines_exec*') split_all_obs(prefix+'_SB1_lines_shift.ms', prefix+'_SB1_lines_exec') split_all_obs(prefix+'_SB2_lines_shift.ms', prefix+'_SB2_lines_exec') split_all_obs(prefix+'_SB3_lines_shift.ms', prefix+'_SB3_lines_exec') split_all_obs(prefix+'_LB1_lines_shift.ms', prefix+'_LB1_lines_exec') os.system('rm -rf *rescaled.*') rescale_flux(prefix+'_SB1_lines_exec0.ms', [1.062]) rescale_flux(prefix+'_SB1_lines_exec1.ms', [1.023]) rescale_flux(prefix+'_SB2_lines_exec0.ms', [1.002]) rescale_flux(prefix+'_SB3_lines_exec1.ms', [1.038]) rescale_flux(prefix+'_SB3_lines_exec2.ms', [1.073]) rescale_flux(prefix+'_LB1_lines_exec0.ms', [1.080]) rescale_flux(prefix+'_LB1_lines_exec1.ms', [0.961]) """ Apply the self-calibration solutions to the SB data """ SB_CO = prefix+'_SB_CO' os.system('rm -rf %s*' % SB_CO) concat(vis=[prefix+'_SB1_lines_exec0_rescaled.ms', prefix+'_SB1_lines_exec1_rescaled.ms', prefix+'_SB2_lines_exec0_rescaled.ms', prefix+'_SB3_lines_exec0.ms', prefix+'_SB3_lines_exec1_rescaled.ms', prefix+'_SB3_lines_exec2_rescaled.ms'], concatvis=SB_CO+'.ms', dirtol='0.1arcsec', copypointing=False) # First selfcal had individual SPW solutions applycal(vis=SB_CO+'.ms', spw='0~20', gaintable=[prefix+'_SB.p1'], interp='linearPD', calwt=True, flagbackup=False) # Other selfcal had combined SPW solutions spwmap = [0,0,0,3,3,3,6,6,6,9,9,9,9,13,13,13,13,17,17,17,17] applycal(vis=SB_CO+'.ms', spw='0~20', gaintable=[prefix+'_SB.p2', prefix+'_SB.ap'], spwmap=[spwmap]*2, interp='linearPD', calwt=True, flagbackup=False) os.system('rm -rf '+SB_CO+'_selfcal.ms*') split(vis=SB_CO+'.ms', outputvis=SB_CO+'_selfcal.ms', datacolumn='corrected', keepflags=False) # save some space by deleting intermediate files os.system('rm -rf '+SB_CO+'.ms*') os.system('rm -rf '+prefix+'_SB1_lines_exec*ms*') """ Now apply the self-calibration solutions to the combined dataset """ combined_CO = prefix+'_combined_CO' os.system('rm -rf '+combined_CO+'.ms*') concat(vis=[SB_CO+'_selfcal.ms', prefix+'_LB1_lines_exec0_rescaled.ms', prefix+'_LB1_lines_exec1_rescaled.ms'], concatvis=combined_CO+'.ms', dirtol='0.1arcsec', copypointing=False) spwmap = [0,0,0,3,3,3,6,6,6,9,9,9,9,13,13,13,13,17,17,17,17,21,21,21,21,25,25,25,25] applycal(vis=combined_CO+'.ms', spw='0~28', spwmap=[spwmap]*6, gaintable=[prefix+'_combined.p1', prefix+'_combined.p2', prefix+'_combined.p3', prefix+'_combined.p4', prefix+'_combined.p5', prefix+'_combined.ap'], interp='linearPD', calwt=True, applymode='calonly') os.system('rm -rf '+combined_CO+'_selfcal.ms*') split(vis=combined_CO+'.ms', outputvis=combined_CO+'_selfcal.ms', field='', spw='0,3,6,9,13,17,24,28', datacolumn='corrected', keepflags=False) # save some space by deleting intermediate files os.system('rm -rf '+combined_CO+'.ms*') os.system('rm -rf '+prefix+'_LB1_lines_exec*ms*') """ Continuum subtraction """ fitspw = '0:0~150;300~450, 1:0~150;300~450, 2:0~150;300~450, 3:0~75;150~225, 4:0~75;150~225, 5:0~75;150~225, 6:0~75;150~225, 7:0~75;150~225' os.system('rm -rf '+combined_CO+'_selfcal.ms.contsub*') uvcontsub(vis=combined_CO+'_selfcal.ms', spw='0~7', fitspw=fitspw, excludechans=False, solint='int', fitorder=1, want_cont=False) """ Define the channels of interest """ chanstart = '-3km/s' chanwidth = '0.32km/s' nchan = 65 """ Split and regrid into the channels of interest """ # Continuum-subtracted SB_only1 = prefix+'_CO_SBonly1.ms.contsub' os.system('rm -rf '+SB_only1+'*') split(vis=combined_CO+'_selfcal.ms.contsub', outputvis=SB_only1, spw='0,1,2', datacolumn='data') SB_cvel1 = SB_only1+'.cvel' os.system('rm -rf '+SB_cvel1+'*') mstransform(vis=SB_only1, outputvis=SB_cvel1, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') SB_only2 = prefix+'_CO_SBonly2.ms.contsub' os.system('rm -rf '+SB_only2+'*') split(vis=combined_CO+'_selfcal.ms.contsub', outputvis=SB_only2, spw='3,4,5', datacolumn='data') SB_cvel2 = SB_only2+'.cvel' os.system('rm -rf '+SB_cvel2+'*') mstransform(vis=SB_only2, outputvis=SB_cvel2, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') LB_only = prefix+'_CO_LBonly.ms.contsub' os.system('rm -rf '+LB_only+'*') split(vis=combined_CO+'_selfcal.ms.contsub', outputvis=LB_only, spw='6,7', datacolumn='data') LB_cvel = LB_only+'.cvel' os.system('rm -rf '+LB_cvel+'*') mstransform(vis=LB_only, outputvis=LB_cvel, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') contsub_concat = prefix+'_COcube.ms.contsub' os.system('rm -rf '+contsub_concat) concat(vis=[SB_cvel1, SB_cvel2, LB_cvel], concatvis=contsub_concat, dirtol='0.1arcsec', copypointing=False) # With continuum cSB_only1 = prefix+'_CO_SBonly1.ms' os.system('rm -rf '+cSB_only1) split(vis=combined_CO+'_selfcal.ms', outputvis=cSB_only1, spw='0,1,2', datacolumn='data') cSB_cvel1 = cSB_only1+'.cvel' os.system('rm -rf '+cSB_cvel1) mstransform(vis=cSB_only1, outputvis=cSB_cvel1, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') cSB_only2 = prefix+'_CO_SBonly2.ms' os.system('rm -rf '+cSB_only2) split(vis=combined_CO+'_selfcal.ms', outputvis=cSB_only2, spw='3,4,5', datacolumn='data') cSB_cvel2 = cSB_only2+'.cvel' os.system('rm -rf '+cSB_cvel2) mstransform(vis=cSB_only2, outputvis=cSB_cvel2, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') cLB_only = prefix+'_CO_LBonly.ms' os.system('rm -rf '+cLB_only) split(vis=combined_CO+'_selfcal.ms', outputvis=cLB_only, spw='6,7', datacolumn='data') cLB_cvel = cLB_only+'.cvel' os.system('rm -rf '+cLB_cvel) mstransform(vis=cLB_only, outputvis=cLB_cvel, keepflags=False, datacolumn='data', regridms=True, mode='velocity', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz') cont_concat = prefix+'_COcube.ms' os.system('rm -rf '+cont_concat) concat(vis=[cSB_cvel1, cSB_cvel2, cLB_cvel], concatvis=cont_concat, dirtol='0.1arcsec', copypointing=False) """ Imaging (continuum-subtracted only) """ imagename = prefix+'_CO' for ext in ['.image','.mask','.model','.pb','.psf','.residual','.sumwt']: os.system('rm -rf '+ imagename + ext) tclean(vis=prefix+'_COcube.ms.contsub', imagename=imagename, specmode='cube', imsize=2250, deconvolver='multiscale', start=chanstart, width=chanwidth, nchan=nchan, outframe='LSRK', veltype='radio', restfreq='230.538GHz', cell='0.004arcsec', scales = [0,17,50,83], gain=0.3, niter=50000, weighting='briggs', robust=0.8, threshold='2mJy', uvtaper='0.02arcsec', mask='', interactive=True, nterms=1, restoringbeam='common') # Note that masking is done interactively """ Final outputs """ """ Save the final MS files """ os.system('tar cvzf '+prefix+'_CO.ms.tgz '+prefix+'_COcube.ms.contsub') os.system('tar cvzf '+prefix+'_COcont.ms.tgz '+prefix+'_COcube.ms') """ Save the imaged datacube """ exportfits(imagename+'.image', imagename+'.fits')