import numpy as np import os #====== FUNCTION DEFINITION ====# def getTargetSciSPW(thisMS): """ Function to get the science spectral windows (so ignoring #CH_AVG """ """ SPWs) for ALMA data. """ """ Returns science spectral window IDs as a comma separated string """ """ ready for use in e.g. split """ """ To be run in CASA """ """ Copy this function definition into your CASA script and run as """ """ mySPWstring = getTargetSciSPW(myMS) """ msmd.open(thisMS) targSPWs = msmd.spwsforintent('OBSERVE_TARGET#ON_SOURCE') chavgSPWs = msmd.almaspws(chavg=True,wvr=True) msmd.close() msmd.done() sciSPWstring = ','.join(map(str,np.setdiff1d(targSPWs,chavgSPWs))) return sciSPWstring #===== END FUNCTION DEFINITION ====# #===== MAIN CODE ====# #--- Initial parameters ---# myMS = 'uid___A002_Xd98580_X354.ms.split.cal' #--- This is a final product produced by re-running script for PI target = 'Z_CMa' #--- This is the target name as listed in the data. #--- Parameters for tutees to fill in ---# contSPW ='0,1,2,3,4,5,6,7,8:300~959,9:0~85;115~959,10' myCell = '0.152arcsec' #--- cell/pixel size of your images. A string of format 'X.xxxarcsec' myImsize = 320 myThreshold = '0.226mJy' #--- 3x your calculated theoretical threshold. A string of format 'Y.ymJy' nChanCO = 160 #--- Number of channels to image in CO cube. Format integer startCO = 40 #--- Start Channel of CO cube. Format integer myThreshCO = '32.0mJy' #--- CO channel noise threshold. A string of format 'Y.ymJy' """nChanCS = '' #--- Number of channels to image in 13CS cube startCS = '' #--- Start Channel of CO cube . Format integer myThreshCS = '' #--- 13CS channel noise threshold. A string of format 'Y.ymJy """ myChansRed = '' #--- Redshifted emission channel range. A string of format 'X~Y' myChansBlue = '' #--- Redshifted emission channel range. A string of format 'X~Y' #--- Logic switches ---# """ Set the step you wish to run to = True and all others to = False """ step0 = False #--- Run a listobs on the measurementSet myMS step1 = False #--- Find and split out science Spectral Windows (SPWs) step2 = False #--- Inpsect the data inPlotsMS step3 = False #--- Image the continuum. step4 = True #--- Image the CO and 13^CS line step5 = False #--- Moment maps step6 = False #--- Make nicer momments #-------------------# #---- THE STEPS ----# #-------------------------------------------------# #--- STEP 0: Make listobs of the calibrated MS ---# if step0: print "\n!!!STEP 0: Make listobs of the calibrated measurement set (MS).!!!\n" """ This step makes a listobs file of our calibrated MS. The listobs """ """ file contains useful information about the observation setup e.g. """ """ frequency settings, time on source, number of antennas and more. """ """ Making a listob at this stages means we can see what data in this """ """ MS we wish to retain. """ if os.path.exists(myMS): print '>> Generating listobs file of '+myMS listobs(myMS,listfile=myMS+'.listobs') """>> You can now inspect the listobs file in a standard text editor. """ else: print ">> You do not seem to have "+myMS+" in the current directory, " print ">> skipping the rest of this step. """ #---------------------------------------------------------# #--- STEP 1: Split out science target and science SPWs ---# if step1: print "\n!!! STEP 1: Split out science target and science Spectral Windows (SPWs).!!!\n" """ myMS contains lots of data we won't actually need during imaging the""" """ science target. For example scans on the observed calibrators and """ """ spectral windows (SPWs) which are used only for calibration. """ """ To give us a smaller MS to work with we will first split out from """ """ myMS, the Target source and the science SPWs. """ scispws = getTargetSciSPW(myMS) print '>> scispws: '+scispws split(vis=myMS, intent = 'OBSERVE_TARGET#ON_SOURCE', datacolumn='corrected', field=target, spw = scispws, outputvis=target+'_TM2.split.cal') print '>> New measurement set: '+target+'_TM2.split.cal' print '>> Generating listobs file of new MS' """ Now we generate a listobs file to see the contents of our new """ """ measurement set """ listobs(target+'_TM2.split.cal',listfile=target+'_TM2.split.cal.listobs') #------------------------------------------# #--- STEP 2: Inspect the data in plotMS ---# if step2: print "\n!!! STEP 2: Inspect the data in plotMS !!!\n" print ">> This is a hands-on step, please see the instructions on the tutorial webpage." print ">> Once you have inspected the data please fill in the contSPW parameter in" print ">> the 'Initial parameters' section of this script." #-------------------------------# #--- STEP 3: Image continuum ---# if step3: print "\n!!!STEP 3: Image continuum!!!\n" print ">> Here we will image the spectral line free continuum emission of " print ">> the target source Z-CMa.\n" print ">> Please read STEP3 of the tutorial, calculate and input the required" print ">> parameters into the 'Initial parameters' section of this script.\n" print ">> Before proceeding, please confirm you have set contSPW, myCell," print ">> myImsize and myThreshold? [y/n]" answerS3 = raw_input() if answerS3 == 'y': doContClean = True else: print " >> Please fill in these parameters before proceeding." doContClean = False if doContClean: print " >> Cleaning Z-CMa continuum." tclean(vis = target+'_TM2.split.cal', imagename = target+'_ALL_SPW_CONTINUUM', field = target, spw=contSPW, threshold=myThreshold, cell=[myCell], imsize=[myImsize,myImsize], #should cover about the FoV outframe='LSRK', niter=10000, deconvolver='multiscale', scales=[0,5,15],# inpixels equivalent to point source and 1, 3 times the beam. weighting = 'briggs', robust=0.5, pbcor=True, specmode='mfs', restoringbeam='common', nterms=1, chanchunks=-1, gridder='standard', interactive=True ) #----------------------------------------------------# #--- STEP 4: Subtract the continuum & Image lines ---# if step4: print "\n!!!STEP 4: Subtract the continuum & Image lines!!!\n" """ First we need to generate a continuum subtracted measurement set. """ """ For this we use the same SPW & channel ranges from the continuum """ """ imaging to give CASA a range of channels to fit the continuum to """ """ within the task uvcontsub""" print ">> Generating continuum subtracted measurement set." uvcontsub(vis = target+'_TM2.split.cal', fitspw = contSPW) print ">> "+target+'_TM2.split.cal now created.' """ This creates a new MS called target+'_TM2.split.cal.contsub which we""" """ should use for line imaging """ """ Now we image the CO line, we will use CASA's automasking capability.""" """ CLEANING OF THE CO CUBE """ #---uncommeny below lines if you want to delete existing image files ---# #--- with same file name before imaging ---# #os.system('rm -rf '+target+'_CO_cube.image') #os.system('rm -rf '+target+'_CO_cube.image.pbcor') #os.system('rm -rf '+target+'_CO_cube.psf') #os.system('rm -rf '+target+'_CO_cube.mask') #os.system('rm -rf '+target+'_CO_cube.model') #os.system('rm -rf '+target+'_CO_cube.residual') #os.system('rm -rf '+target+'_CO_cube.sumwt') #os.system('rm -rf '+target+'_CO_cube.pb') print " >> Cleaning CO emission in Z-CMa." tclean(vis = target+'_TM2.split.cal.contsub', imagename = target+'_CO_cube', field = target, stokes = 'I', spw = '8', outframe = 'LSRK', restfreq = '230.538GHz', specmode = 'cube', width = 1, nchan = nChanCO, start = startCO, cell=[myCell], imsize=[myImsize,myImsize], #covers about the FoV deconvolver = 'multiscale', scales=[0,5,15],# inpixels equivalent to point source and 1, 3 times the beam. niter = 10000000, weighting = 'briggs', robust = 0.5, usemask = 'auto-multithresh', sidelobethreshold = 2.0, noisethreshold = 4.25, minbeamfrac = 0.3, lownoisethreshold = 1.5, negativethreshold = 0.0, gridder = 'standard', pbcor = True, threshold = myThreshCO, interactive = True, restoringbeam = '' ) """ CLEANING OF THE 13^ CS CUBE """ #---uncommeny below lines if you want to delete existing image files ---# #--- with same file name before imaging ---# #os.system('rm -rf '+target+'_13CS_cube.image') #os.system('rm -rf '+target+'_13CS_cube.image.pbcor') #os.system('rm -rf '+target+'_13CS_cube.psf') #os.system('rm -rf '+target+'_13CS_cube.mask') #os.system('rm -rf '+target+'_13CS_cube.model') #os.system('rm -rf '+target+'_13CS_cube.residual') #os.system('rm -rf '+target+'_13CS_cube.sumwt') #os.system('rm -rf '+target+'_13CS_cube.pb') """print " >> Now Cleaning 13CS emission in Z-CMa." tclean(vis = target+'_TM2.split.cal.contsub', imagename = target+'_13CS_cube', field = target, stokes = 'I', spw = '9', outframe = 'LSRK', restfreq = '231.2206GHz', specmode = 'cube', width = 1, nchan = nChanCS, start = startCS, cell=[myCell], imsize=[myImsize,myImsize], #covers about the FoV deconvolver = 'multiscale', scales=[0,5,15],# inpixels equivalent to point source and 1, 3 times the beam. niter = 10000000, weighting = 'briggs', robust = 0.5, usemask = 'auto-multithresh', sidelobethreshold = 2.0, noisethreshold = 4.25, minbeamfrac = 0.3, lownoisethreshold = 1.5, negativethreshold = 0.0, gridder = 'standard', pbcor = True, threshold = myThreshCS, interactive = True, restoringbeam = '' )""" #-------------------------------# #--- STEP 5: Moment Analysis ---# if step5: print '\n!!! STEP 5: Moment Analysis !!!\n' """ Finally we generate 0th, 1st and 2nd order image moments of the CO """ """ emission from Z-CMa. This gives us integrated line intensity, """ """ the velocity fields and the velocity distributions of the CO. """ print " >> Now Generating redshifted CO moments 0,1 & 2 for Z-CMa." immoments(imagename = 'Z_CMa_CO_cube.image', moments =[0,1,2], axis='spectral', chans =myChansRed, excludepix = [-100000.0,0.2], outfile = 'Z_CMa_CO_cube.image_red' ) print " >> Now Generating blueshifted CO moments 0,1 & 2 for Z-CMa." immoments(imagename = 'Z_CMa_CO_cube.image', moments =[0,1,2], axis='spectral', chans =myChansBlue, excludepix = [-100000.0,0.2], outfile = 'Z_CMa_CO_cube.image_blue' ) #-------------------------------# #--- STEP 6: Moment Blanking ---# #---- Code by Ana Karla :D -----# if step6: imagename = 'Z_CMa_CO_cube.image_red' threshold=0.07*5.0 mask = "'"+imagename+".integrated'>="+str(threshold) #masking using integrated moment print ">> Blanking moments using: "+imagename+".integrated" for end_name in ['.integrated','.weighted_coord','.weighted_dispersion_coord']: immath(imagename=imagename+end_name, expr='IM0',mask=mask, outfile=imagename+end_name+'_B') imagename = 'Z_CMa_CO_cube.image_blue' threshold=0.07*5.0 mask = "'"+imagename+".integrated'>="+str(threshold) #masking using integrated moment print ">> Blanking moments using: "+imagename+".integrated" for end_name in ['.integrated','.weighted_coord','.weighted_dispersion_coord']: immath(imagename=imagename+end_name, expr='IM0',mask=mask, outfile=imagename+end_name+'_B') #--- END THE STEPS ---# #---------------------#