#!/usr/bin/env """ split, clean, and self-cal continuum and line data NOTE: this is intended to be an interactive, iterative process so this is more a log that should be run by cutting and pasting into casa rather than as an executable script search "CHANGEME" for variables to be changed 10/9/15 MCA """ # ======================== Setup =========================== # LupusIII_50 16:08:30.700 -38:28:26.800 # Class II, K1 field = 9 # CHANGEME file_ms = '../science_calibrated.ms' contspw = '2,3,4' contspw_w = [128,3840,1920] robust = 0.5 imsize = [640,640] cell = '0.03arcsec' imagermode = 'csclean' refant = 'DA52' # CHANGEME xc = 326 # CHANGEME yc = 306 # CHANGEME in_a = 80 out_a = 120 aper = 1.25 boxwidth = 300. box = rg.box([xc-boxwidth,yc-boxwidth],[xc+boxwidth,yc+boxwidth]) # ======================= Split Off Continuum ======================== # split off field from full ms split(vis = file_ms, outputvis = 'f'+str(field)+'.vis', field = field, datacolumn = 'data') # split off continuum (take the large bw spw and average split(vis = 'f'+str(field)+'.vis', outputvis = 'f'+str(field)+'_cont.vis', spw = contspw, width = contspw_w, datacolumn = 'data') # plot uv-distance vs. amplitude plotms(vis='f'+str(field)+'_cont.vis', xaxis='uvdist',yaxis='amp', coloraxis='spw') # plotfile='f'+str(field)+'_ampuv_orig.png' # showgui=False, # highres=True, # overwrite=True) # source is resolved # find antenna close to center of configuration # check pipeline log that this ant is OK plotants(vis='f'+str(field)+'_cont.vis') #, figfile='f'+str(field)+'_ants.png') # ================== Clean continuum before selfcal ================== # os.system('rm -rf f'+str(field)+'_cont_b4sc*') # light clean (100 iterations) to set the mask around the main peaks clean(vis = 'f'+str(field)+'_cont.vis', imagename = 'f'+str(field)+'_cont_b4sc', mode = 'mfs', psfmode = 'clark', niter = 100, threshold = '0.0mJy', interactive = True, mask = '', cell = cell, imsize = imsize, weighting = 'briggs', robust = robust, imagermode = imagermode) im_max = imstat(imagename = 'f'+str(field)+'_cont_b4sc.image')['max'][0] im_rms = imstat(imagename = 'f'+str(field)+'_cont_b4sc.image', region='annulus[['+str(xc)+'pix,'+str(yc)+'pix],['+str(in_a)+'pix,'+str(out_a)+'pix]]')['rms'][0] print 'Peak = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_max, 1000*im_rms, im_max/im_rms) # Peak = 29.86 mJy, rms = 0.41 mJy, S/N = 73.7 # ======================== Self-Calibrate 1 ================== # first combine all the data by time (solint = inf) # i.e., phase self-cal over entire integration time gaincal(vis = 'f'+str(field)+'_cont.vis', caltable = 'f'+str(field)+'_cont_pcal1', refant = refant, solint = 'inf', combine = 'spw', gaintype = 'T', spw = '', calmode = 'p', minblperant = 4, minsnr = 3) # plot phase for each antenna plotcal(caltable = 'f'+str(field)+'_cont_pcal1', xaxis = 'time', yaxis = 'phase', spw = '', iteration = 'antenna', subplot = 421, plotrange = [0,0,-20,20]) # narrow yrange phases close to zero # no variation between integration times (expected since solin=inf) # note DV02 has no data (100% flagged in pipeline calibration) # apply calibration to data applycal(vis = 'f'+str(field)+'_cont.vis', spw = '', gaintable = ['f'+str(field)+'_cont_pcal1'], spwmap = [0,0,0], calwt = T, flagbackup = F) # clean self-calibrated data clean(vis = 'f'+str(field)+'_cont.vis', imagename = 'f'+str(field)+'_cont_pcal1_clean', mode = 'mfs', psfmode = 'clark', niter = 100, threshold = '0.0mJy', interactive = False, mask = 'f'+str(field)+'_cont_b4sc.mask', cell = cell, imsize = imsize, weighting = 'briggs', robust = robust, imagermode = imagermode) im_max = imstat(imagename = 'f'+str(field)+'_cont_pcal1_clean.image')['max'][0] im_rms = imstat(imagename = 'f'+str(field)+'_cont_pcal1_clean.image', region='annulus[['+str(xc)+'pix,'+str(yc)+'pix],['+str(in_a)+'pix,'+str(out_a)+'pix]]')['rms'][0] print 'Peak = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_max, 1000*im_rms, im_max/im_rms) #Peak = 30.31 mJy, rms = 0.40 mJy, S/N = 76.6 (slight improvement) # inspect images imview(raster=[{'file':'f'+str(field)+'_cont_b4sc.image'}, {'file':'f'+str(field)+'_cont_pcal1_clean.image'}]) # second peak has strengthened a little # peaks are slightly more defined # ======================== Self-Calibrate 2 ================== # decrease phase self-cal solution interval to a few times integration time # int = 6s (from X125.log) gaincal(vis = 'f'+str(field)+'_cont.vis', caltable = 'f'+str(field)+'_cont_pcal2', refant = refant, solint = '20s', # CHANGEME combine = 'spw', gaintype = 'T', spw = '', calmode = 'p', minblperant = 4, minsnr = 3) plotcal(caltable = 'f'+str(field)+'_cont_pcal2', xaxis = 'time', yaxis = 'phase', spw = '', iteration = 'antenna', subplot = 421, plotrange = [0,0,-20,20]) applycal(vis = 'f'+str(field)+'_cont.vis', spw = '', gaintable = ['f'+str(field)+'_cont_pcal2'], spwmap = [0,0,0], calwt = T, flagbackup = F) clean(vis = 'f'+str(field)+'_cont.vis', imagename = 'f'+str(field)+'_cont_pcal2_clean', mode = 'mfs', psfmode = 'clark', niter = 100, threshold = '0.0mJy', interactive = False, mask = 'f'+str(field)+'_cont_b4sc.mask', cell = cell, imsize = imsize, weighting = 'briggs', robust = robust, imagermode = imagermode) im_max = imstat(imagename = 'f'+str(field)+'_cont_pcal2_clean.image')['max'][0] im_rms = imstat(imagename = 'f'+str(field)+'_cont_pcal2_clean.image', region='annulus[[320pix,320pix],[60pix,100pix]]')['rms'][0] print 'Peak = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_max, 1000*im_rms, im_max/im_rms) # Peak = 30.32 mJy, rms = 0.39 mJy, S/N = 78.00 (slight change from pcal1) # inspection of the image shows no change from pcal1 imview(raster=[{'file':'f'+str(field)+'_cont_b4sc.image'}, {'file':'f'+str(field)+'_cont_pcal1_clean.image'}, {'file':'f'+str(field)+'_cont_pcal2_clean.image'}]) # very slightly more defined # ======================== Self-Calibrate 3 ================== # try smallest possible solution interval = integration time # 1 of 40 solutions flagged due to SNR < 3 in spw=0 gaincal(vis = 'f'+str(field)+'_cont.vis', caltable = 'f'+str(field)+'_cont_pcal3', refant = refant, solint = 'int', combine = 'spw', gaintype = 'T', spw = '', calmode = 'p', minblperant = 4, minsnr = 3) # can now see the individual integrations within the two visits to target # and there is scatter from one point to the next => intrinsic phase noise plotcal(caltable = 'f'+str(field)+'_cont_pcal3', xaxis = 'time', yaxis = 'phase', spw = '', iteration = 'antenna', subplot = 421, plotrange = [0,0,-20,20]) applycal(vis = 'f'+str(field)+'_cont.vis', spw = '', gaintable = ['f'+str(field)+'_cont_pcal3'], spwmap = [0,0,0], calwt = T, flagbackup = F) clean(vis = 'f'+str(field)+'_cont.vis', imagename = 'f'+str(field)+'_cont_pcal3_clean', mode = 'mfs', psfmode = 'clark', niter = 100, threshold = '0.0mJy', interactive = False, mask = 'f'+str(field)+'_cont_b4sc.mask', cell = cell, imsize = imsize, weighting = 'briggs', robust = robust, imagermode = imagermode) im_max = imstat(imagename = 'f'+str(field)+'_cont_pcal3_clean.image')['max'][0] im_rms = imstat(imagename = 'f'+str(field)+'_cont_pcal3_clean.image', region='annulus[[320pix,320pix],[60pix,100pix]]')['rms'][0] print 'Peak = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_max, 1000*im_rms, im_max/im_rms) # Peak = 30.67 mJy, rms = 0.39 mJy, S/N = 78.6 # inspection of the image shows first peak strengthening again, so will use this one imview(raster=[{'file':'f'+str(field)+'_cont_b4sc.image'}, {'file':'f'+str(field)+'_cont_pcal1_clean.image'}, {'file':'f'+str(field)+'_cont_pcal2_clean.image'}, {'file':'f'+str(field)+'_cont_pcal3_clean.image'}]) # all pcal looks similar so will use first # ======================== Best Continuum Map ================== # so now run the same applycal but with flagbackup = T, applycal(vis = 'f'+str(field)+'_cont.vis', spw = '', gaintable = ['f'+str(field)+'_cont_pcal1'], # CHANGEME spwmap = [0,0,0], calwt = T, flagbackup = T) # deep clean, trying different robust weights clean(vis = 'f'+str(field)+'_cont.vis', imagename = 'f'+str(field)+'_cont_best', mode = 'mfs', psfmode = 'clark', niter = 2000, threshold = '0.0mJy', interactive = True, mask = '', cell = cell, imsize = imsize, weighting = 'briggs', robust = -1, # CHANGEME imagermode = imagermode) # placed mask around outer continuum contour # stopped after 700 iterations once the inside became green im_max = imstat(imagename = 'f'+str(field)+'_cont_best.image')['max'][0] im_rms = imstat(imagename = 'f'+str(field)+'_cont_best.image', region='annulus[['+str(xc)+'pix,'+str(yc)+'pix],['+str(in_a)+'pix,'+str(out_a)+'pix]]')['rms'][0] bmaj = imhead(imagename = 'f'+str(field)+'_cont_best.image', mode="get", hdkey="beammajor") bmin = imhead(imagename = 'f'+str(field)+'_cont_best.image', mode="get", hdkey="beamminor") print 'Peak = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_max, 1000*im_rms, im_max/im_rms) print 'Beam = {0:.2f} x {1:.2f} arcsec'.format(bmaj.get('value'),bmin.get('value')) # robust = -1 # Peak = 29.54 mJy, rms = 0.48 mJy, S/N = 61.7 # Beam = 0.31 x 0.28 arcsec # robust = -1 gives lower SNR but higher resolution # clearer view of the inner cavity so is worth the payoff # save this to a fits file exportfits(imagename='f'+str(field)+'_cont_best.image', fitsimage='f'+str(field)+'_cont.fits') # compare to before self-cal imview(raster=[{'file':'f'+str(field)+'_cont_b4sc.image'}, {'file':'f'+str(field)+'_cont_best.image'}]) # measure flux im_rms = imstat(imagename = 'f'+str(field)+'_cont_best.image', region='annulus[['+str(xc)+'pix,'+str(yc)+'pix],['+str(in_a)+'pix,'+str(out_a)+'pix]]')['rms'][0] im_flux = imstat(imagename = 'f'+str(field)+'_cont_best.image', region='circle[['+str(xc)+'pix,'+str(yc)+'pix],'+str(aper)+'arcsec]')['flux'][0] print 'Flux = {0:.2f} mJy, rms = {1:.2f} mJy, S/N = {2:.1f}'.format(1000*im_flux, 1000*im_rms, im_flux/im_rms) # Flux = 135.14 mJy, rms = 0.48 mJy, S/N = 283.1 # re-center image on source and use measure.py to get COG flux ia.fromimage(outfile = 'f'+str(field)+'_cont_cropped.image', infile = 'f'+str(field)+'_cont.fits', region = box ) ia.close() exportfits(imagename = 'f'+str(field)+'_cont_cropped.image', fitsimage = 'f'+str(field)+'_cont_cropped.fits') ''' Measuring COG for G/f9_cont_cropped.fits Assuming object center (300.0,300.0) Background: 0.00 mJy/beam km/s RMS in annulus 4.0-9.0 arcsec = 0.48 mJy/beam km/s i radius flux err snr (asec) (mJy) (mJy) 0 0.10 5.43 0.15 36.3 1 0.20 20.85 0.37 57.1 2 0.30 46.48 0.61 75.8 3 0.40 73.23 0.62 117.2 4 0.50 98.43 0.71 138.7 5 0.60 115.99 0.75 154.7 6 0.70 126.27 0.97 129.5 7 0.80 130.15 1.14 113.8 8 0.90 131.64 1.09 121.1 9 1.00 133.41 1.03 129.4 10 1.10 134.88 1.37 98.7 11 1.20 135.31 1.12 120.7<--lvels off 12 1.30 135.07 1.42 95.4 13 1.40 135.46 1.67 81.0 14 1.50 136.39 1.66 82.0 15 1.60 137.09 1.52 90.1 F = 137.09 mJy E = 1.52 mJy S = 90.15 D = 3.20 arcsec ''' # ======================== Measure flux with UVMODELFIT ================== # calculate offset from phase center in arcsec pixscale = 0.03 # must match 'cell' dx = pixscale*(320.0-xc) # offset to east (left) dy = pixscale*(yc-320.0) # offset to north (up) # measure flux as gaussian uvmodelfit(vis = 'f'+str(field)+'_cont.vis', comptype = 'G', sourcepar = [im_flux,dx,dy,0.5,0.5,-80.0], varypar = [T,T,T,T,T,T], niter = 10) ''' reduced chi2=1.47494 I = 0.151823 +/- 0.00108974 x = -0.180243 +/- 0.00338909 arcsec y = -0.416658 +/- 0.00149899 arcsec a = 1.14263 +/- 0.00870552 arcsec r = 0.273587 +/- 0.00354017 p = -72.6574 +/- 0.221338 deg consistent with aperture method 16:08:30.685 -38:28:27.217 '''