######################################################################## # Data reduction script for BRI1202 Band 7: # - Imaging script: "BR1202_Band7_Imaging.py" - # # Written for CASA version 3.3 ... For version 3.4 you need to search and # replace CALREADY->USESCRATCH (and possibly other changes). ######################################################################## """ See accompanying README file for details of this script, the necessary input files and comments on the data """ # # ---------------------------------------------------------------------------- # CONTROL PROGRAM FLOW # ---------------------------------------------------------------------------- # Input file (contains only calibrated data for our source) ms_basename = 'uid___A002_X36ffba_X3ef.cal' # Reference antenna (used for selfcal) ref_ant = 'DV06' # Define an empty box (used for noise estimates) empty_box = '10,150,240,240' # Define the clean box clean_box = 'box [ [ 120pix , 120pix] , [140pix, 145pix ] ]' clean_box_simple = '120,120,140,145' # Run CLEAN interactively? Used by all CLEAN runs. interactive_clean = False # -=-=-=-=--=-=-=-=--=-=-=-=- # Step #1: Inspection # -=-=-=-=--=-=-=-=--=-=-=-=- save_orig = True # ... save original flags? do_reset = True # ... reset the flagging to its original state? do_listobs = True # ... list content of observations? # -=-=-=-=--=-=-=-=--=-=-=-=-=-= # Step #2: Continum Subtraction # -=-=-=-=--=-=-=-=--=-=-=-=-=-= do_contsub = True # ... do continuum subtraction do_plot_spw0 = False # ... plot spw0 for line-free channel selection do_uvcontsub = True # ... actually carry out the uv-continuum subtraction do_plot_contfit = False # ... plot the results of the continuum fit # -=-=-=-=--=-=-=-=--=-=-=-=-=-= # Dirty Image # -=-=-=-=--=-=-=-=--=-=-=-=-=-= do_dirtyimage = True # ... make a dirty image dirty_cont = True # ... make a dirty continuum image dirty_line = True # ... make a dirty line image inspect_dirty = False # ... inspect the results? # -=-=-=-=--=-=-=-=--=-=-=-=-=-= # Clean Image # -=-=-=-=--=-=-=-=--=-=-=-=-=-= do_cleanimage = True # ... make a clean image clean_cont = True # ... make a clean continuum image clean_line = True # ... make a clean dirty image inspect_clean = False # ... inspect the results? # -=-=-=-=--=-=-=-=--=-=-=-=-=-= # Self calibration # -=-=-=-=--=-=-=-=--=-=-=-=-=-= do_selfcal = True # ... carry out a self-calibration do_selfcal_line = True # ... self-calibrate the line data do_selfcal_cont = True # ... self-calibrate the continuum do_apply_selfcal = True # ... apply the self-calibration do_image_selfcal = True # ... image the self-calbirated data do_inspect_selfcal = False # ... inspect the results # -=-=-=-=--=-=-=-=--=-=-=-=-=-= # Analysis # -=-=-=-=--=-=-=-=--=-=-=-=-=-= do_stats = True # ... run statistics on the result do_export = True # ... export to fits # ---------------------------------------------------------------------------- # SCRIPT (EXECUTES BASED ON CONTROL FLOW CHOICES ABOVE) # ---------------------------------------------------------------------------- #------------------------------- # FLAGGING AND RESET #------------------------------- if save_orig: print "Saving original flags. Run this once before any selfcal." flagmanager(ms_basename+'.ms', mode='save', versionname='original', comment='Flags before any selfcal.') if do_reset: print "Restoring original flag state and reseting scratch columns." flagmanager(ms_basename+'.ms', mode='restore', versionname='original') clearcal(ms_basename+'.ms') #------------------------------- # RUN LISTOBS TO LOOK AT CONTENT #------------------------------- if do_listobs: print "Running listobs." os.system('rm -rf '+ms_basename+'.listobs.txt') listobs(vis=ms_basename+'.ms', verbose=True, listfile=ms_basename+'.listobs.txt') #------------------------------- # DO CONTINUUM SUBTRACTION #------------------------------- if do_contsub: print "Carrying out continuum subtraction." # Plot up spw 0 to find good channels if do_plot_spw0: plotms(vis=ms_basename+'.ms', xaxis='freq', yaxis='amp', spw='0:44~101', avgtime='1e8', avgscan=True) dummy = raw_input("Hit to continue.") plotms(vis=ms_basename+'.ms', xaxis='freq', yaxis='amp', spw='0:0~43;102~127,1', avgtime='1e8', avgscan=True) # Do uv continuum subtraction if do_uvcontsub: uvcontsub2(vis=ms_basename+'.ms', fitspw='0:0~43;102~127,1', combine='spw', want_cont=True, spw='0~3') # Plot results if do_plot_contfit: plotms(vis=ms_basename+'.ms', xaxis='freq', yaxis='amp', spw='0,1', avgtime='1e8', avgscan=True) dummy = raw_input("Hit to continue.") plotms(vis=ms_basename+'.ms.contsub', xaxis='freq', yaxis='amp', spw='0,1', avgtime='1e8', avgscan=True) dummy = raw_input("Hit to continue.") plotms(vis=ms_basename+'.ms.cont', xaxis='freq', yaxis='amp', spw='0,1', avgtime='1e8', avgscan=True) #------------------------------- # MAKE DIRTY IMAGE #------------------------------- if do_dirtyimage: print "Making dirty images." # Continuum if dirty_cont: os.system('rm -rf BR1202.cont.dirty.*') clean(vis=ms_basename+'.ms', imagename='BR1202.cont.dirty', mode='mfs', calready=True, spw='1,2,3', psfmode='clark', imagermode='csclean', imsize=256, cell=['0.25arcsec'], niter=0, weighting='briggs', robust=0.5, mask=clean_box, interactive=interactive_clean, threshold='0.5mJy') if inspect_dirty: viewer('BR1202.cont.dirty.image') dirty_cont_rms = (imstat('BR1202.cont.dirty.image',box=empty_box))['rms'] print("Continuum dirty image. RMS appears to be "+str(dirty_cont_rms*1e3)+" mJy") dummy = raw_input("Hit to continue.") # Line if dirty_line: os.system('rm -rf BR1202.line.dirty.*') clean(vis=ms_basename+'.ms.contsub', imagename='BR1202.line.dirty', mode='channel', calready=True, spw='0', psfmode='clark', imagermode='csclean', imsize=256, cell=['0.25arcsec'], niter=0, weighting='briggs', robust=0.5, mask=clean_box, interactive=interactive_clean, threshold='2.0mJy', interpolation='nearest', outframe='TOPO', restfreq='333.829GHz') if inspect_dirty: viewer('BR1202.line.dirty.image') dirty_line_rms = (imstat('BR1202.line.dirty.image',box=empty_box))['rms'] print("Line dirty image. RMS appears to be "+str(dirty_line_rms*1e3)+" mJy") dummy = raw_input("Hit to continue.") #------------------------------- # MAKE CLEAN IMAGE #------------------------------- if do_cleanimage: print "Making CLEAN images." # Continuum if clean_cont: os.system('rm -rf BR1202.cont.clean.*') clean(vis=ms_basename+'.ms', imagename='BR1202.cont.clean', mode='mfs', calready=True, spw='1,2,3', psfmode='clark', imagermode='csclean', imsize=256, cell=['0.25arcsec'], niter=10000, weighting='briggs', robust=0.5, mask=clean_box, interactive=interactive_clean, threshold='1.15mJy') if inspect_clean: viewer('BR1202.cont.clean.image') clean_cont_rms = (imstat('BR1202.cont.clean.image',box=empty_box))['rms'] print("Contintuum clean image. RMS appears to be "+ str(clean_cont_rms*1e3)+" mJy") dummy = raw_input("Hit to continue.") # Line if clean_line: os.system('rm -rf BR1202.line.clean.*') clean(vis=ms_basename+'.ms.contsub', imagename='BR1202.line.clean', mode='channel', calready=True, spw='0', psfmode='clark', imagermode='csclean', imsize=256, cell=['0.25arcsec'], niter=10000, weighting='briggs', robust=0.5, mask=clean_box, interactive=interactive_clean, threshold='6.6mJy', interpolation='nearest', outframe='TOPO', restfreq='333.829GHz') if inspect_clean: viewer('BR1202.line.clean.image') clean_line_rms = (imstat('BR1202.line.clean.image',box=empty_box))['rms'] print("Line clean image. RMS appears to be "+str(clean_line_rms*1e3)+" mJy") dummy = raw_input("Hit to continue.") #---------------------------------------- # RUN A SELF CALIBRATION ON THE LINE DATA #---------------------------------------- if do_selfcal: print "Self-calibrating the data." os.system('rm -rf BR1202.cont.clean.*') clean(vis=ms_basename+'.ms', imagename='BR1202.cont.clean', mode='mfs', calready=True, spw='0,1,2,3', psfmode='clark', imagermode='csclean', imsize=256, cell=['0.25arcsec'], niter=10000, weighting='briggs', robust=0.5, mask=clean_box, interactive=interactive_clean, threshold='1.5mJy') os.system('rm -rf phase.selfcal') gaincal(vis=ms_basename+'.ms', caltable='phase.selfcal', calmode='p', spw='0,1,2,3', solint='900s', combine='scan', gaintype='G', refant=ref_ant, minblperant=4, minsnr=2) plotcal(caltable='phase.cont.cal', xaxis='time', yaxis='phase', iteration='antenna', subplot=331, plotrange=[0,0,-180,180]) if do_apply_selfcal: if do_selfcal_line: applycal(vis=ms_basename+'.ms.contsub', gaintable='phase.selfcal', interp='nearest', flagbackup=False) if do_selfcal_cont: applycal(vis=ms_basename+'.ms', gaintable='phase.selfcal', interp='nearest', flagbackup=False) if do_image_selfcal: if do_selfcal_line: os.system('rm -rf BR1202.line.selfcal.clean.*') clean(vis=ms_basename+'.ms.contsub', imagename='BR1202.line.selfcal.clean', mode='channel', calready=True, spw='0', psfmode='clark', imagermode='csclean', imsize=256, cell=['0.25arcsec'], niter=10000, weighting='briggs', robust=0.5, mask=clean_box, interactive=interactive_clean, threshold='6.6mJy', interpolation='nearest', outframe='TOPO', restfreq='333.829GHz') if do_selfcal_cont: os.system('rm -rf BR1202.cont.selfcal.clean.*') clean(vis=ms_basename+'.ms', imagename='BR1202.cont.selfcal.clean', mode='mfs', calready=True, spw='1,2,3', psfmode='clark', imagermode='csclean', imsize=256, cell=['0.25arcsec'], niter=10000, weighting='briggs', robust=0.5, mask=clean_box, interactive=interactive_clean, threshold='0.5mJy') if do_inspect_selfcal: if do_selfcal_line: viewer('BR1202.line.selfcal.clean.image') selfcal_line_rms = \ (imstat('BR1202.line.selfcal.clean.image',box=empty_box))['rms'] print("Line selfcal image. RMS appears to be "+ \ str(selfcal_line_rms*1e3)+" mJy") dummy = raw_input("Hit to continue.") if do_selfcal_cont: viewer('BR1202.cont.selfcal.clean.image') selfcal_cont_rms = \ (imstat('BR1202.cont.selfcal.clean.image',box=empty_box))['rms'] print("Continuum selfcal image. RMS appears to be "+ \ str(selfcal_cont_rms*1e3)+" mJy") dummy = raw_input("Hit to continue.") #--------------------------------------- # Get basic statistics on cleaned images #--------------------------------------- if do_stats: print "Running statistics and making moment maps." print "----------------------------" print "Continuum Image" print "----------------------------" dirty_cont_rms = \ (imstat('BR1202.cont.dirty.image',box=empty_box))['rms'] dirty_cont_flux = \ (imstat('BR1202.cont.dirty.image',box=clean_box_simple))['flux'] clean_cont_rms = \ (imstat('BR1202.cont.clean.image',box=empty_box))['rms'] clean_cont_flux = \ (imstat('BR1202.cont.clean.image',box=clean_box_simple))['flux'] selfcal_cont_rms = \ (imstat('BR1202.cont.selfcal.clean.image',box=empty_box))['rms'] selfcal_cont_flux = \ (imstat('BR1202.cont.selfcal.clean.image',box=clean_box_simple))['flux'] print "RMS in " print "... dirty map "+str(dirty_cont_rms*1e3)+" mJy" print "... clean map "+str(clean_cont_rms*1e3)+" mJy" print "... selfcal map "+str(selfcal_cont_rms*1e3)+" mJy" print "Flux in " print "... dirty map "+str(dirty_cont_flux*1e3)+" mJy" print "... clean map "+str(clean_cont_flux*1e3)+" mJy" print "... selfcal map "+str(selfcal_cont_flux*1e3)+" mJy" print "----------------------------" print "Line Image" print "----------------------------" dirty_line_rms = \ (imstat('BR1202.line.dirty.image',box=empty_box))['rms'] dirty_line_flux = \ (imstat('BR1202.line.dirty.image',box=clean_box_simple))['flux'] clean_line_rms = \ (imstat('BR1202.line.clean.image',box=empty_box))['rms'] clean_line_flux = \ (imstat('BR1202.line.clean.image',box=clean_box_simple))['flux'] selfcal_line_rms = \ (imstat('BR1202.line.selfcal.clean.image',box=empty_box))['rms'] selfcal_line_flux = \ (imstat('BR1202.line.selfcal.clean.image',box=clean_box_simple))['flux'] print "RMS in " print "... dirty map "+str(dirty_line_rms*1e3)+" mJy" print "... clean map "+str(clean_line_rms*1e3)+" mJy" print "... selfcal map "+str(selfcal_line_rms*1e3)+" mJy" print "Flux in " print "... dirty map "+str(dirty_line_flux*1e3)+" mJy*channel" print "... clean map "+str(clean_line_flux*1e3)+" mJy*channel" print "... selfcal map "+str(selfcal_line_flux*1e3)+" mJy*channel" # Use imview to make image of continuum map imview(raster={'file': 'BR1202.cont.dirty.image', 'range': [-0.0009,0.007], 'colorwedge': True}, zoom=1, out='BR1202.cont.dirty.image.png') imview(raster={'file': 'BR1202.cont.clean.image', 'range': [-0.0009,0.007], 'colorwedge': True}, zoom=1, out='BR1202.cont.clean.image.png') imview(raster={'file': 'BR1202.cont.selfcal.clean.image', 'range': [-0.0009,0.007], 'colorwedge': True}, zoom=1, out='BR1202.cont.selfcal.image.png') # Use immoments to make momoment 0 and 1 map of [CII] os.system('rm -rf BR1202.line.dirty.image.mom0') immoments(imagename='BR1202.line.dirty.image', moments=0, chans='44~101', outfile='BR1202.line.dirty.image.mom0', includepix=[3.0*selfcal_line_rms[0], 1e3]) os.system('rm -rf BR1202.line.clean.image.mom0') immoments(imagename='BR1202.line.clean.image', moments=0, chans='44~101', outfile='BR1202.line.clean.image.mom0', includepix=[3.0*selfcal_line_rms[0], 1e3]) os.system('rm -rf BR1202.line.selfcal.image.mom0') immoments(imagename='BR1202.line.selfcal.clean.image', moments=0, chans='44~101', outfile='BR1202.line.selfcal.image.mom0', includepix=[3.0*selfcal_line_rms[0], 1e3]) os.system('rm -rf BR1202.line.selfcal.image.mom1') immoments(imagename='BR1202.line.selfcal.clean.image', moments=1, chans='44~101', outfile='BR1202.line.selfcal.image.mom1', includepix=[5.0*selfcal_line_rms[0], 1e3]) # 3sigma: 44~101 # 2sigma: 52~93 # Use imview to make image of line map imview(raster={'file': 'BR1202.line.dirty.image.mom0', 'range': [0.0068,4.7], 'colorwedge': False}, zoom=2, contour={'file': 'BR1202.cont.selfcal.clean.image', 'levels': [5], 'base': 0, 'unit': 0.25e-3}, out='BR1202.line.dirty.image.png') imview(raster={'file': 'BR1202.line.clean.image.mom0', 'range': [0.0068,4.7], 'colorwedge': False}, zoom=2, contour={'file': 'BR1202.cont.selfcal.clean.image', 'levels': [5], 'base': 0, 'unit': 0.25e-3}, out='BR1202.line.clean.image.png') imview(raster={'file': 'BR1202.line.selfcal.image.mom0', 'range': [0.0068,4.7], 'colorwedge': False}, zoom=2, contour={'file': 'BR1202.cont.selfcal.clean.image', 'levels': [5], 'base': 0, 'unit': 0.25e-3}, out='BR1202.line.selfcal.image.png') #--------------------------------------- # Export FITS Files #--------------------------------------- if do_export: print "Exporting FITS files." exportfits(imagename='BR1202.line.selfcal.clean.image', fitsimage='BR1202.line.fits', overwrite=True) exportfits(imagename='BR1202.line.selfcal.image.mom0', fitsimage='BR1202.line.mom0.fits', overwrite=True) exportfits(imagename='BR1202.cont.selfcal.clean.image', fitsimage='BR1202.cont.fits', overwrite=True) #---------------------------------------------------------------------------------- #----- End of imaging script. #----------------------------------------------------------------------------------