# ALMA Data Reduction Script # Calibration thesteps = [] step_title = {0: 'Stokes dirty images of Polarization calibrator', 1: 'Stokes clean images of Polarization calibrator', 2: 'Stokes images of Phase calibrator (optional)', 3: 'IQU Stokes clean images of the Target', 4: 'IQU Stokes (LINE) clean images of the Target', 5: 'Total intensity imaging of the check source'} try: print 'List of steps to be executed ...', mysteps thesteps = mysteps except: print 'global variable mysteps not set.' if (thesteps==[]): thesteps = range(0,len(step_title)) print 'Executing all steps: ', thesteps # The Python variable 'mysteps' will control which steps # are executed when you start the script using # execfile('scriptForCalibration.py') # e.g. setting # mysteps = [2,3,4]# before starting the script will make the script execute # only steps 2, 3, and 4 # Setting mysteps = [] will make it execute all steps. import casadef import re if re.search('^5.1.1', casadef.casa_version) == None: sys.exit('ERROR: PLEASE USE THE SAME VERSION OF CASA THAT YOU USED FOR GENERATING THE SCRIPT: 5.1.1') # To prevent a stray CR character in exported FITS #os.unsetenv('CASAVERSION') msname = 'calibrated.ms.cal' total='calibrated.ms' #msname2='calibrated.ms' #bandpass='1' # J0510+1800 phasecal = '1' # J0541-0211 polcalib = '0' # J0522-3627 check1='4' #J0532-0307 target = 'HH_212' # HH_212 refant = 'DV07' #listobs(vis=msname,listfile=msname+'.listobs') # CONTINUUM and line channel selection, using continuum subtraction flagmanager(vis=msname, mode='save', versionname='beforeflags') flagdata(vis=msname, field='HH_212', spw='3:60~100;180~220;280~340;440~470;500~540;680~720;840~900', flagbackup=False) split(vis=msname, width=[1,1,1,10], outputvis=msname+'.forcont', field='', datacolumn='data') flagmanager(vis=msname, mode='restore', versionname='beforeflags') msname2=msname+'.forcont' #uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2.mfs.IQUV.manual ### Stokes dirty images of the polarization calibrator mystep = 0 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] ''' Os.system('rm -rf Pol_Cal.Stokes.dirty*') clean(vis=msname, field=polcalib, imagename='Pol_Cal.Stokes.dirty', cell=['0.015arcsec'], imsize=[2048,2048], spw='0,1,2', stokes='IQUV', psfmode='clarkstokes', interactive=False,niter=0) exportfits(imagename='Pol_Cal.Stokes.dirty.image', fitsimage='Pol_Cal.Stokes.dirty.image.fits') # export the corrected image s ''' ### Stokes clean images of the polarization calibrator mystep = 1 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2_3.mfs.IQUV.manual*') clean(vis=msname2, field=polcalib, imagename='uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2_3.mfs.IQUV.manual', cell=['0.012arcsec'], imsize=[2048,2048], stokes='IQUV', psfmode='clarkstokes', spw='', robust=0.5, threshold='0.015mJy', interactive=False,niter=500) exportfits(imagename='uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2_3.mfs.IQUV.manual.image',fitsimage='uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2_3.mfs.IQUV.manual.image.fits') os.system('rm -rf Pol_Cal.StokesI.clean.image*') immath(imagename='uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2_3.mfs.IQUV.manual.image',outfile='Pol_Cal.StokesI.clean.image',expr='IM0',stokes='I') os.system('rm -rf Pol_Cal.StokesQ.clean.image*') immath(imagename='uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2_3.mfs.IQUV.manual.image',outfile='Pol_Cal.StokesQ.clean.image',expr='IM0',stokes='Q') os.system('rm -rf Pol_Cal.StokesU.clean.image*') immath(imagename='uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2_3.mfs.IQUV.manual.image',outfile='Pol_Cal.StokesU.clean.image',expr='IM0',stokes='U') os.system('rm -rf Pol_Cal.StokesV.clean.image*') immath(imagename='uid___A001_X1273_X1f9.J0522-3627_polleak.spw0_1_2_3.mfs.IQUV.manual.image',outfile='Pol_Cal.StokesV.clean.image',expr='IM0',stokes='V') resI=imfit(imagename = 'Pol_Cal.StokesI.clean.image', box = '1010,1010,1030,1030') resQ=imfit(imagename = 'Pol_Cal.StokesQ.clean.image', box = '1010,1010,1030,1030') resU=imfit(imagename = 'Pol_Cal.StokesU.clean.image', box = '1010,1010,1030,1030') # and then we extract the flux and error values for each Stokes fluxI=resI['results']['component0']['flux']['value'][0] errorI=resI['results']['component0']['flux']['error'][0] fluxQ=resQ['results']['component0']['flux']['value'][1] errorQ=resQ['results']['component0']['flux']['error'][1] fluxU=resU['results']['component0']['flux']['value'][2] errorU=resU['results']['component0']['flux']['error'][2] #Now we use these values to compute polarization intensity, angle and ratio, and their errors: fluxPI = sqrt( fluxQ**2 + fluxU**2 ) errorPI = sqrt( (fluxQ*errorU)**2 + (fluxU*errorQ)**2 ) / fluxPI fluxPImjy = 1000*fluxPI errPImjy = 1000*errorPI polRatio = fluxPI / fluxI errPRat = polRatio * sqrt( (errorPI/fluxPI)**2 + (errorI/fluxI)**2 ) polAngle = 0.5 * degrees( atan2(fluxU,fluxQ) ) errPA = 0.5 * degrees( errorPI / fluxPI ) print 'Pol ratio of Polarization calibrator: ', polRatio print 'Pol angle of Polarization calibrator: ', polAngle #Pol ratio of Polarization calibrator: 0.0423581616031 #Pol angle of Polarization calibrator: -74.5517197748 # Polarization fraction reported by last ALMA measurement is ~ . # uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual ### Stokes images of the phase calibrator mystep = 2 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual*') clean(vis=msname2, field=phasecal, imagename='uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual', cell=['0.012arcsec'], imsize=[2048,2048], stokes='IQUV', psfmode='clarkstokes', spw='', robust=0.5, threshold='0.05mJy', interactive=False,niter=100) exportfits(imagename='uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual.image',fitsimage='uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual.image.fits') os.system('rm -rf Phase_Cal.StokesI.clean.image*') immath(imagename='uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual.image',outfile='Phase_Cal.StokesI.clean.image',expr='IM0',stokes='I') os.system('rm -rf Phase_Cal.StokesQ.clean.image*') immath(imagename='uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual.image',outfile='Phase_Cal.StokesQ.clean.image',expr='IM0',stokes='Q') os.system('rm -rf Phase_Cal.StokesU.clean.image*') immath(imagename='uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual.image',outfile='Phase_Cal.StokesU.clean.image',expr='IM0',stokes='U') os.system('rm -rf Phase_Cal.StokesV.clean.image*') immath(imagename='uid___A001_X1273_X1f9.J0541-0211_ph.spw0_1_2_3.mfs.IQUV.manual.image',outfile='Phase_Cal.StokesV.clean.image',expr='IM0',stokes='V') resI=imfit(imagename = 'Phase_Cal.StokesI.clean.image', box = '1010,1010,1030,1030') resQ=imfit(imagename = 'Phase_Cal.StokesQ.clean.image', box = '1010,1010,1030,1030') resU=imfit(imagename = 'Phase_Cal.StokesU.clean.image', box = '1010,1010,1030,1030') # and then we extract the flux and error values for each Stokes fluxI=resI['results']['component0']['flux']['value'][0] errorI=resI['results']['component0']['flux']['error'][0] fluxQ=resQ['results']['component0']['flux']['value'][1] errorQ=resQ['results']['component0']['flux']['error'][1] fluxU=resU['results']['component0']['flux']['value'][2] errorU=resU['results']['component0']['flux']['error'][2] #Now we use these values to compute polarization intensity, angle and ratio, and their errors: fluxPI = sqrt( fluxQ**2 + fluxU**2 ) errorPI = sqrt( (fluxQ*errorU)**2 + (fluxU*errorQ)**2 )/fluxPI fluxPImjy = 1000*fluxPI errPImjy = 1000*errorPI polRatio = fluxPI / fluxI errPRat = polRatio * sqrt((errorPI/fluxPI)**2+(errorI/fluxI)**2) polAngle = 0.5 * degrees(atan2(fluxU,fluxQ)) errPA = 0.5 * degrees(errorPI/fluxPI) print 'Pol ratio of Phase calibrator: ', polRatio print 'Pol angle of Phase calibrator: ', polAngle #Pol ratio of Phase calibrator: 0.0291542513402 #Pol angle of Phase calibrator: -75.588994553 # uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual ### IQU Stokes clean images of the Target mystep = 3 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] # continuum using spw 4 os.system('rm -rf uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual*') clean(vis=msname2, field=target, imagename='uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual', cell=['0.012arcsec'], imsize=[2048,2048], stokes='IQU', psfmode='clarkstokes', outframe='LSRK', interactive=False, niter=1000, weighting='briggs', spw='', threshold='0.05mJy', # uvtaper=False, # outertaper=['0.09arcsec'], robust=0.5) os.system('rm -rf HH212.StokesI.clean.image*') immath(imagename='uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual.image',outfile='HH212.StokesI.clean.image',expr='IM0',stokes='I') os.system('rm -rf HH212.StokesQ.clean.image*') immath(imagename='uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual.image',outfile='HH212.StokesQ.clean.image',expr='IM0',stokes='Q') os.system('rm -rf HH212.StokesU.clean.image*') immath(imagename='uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual.image',outfile='HH212.StokesU.clean.image',expr='IM0',stokes='U') # Primary beam corrected Stokes image and fits convertion myimagebase = 'uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual' impbcor(imagename=myimagebase+'.image', pbimage=myimagebase+'.flux', outfile=myimagebase+'.image.pbcor', overwrite=True) # perform PBcorr exportfits(imagename=myimagebase+'.image.pbcor', fitsimage=myimagebase+'.image.pbcor.fits',overwrite=True) # export the corrected image exportfits(imagename=myimagebase+'.flux', fitsimage=myimagebase+'.flux.fits',overwrite=True) # export the PB image stokess=['I','Q','U'] for i in stokess: myimagebase = 'uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual_stokes_%s'%str(i) os.system('rm -rf uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual_stokes_%s.image'%str(i)) immath(imagename='uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual.image',outfile='uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual_stokes_%s.image'%str(i),expr='IM0',stokes=str(i)) # exportfits(imagename=myimagebase+'.image', fitsimage=myimagebase+'.image.fits',overwrite=True) # export the corrected image ### Polarization images production os.system('rm -rf HH212.POLI') immath(outfile='HH212.POLI', mode='poli', imagename=['uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual_stokes_Q.image','uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual_stokes_U.image'], sigma='0.0Jy/beam') os.system('rm -rf HH212.POLA') immath(outfile='HH212.POLA', mode='pola', imagename=['uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual_stokes_Q.image','uid___A001_X1273_X1f9.HH_212_sci.spw0_1_2_3.mfs.IQU.manual_stokes_U.image'], polithresh='0.15mJy/beam') # uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual ### (LINE) IQU Stokes clean images of the Target mystep = 4 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual*') clean(vis=msname, field=target, imagename='uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual', cell=['0.012arcsec'], mode='velocity', restfreq='346.559GHz', width='3.4km/s', nchan=100, start='-170.0km/s', imsize=[512,512], stokes='IQU', psfmode='clarkstokes', outframe='LSRK', interactive=False, niter=2000, weighting='briggs', spw='3', threshold='2mJy', robust=0.5) # uvtaper=T, # mask='HH212.spw3.clean.mask', # outertaper=['0.09arcsec']) # exporting images to fits myimagebase = 'uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual' impbcor(imagename=myimagebase+'.image', pbimage=myimagebase+'.flux', outfile=myimagebase+'.image.pbcor', overwrite=True) # perform PBcorr exportfits(imagename=myimagebase+'.image.pbcor', fitsimage=myimagebase+'.image.pbcor.fits',overwrite=True) # export the corrected image exportfits(imagename=myimagebase+'.flux', fitsimage=myimagebase+'.flux.fits',overwrite=True) # export the PB image stokess=['I','Q','U'] for i in stokess: # myimagebase = 'uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual' os.system('rm -rf uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual_stokes_%s.image'%str(i)) immath(imagename='uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual.image',outfile='uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual_stokes_%s.image'%str(i),expr='IM0',stokes=str(i)) # exportfits(imagename=myimagebase+'.image', fitsimage=myimagebase+'.image.fits',overwrite=True) # export the corrected image os.system('rm -rf HH212.spw3.clean.POLI') immath(outfile='HH212.spw3.clean.POLI', mode='poli', imagename=['uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual_stokes_Q.image','uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual_stokes_U.image'], sigma='0.0Jy/beam') os.system('rm -rf HH212.spw3.clean.POLA') immath(outfile='HH212.spw3.clean.POLA', mode='pola', imagename=['uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual_stokes_Q.image','uid___A001_X1273_X1f9.HH_212_sci.spw3.cube.IQU.manual_stokes_U.image'], polithresh='0.15mJy/beam') # uid___A001_X1273_X1f9.J0532-0307_chk.spw0_1_2_3.mfs.I.manual ### Total intensity imaging of the check source (Continuum only) mystep = 5 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print 'Step ', mystep, step_title[mystep] os.system('rm -rf uid___A001_X1273_X1f9.J0532-0307_chk.spw0_1_2_3.mfs.I.manual*'), tclean(vis = total, imagename = 'uid___A001_X1273_X1f9.J0532-0307_chk.spw0_1_2_3.mfs.I.manual', field = check1, # specmode = 'mfs', # restfreq='102.9GHz', # width='96.44km/s', # start='-200.0km/s', # nchan=100, outframe='LSRK', restoringbeam='common', spw='', interactive = False, niter = 100, imsize = [1024,1024], cell = '0.01arcsec', phasecenter = 'J2000 05h32m07.519317s -03d07m07.03705s', weighting = 'briggs', threshold='0.6mJy', robust = 0.5, datacolumn='data', deconvolver='hogbom', nterms=2, pbcor=False, chanchunks=-1, gridder='standard') # uvtaper=['0.5arcsec','0.55arcsec']) # impbcor(imagename='uid___A001_X1273_X1f9.J0532-0307_chk.spw0_1_2_3.mfs.I.manual.image',pbimage='uid___A001_X1273_X1f9.J0532-0307_chk.spw0_1_2_3.mfs.I.manual.pb',outfile='uid___A001_X1273_X1f9.J0532-0307_chk.spw0_1_2_3.mfs.I.manual.image.pbcor') exportfits(imagename='uid___A001_X1273_X1f9.J0532-0307_chk.spw0_1_2_3.mfs.I.manual.image',fitsimage='uid___A001_X1273_X1f9.J0532-0307_chk.spw0_1_2_3.mfs.I.manual.image.fits')