### For running a selection of a few steps use: #my_steps=[0,1,2,3,4,5,6,7,8,9,10] ### For running a row of many steps use: #my_steps=range(0,10) # from 0 to 9 my_vis = 'uid___A002_X8a5fcf_X125f.ms.split.cal' # ms-file name targetname = 'ngc_1614' # choose name for image praefix fld=targetname if 0 in my_steps : ######################################## ### Inspection for spectral line ### ######################################## ## Visually inspect amplitudes vs. frequency of our scientific target in order to identify ## emission line feature(s) os.system('mkdir '+my_vis+'.split.inspect.plots') # make new folder for plots plotms(vis=my_vis, # ms name spw='', xaxis='channel', yaxis='amp', field=fld, # target avgtime='1e8', avgscan=True, showatm=True, # Show the atmospheric transmission? iteraxis='spw', # plot for each SpW (iteration) - redo for iteration over scan (set also avgscan=False!) coloraxis='corr', # colour code for data points showgui=True, # open user interface plotfile=my_vis+'.split.inspect.plots/field5-amp-channel.png', # export the plot - iteration name will be inserted automatically - change field name manually! exprange='all', # make plot for every iteration overwrite=True) ## user_check=raw_input('When you are done with the window, close it and press enter to continue:') ## CO line: spw 0:850 - 1200 ## unknown: spw1&2 - noise? ## Spectral window 3 is devoid of any line emission -> pure continuum; let's check if continuum ## source is extended plotms(vis=my_vis, # ms name spw='', xaxis='uvwave', # uvdist[m]/wavelength yaxis='amp', field=fld, # target avgtime='1e8', avgchannel='1e8', coloraxis='baseline', iteraxis='spw', # plot for each SpW (iteration) - redo for iteration over scan (set also avgscan=False!) #coloraxis='corr', # colour code for data points showgui=True, # open user interface plotfile=my_vis+'.split.inspect.plots/field5-amp-uvdist.png', # export the plot - iteration name will be inserted automatically - change field name manually! exprange='all', # make plot for every iteration overwrite=True) ## In addition, note the largest baseline in the plot in klambda ## -> corresponds to an angular resolution of theta[arcsec] ~ 206/baseline[klambda] ## ## b = 352 klambda ---> theta = 0.59'' (spw3) ## ## We will need this information to set up the pixel scale of the image. ## user_check=raw_input('When you are done with the window, close it and press enter to continue:') if 1 in my_steps : #################################### ### IMAGING: DIRTY continuum ### #################################### ###################### dirty continuum map ################################ ## ## Let's obtain a dirty map of the continuum emission, excluding all line emission and ## avoiding higher noise channels os.system('rm -rf '+targetname+'_cont_dirty.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_dirty', # image name spw='0:0~850;1200~1919,1,2,3', # channel selection - avoid spectral features and noisy channels specmode = 'mfs', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=512, # number of pixels per image axis - depends primary beam size (field-of-view) - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter interactive=False) # work in a GUI or not? ## logger: Fitted beam used in restoration: 0.808463 by 0.486767 (arcsec) at pa -86.7345 (deg) ## ## Use the viewer to inspect the dirty continuum image and find an emission-free region viewer() if 2 in my_steps : #################################### ### IMAGING: CLEAN continuum ### #################################### ############## Obtain the statistics for the dirty continuum map ############# ## ## Get an idea of the rms-noise in an emission-free corner of the image ## ## Viewer: Use the region marker to select a part of your image (box, ellipse, polygon), ## and look at 'regions panel' (might have to activate it in 'view -> regions') ## ## rms from viewer: 0.74mJy ## ## ## - or - ## Script/Terminal: A more reproducible approach: stats_cont_dirty = imstat(imagename=targetname+'_cont_dirty.image', # image name box='10,10,100,100') # define a box-shaped region ('x_min,y_min,x_max,y_max') where to obtain the statistics from ## imstat, like many analysis tasks in CASA, creates a dictionary of parameters - skim through it to get an impression of it (beware units): print('\n Dictionary of the statistics of the dirty continuum image: \n\n', stats_cont_dirty) # \n means 'new line' ## Get the rms-noise of dirty continuum image in mJy/beam: ## Select the first entry of the 'rms'-list rms_cont_dirty=stats_cont_dirty["rms"][0]*1e3 ## define a cleaning threshold = maximum remaining residual when cleaning clthr_cont=3*rms_cont_dirty # a source is called 'detected' when its flux > 3*sigma - still noise peaks till 5*sigma possible - use conservative limit for extended emission! clthr_cont="{:.1f}".format(clthr_cont)+'mJy' # reformat the variable entry into a CASA-digestible shape ###################### clean continuum map ################################ ## ## Let's obtain a clean map of the continuum emission, excluding all line emission and ## avoiding higher noise channels. We will use a rough clean mask os.system('rm -rf '+targetname+'_cont.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont', # image name spw='0:0~850;1200~1919,1,2,3', # channel selection - avoid spectral features and noisy channels field=targetname, # target specmode = 'mfs', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% niter=10000, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold=clthr_cont, # clean threshold - maximum remaining residual when cleaning mask='box[[245pix,245pix],[277pix,277pix]]', # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=512, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter interactive=False) # work in a GUI or not? ############## Obtain the statistics for the clean continuum map ############# ## ## Get the rms: stats_cont = imstat(imagename=targetname+'_cont.image', box='10,10,100,100') # emission-free region rms_cont=stats_cont["rms"][0] # rms-noise of continuum image in mJy/beam ## Get the peak-flux in the image: stats_cont = imstat(imagename=targetname+'_cont.image', box='245,245,277,277') # reuse the clean mask is sufficient peak_cont=stats_cont["max"][0] # peak-flux of continuum image in mJy/beam ## Derive the dynamic range: dyn_cont=peak_cont/rms_cont # dynamic range in continuum image print('\nStatistics of the clean continuum image: \nRMS:', rms_cont*1e3, 'mJy - peak:', peak_cont*1e3, 'mJy - dynamic range:', dyn_cont) # see your requested statistics in the terminal ## Output a png-image of the clean continuum maps imview(raster={'file': targetname+'_cont.image', # image name 'colorwedge':True, # color bar 'range':[-rms_cont, peak_cont], # flux values from ... to ... 'scaling':0, # brightness scaling 'colormap':'Rainbow 2'}, # map color scheme out=targetname+'_cont.png', # output image zoom=3) # zoom factor if 3 in my_steps : ###################### dirty continuum map at different uv weighting (natural) ################################ ## ## Let's obtain a lower resolution dirty map of the continuum emission (using close to natural ## weighting, excluding all line emission and avoiding higher noise channels. os.system('rm -rf '+targetname+'_cont_dirty_nw*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_dirty_nw', # image name spw='0:0~850;1200~1919,1,2,3', # channel selection - avoid spectral features and noisy channels specmode = 'mfs', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=512, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=2.0, # briggs weighting scheme parameter to 'natural' interactive=False) # work in a GUI or not? ############## Obtain the statistics for the newly weighted dirty continuum map ############# ## ## Get an idea of the rms-noise in an emission-free corner of the image stats_cont_dirty_nw = imstat(imagename=targetname+'_cont_dirty_nw.image', box='10,10,100,100') rms_cont_dirty_nw=stats_cont_dirty_nw["rms"][0]*1e3 # rms-noise of dirty continuum image in mJy/beam clthr_cont_nw=3*rms_cont_dirty_nw # a source is called 'detected' when its flux > 3*sigma - still noise peaks till 5*sigma possible - use conservative limit for extended emission! clthr_cont_nw="{:.1f}".format(clthr_cont_nw)+'mJy' ###################### clean continuum map at different uv weighting (natural)################################ ## ## Let's obtain a lower resolution clean map of the continuum emission, excluding all line emission and ## avoiding higher noise channels. We will use a rough clean mask. We will use a robust parameter close ## to natural weighting os.system('rm -rf '+targetname+'_cont_nw.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_nw', # image name spw='0:0~850;1200~1919,1,2,3', # channel selection - avoid spectral features and noisy channels specmode = 'mfs', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=500, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold=clthr_cont_nw, # clean threshold - maximum remaining residual when cleaning mask='box[[245pix,245pix],[277pix,277pix]]', # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=512, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=2.0, # weighting scheme parameter to 'natural' interactive=False) # work in a GUI or not? ############## Obtain the statistics for the newly weighted clean continuum map ############# ## ## Get the rms: stats_cont_nw = imstat(imagename=targetname+'_cont_nw.image', box='10,10,100,100') rms_cont_nw=stats_cont_nw["rms"][0] # rms-noise of continuum image in mJy/beam ## Get the peak-flux in the image: stats_cont_nw = imstat(imagename=targetname+'_cont_nw.image', box='245,245,277,277') peak_cont_nw=stats_cont_nw["max"][0] # peak-flux of continuum image in mJy/beam ## Derive the dynamic range: dyn_cont_nw=peak_cont_nw/rms_cont_nw # dynamic range in continuum image print('\nStatistics of the naturally weighted clean continuum image: \nRMS:', rms_cont_nw*1e3, 'mJy - peak:', peak_cont_nw*1e3, 'mJy - dynamic range:', dyn_cont_nw) # see your requested statistics in the termianl ## Output a png-image of the clean continuum maps imview(raster={'file': targetname+'_cont_nw.image', # image name 'colorwedge':True, # color bar 'range':[-rms_cont_nw, peak_cont_nw], # flux values from ... to ... 'scaling':0, # brightness scaling 'colormap':'Rainbow 2'}, # map color scheme out=targetname+'_cont_nw.png', # output image zoom=3) # zoom factor ###################### dirty continuum map at different uv weighting (uniform) ################################ ## ## Let's obtain a lower resolution dirty map of the continuum emission (using close to uniform ## weighting, excluding all line emission and avoiding higher noise channels. os.system('rm -rf '+targetname+'_cont_dirty_uw*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_dirty_uw', # image name spw='0:0~850;1200~1919,1,2,3', # channel selection - avoid spectral features and noisy channels specmode = 'mfs', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=512, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=-2.0, # briggs weighting scheme parameter to 'natural' interactive=False) # work in a GUI or not? ############## Obtain the statistics for the newly weighted dirty continuum map ############# ## ## Get an idea of the rms-noise in an emission-free corner of the image stats_cont_dirty_uw = imstat(imagename=targetname+'_cont_dirty_uw.image', box='10,10,100,100') rms_cont_dirty_uw=stats_cont_dirty_uw["rms"][0]*1e3 # rms-noise of dirty continuum image in mJy/beam clthr_cont_uw=3*rms_cont_dirty_uw # a source is called 'detected' when its flux > 3*sigma - still noise peaks till 5*sigma possible - use conservative limit for extended emission! clthr_cont_uw="{:.1f}".format(clthr_cont_uw)+'mJy' ###################### clean continuum map at different uv weighting (uniform) ################################ ## ## Let's obtain a lower resolution clean map of the continuum emission, excluding all line emission and ## avoiding higher noise channels. We will use a rough clean mask. We will use a robust parameter close ## to uniform weighting os.system('rm -rf '+targetname+'_cont_uw.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_uw', # image name spw='0:0~850;1200~1919,1,2,3', # channel selection - avoid spectral features and noisy channels specmode = 'mfs', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=500, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold=clthr_cont_uw, # clean threshold - maximum remaining residual when cleaning mask='box[[245pix,245pix],[277pix,277pix]]', # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=512, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=-2.0, # weighting scheme parameter to 'natural' interactive=False) # work in a GUI or not? ############## Obtain the statistics for the newly weighted clean continuum map ############# ## ## Get the rms: stats_cont_uw = imstat(imagename=targetname+'_cont_uw.image', box='10,10,100,100') rms_cont_uw=stats_cont_uw["rms"][0] # rms-noise of continuum image in mJy/beam ## Get the peak-flux in the image: stats_cont_uw = imstat(imagename=targetname+'_cont_uw.image', box='245,245,277,277') peak_cont_uw=stats_cont_uw["max"][0] # peak-flux of continuum image in mJy/beam ## Derive the dynamic range: dyn_cont_uw=peak_cont_uw/rms_cont_uw # dynamic range in continuum image print('\nStatistics of the uniformly weighted clean continuum image: \nRMS:', rms_cont_uw*1e3, 'mJy - peak:', peak_cont_uw*1e3, 'mJy - dynamic range:', dyn_cont_uw) # see your requested statistics in the termianl ## Output a png-image of the clean continuum maps imview(raster={'file': targetname+'_cont_uw.image', # image name 'colorwedge':True, # color bar 'range':[-rms_cont_uw, peak_cont_uw], # flux values from ... to ... 'scaling':0, # brightness scaling 'colormap':'Rainbow 2'}, # map color scheme out=targetname+'_cont_uw.png', # output image zoom=3) # zoom factor if 4 in my_steps : ########################################################## ### LINE IMAGING Preparation: Continuum Subtraction ### ########################################################## #Subtract continuum emission from line channels in uv-plane os.system('rm -rf '+my_vis+'.contsub') # delete old file to start from scratch, else clean will continue to work on this image ! uvcontsub(vis=my_vis, # ms name fitspw='xx', # line-free channels field=targetname, # target solint ='xx', # solution interval - best: finest fitorder = xx, # 0 for constant, 1 for spectral slope, 2 for curvature (not recommended - SNR ususlly not high enough) combine='spw') # ## Take care, the new MS contsub has its SpWs and Fields re-numbered ## In doubt LISTOBS !! if 5 in my_steps : ##################################### ### IMAGING: line - dirty image ### ##################################### ###################### dirty line cube ###################################### ## ## Let's make a dirty cube of the line emission first. Think about the velocity ## range and channel size you want to image. Have the rest frequency at the target's ## redshift ready. os.system('rm -rf '+targetname+'_line_CO_dirty.*') # delete old file to start from scratch, else clean will continue to work on this image ! os.system('rm -rf '+targetname+'_spw*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_line_CO_dirty', # image name field = targetname, # NGC1614 # target spw = 'xx', # channel selection - the spectral feature you want to image specmode = 'xx', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = xx, # number of channels start = 'xxkm/s', # starting channel width = 'xxkm/s', # channel width restfreq='xxGHz', # CO(1-0) Obs frequency (z=xx) # rest frequency niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=xx, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='xxarcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=xx, # briggs weighting scheme parameter to something intermediate interactive=False) # work in a GUI or not? ##### check for lines in all spw tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_spw0_dirty', # image name field = targetname, # NGC1614 # target spw = '0', # channel selection - the spectral feature you want to image specmode = 'cube', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = -1, # number of channels start = '', # starting channel width = '20km/s', # channel width niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=512, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter to something intermediate interactive=False) # work in a GUI or not? tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_spw1_dirty', # image name field = targetname, # NGC1614 # target spw = '1', # channel selection - the spectral feature you want to image specmode = 'cube', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = -1, # number of channels start = '', # starting channel width = '20km/s', # channel width niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=512, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter to something intermediate interactive=False) # work in a GUI or not? tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_spw2_dirty', # image name field = targetname, # NGC1614 # target spw = '2', # channel selection - the spectral feature you want to image specmode = 'cube', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = -1, # number of channels start = '', # starting channel width = '20km/s', # channel width niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=512, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter to something intermediate interactive=False) # work in a GUI or not? tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_spw3_dirty', # image name field = targetname, # NGC1614 # target spw = '3', # channel selection - the spectral feature you want to image specmode = 'cube', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = -1, # number of channels start = '', # starting channel width = '20km/s', # channel width niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=512, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter to something intermediate interactive=False) # work in a GUI or not? if 6 in my_steps : ########################################## ### IMAGING: CO line - CLEAN image ### ########################################## ## Have a look at it - determine a set of line free channels and check the noise # viewer(targetname+'_line_CO_dirty.image') ## moment 6 gives average rms for a given number of (line-free) channels os.system('rm -rf '+targetname+'_line_CO_dirty.mom6') # delete old immoments(imagename=targetname+'_line_CO_dirty.image', # image name moments=[6], # 'CASA moment' to use chans='xx', # line-free channels to use outfile=targetname+'_line_CO_dirty.mom6') # output name ## get the image statistics of the moment 6 map in a definitely emission free corner (like for the continuum) imst=imstat(imagename=targetname+'_line_CO_dirty.mom6', # input name region='box[[xxpix,xxpix],[xxpix,xxpix]]') # area in corner of the image ## Construct an imaging region mask: set pixels with a residual above a defined threshold to 1, below to 0. ## Depending on the noise and the PSF quality you might have to start with a high threshold to really mark only the target source and no noise peaks os.system('rm -rf '+targetname+'_line_CO_10SN_region') # delete old immath(imagename = targetname+'_line_CO_dirty.residual', # image name outfile = targetname+'_line_CO_10SN_region', # output expr = 'iif(IM0>'+str(imst["rms"][0]*10)+',1,0)') # mathematical operation - here: if pixel value in input > certain flux (here: multiple of rms), set output pixel to 1, else 0 ###################### clean line cube (part 1) ###################################### ## ## Let's start working on the clean cube of the line emission. For this, add the region ## mask (generated above) and a threshold to the parameter list. Depending on the messiness ## of the source morphology, the noise level and the PSF quality you could use a slightly ## lower threshold than used for generating the mask. os.system('rm -rf '+targetname+'_line_CO.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_line_CO', # image name field = targetname, # NGC1614 # target spw = '0', # channel selection - the spectral feature you want to image specmode = 'cube', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = 50, # number of channels start = '-500km/s', # starting channel width = '20km/s', # channel width restfreq='113.462GHz', # CO(1-0) Obs frequency (z=0.01594) # rest frequency niter=10000, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold = str(imst["rms"][0]*7)+'Jy', # clean threshold - maximum remaining residual when cleaning mask=targetname+'_line_CO_10SN_region', # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=512, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter to something intermediate interactive=False) # work in a GUI or not? ## moment 6 gives average rms for a given number of (line-free) channels os.system('rm -rf '+targetname+'_line_CO.mom6') # delete old immoments(imagename=targetname+'_line_CO.image', # image name moments=[6], # 'CASA moment' to use chans='xx', # line-free channels to use outfile=targetname+'_line_CO.mom6') # output name ## get the image statistics of the moment 6 map in a definitely emission free corner (like for the continuum) imst=imstat(imagename=targetname+'_line_CO.mom6', # input name region='box[[xxpix,xxpix],[xxpix,xxpix]]') # area in corner of the image ## For a careful approach you could uncomment the following sub-steps and redo the line imaging step. ## For now, we skip that and jump to the interactive mode below. ##################### Beginning of finer clean steps section #################### ## Construct an imaging region mask: set pixels with a residual above a defined threshold to 1, below to 0. ## Depending on the noise and the PSF quality you might have to start with a high threshold to really mark only the target source and no noise peaks os.system('rm -rf '+targetname+'_line_CO_7SN_region') # delete old immath(imagename = targetname+'_line_CO.residual', # image name outfile = targetname+'_line_CO_7SN_region', # output expr = 'iif(IM0>'+str(imst["rms"][0]*7.0)+',1,0)') # mathematical operation - here: if pixel value in input > certain flux, set output pixel to 1, else 0 ###################### clean line cube (part 2) ###################################### ## ## Continue with additional mask in the mask list parameter. Depending on the messiness ## of the source morphology, the noise level and the PSF quality you could use a slightly ## lower threshold than used for generating the mask. tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_line_CO', # image name field = targetname, # NGC1614 # target spw = '0', # channel selection - the spectral feature you want to image specmode = 'cube', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = 50, # number of channels start = '-500km/s', # starting channel width = '20km/s', # channel width restfreq='113.462GHz', # CO(1-0) Obs frequency (z=0.015938) # rest frequency niter=1000000, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold = str(imst["rms"][0]*5.0)+'Jy', # clean threshold - maximum remaining residual when cleaning mask=[targetname+'_line_CO_10SN_region', targetname+'_line_CO_7SN_region'], # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=512, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter to something intermediate interactive=False) # work in a GUI or not? ## moment 6 gives average rms for a given number of (line-free) channels os.system('rm -rf '+targetname+'_line_CO.mom6') # delete old immoments(imagename=targetname+'_line_CO.image', # image name moments=[6], # 'CASA moment' to use chans='xx', # line-free channels to use outfile=targetname+'_line_CO.mom6') # output name ## get the image statistics of the moment 6 map in a definitely emission free corner (like for the continuum) imst=imstat(imagename=targetname+'_line_CO.mom6', # input name region='box[[xxpix,xxpix],[xxpix,xxpix]]') # area in corner of the image ## Construct an imaging region mask: set pixels with a residual above a defined threshold to 1, below to 0. ## Depending on the noise and the PSF quality you might have to start with a high threshold to really mark only the target source and no noise peaks os.system('rm -rf '+targetname+'_line_CO_5SN_region') # delete old immath(imagename = targetname+'_line_CO.residual', # image name outfile = targetname+'_line_CO_5SN_region', # output expr = 'iif(IM0>'+str(imst["rms"][0]*5.0)+',1,0)') # mathematical operation - here: if pixel value in input > certain flux, set output pixel to 1, else 0 ###################### clean line cube (part 3) ###################################### ## ## Continue with additional mask in the mask list parameter. Depending on the messiness ## of the source morphology, the noise level and the PSF quality you could use a slightly ## lower threshold than used for generating the mask. tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_line_CO', # image name field = targetname, # NGC1614 # target spw = '0', # channel selection - the spectral feature you want to image specmode = 'cube', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = 50, # number of channels start = '-500km/s', # starting channel width = '20km/s', # channel width restfreq='113.462GHz', # CO(1-0) Obs frequency (z=0.015938) # rest frequency niter=1000000, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold = str(imst["rms"][0]*4.0)+'Jy', # clean threshold - maximum remaining residual when cleaning mask=[targetname+'_line_CO_10SN_region', targetname+'_line_CO_7SN_region', targetname+'_line_CO_5SN_region'], # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=512, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter to something intermediate interactive=False) # work in a GUI or not? ## ok, we are hitting a limit - residual is too noisy for another flux clipping mask - time for interactive clean ##################### End of finer clean steps section #################### ###################### clean line cube (part 4) ###################################### ## ## Call clean with the same parameter setup except from activating the interactive mode. ## Set the clean mask manually now. ## Here we do not erase the previsouly created targetname+'_line_CO' file. ## So CLEAN will restart where it finished before, i.e., from the residual map ! tclean(vis=my_vis+'.contsub', # ms name imagename=targetname+'_line_CO', # image name field = targetname, # NGC1614 # target spw = '0', # channel selection - the spectral feature you want to image specmode = 'cube', # continuum or spectral line modes deconvolver = 'hogbom', # minor cycle clean algorithm #outframe='LSRK', # reference frame nchan = 50, # number of channels start = '-500km/s', # starting channel width = '20km/s', # channel width restfreq='113.462GHz', # CO(1-0) Obs frequency (z=0.015938) # rest frequency niter=1000000, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold = str(imst["rms"][0]*3.0)+'Jy', # clean threshold - maximum remaining residual when cleaning #mask=[targetname+'_line_CO_10SN_region'], # clean mask; Since it is a continuation, targetname+'_line_CO'.mask already contains the mask used previously, i.e., '10SN_regions', so no needs to add it here imsize=512, # number of pixels per image axis - depends on source and primary beam size - try it out or estimate cell='0.1arcsec', # size of a pixel - typically a beam axis (use e.g. baseline determined above) should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter to something intermediate interactive=True) # work in a GUI or not? ## If you are sure about a source detection (experience and/or pre-knowledge from other ## observations) you can go below 3 sigma (threshold) to improve the signal from these sources ## - but it's risky: the closer to the noise level, the more noise peaks clean is likely to ## declare as 'real'. Benefit at question, because: image = model (*) PSF + residual. ## So, what is left below the 3 sigma will not be thrown away but be accounted for in your ## final image. ## ## Note the pedestal in the residual corresponding to the source's extent - typical ## (delta-function) clean issue. ## ## If you like your mask, save a copy of the *.mask file. When you start a new clean, ## you can replace its *.mask file by your backup (only if same parameters for dimensions of the image). ##