# ALMA Data Reduction Script # $Id: scriptForImaging_FullPolTemplate.py,v 1.5 2022/10/20 15:40:35 dpetry Exp $ # Imaging thesteps = [] step_title = {0: 'Stokes cleaned images of Polarization calibrator', 1: 'Stokes images of Phase calibrator (optional)', 2: 'IQUV continuum image of the target', 3: 'IQUV cube image of the target at rep freq', 4: 'IQUV cube image of 13CO', 5: 'Polarization images production', 6: 'Export images to FITS format'} 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. thevis = ['concat_S1.ms.FDM.cal', 'concat_S2.ms.FDM.cal'] phasecal = '1' phasecalname = 'J0510+1800' polcal = '0' polcalname = 'J1256-0547' target = 'HH30' originalSPWIDs = '17_19_21_23_25' # Stokes images for target 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 '+polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.IQUV.manual*') tclean(vis = thevis, imagename = polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.IQUV.manual', field = polcal, stokes = 'IQUV', #spw = '0,1,2,3,4,5', spw = '', outframe = 'LSRK', specmode = 'cont', nterms = 1, imsize = [240, 240], # if you change this, you also need to change the subsequent imfit parameters cell = '0.05arcsec', deconvolver = 'clarkstokes', niter = 1000, weighting = 'briggs', robust = 0.5, mask = '', gridder = 'standard', pbcor = True, threshold = '', interactive = True ) myimagebase = polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.IQUV.manual' os.system('rm -rf '+polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.I.image*') immath(imagename=myimagebase+'.image',outfile=polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.I.image',expr='IM0',stokes='I') os.system('rm -rf '+polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.Q.image*') immath(imagename=myimagebase+'.image',outfile=polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.Q.image',expr='IM0',stokes='Q') os.system('rm -rf '+polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.U.image*') immath(imagename=myimagebase+'.image',outfile=polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.U.image',expr='IM0',stokes='U') os.system('rm -rf '+polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.V.image*') immath(imagename=myimagebase+'.image',outfile=polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.V.image',expr='IM0',stokes='V') resI=imfit(imagename = polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.I.image', box = '110,110,145,145') resQ=imfit(imagename = polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.Q.image', box = '115,115,130,130') resU=imfit(imagename = polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.U.image', box = '110,110,145,145') # 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: #Pol angle of Polarization calibrator: calstat=imstat(imagename=myimagebase+'.residual', axes=[0,1]) rms=(calstat['rms']) prms = (rms[1]**2. + rms[2]**2.)**0.5 os.system('rm -rf '+polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.P.manual.image') immath(outfile=polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.P.manual.image', mode='poli', imagename = polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.IQUV.manual.image', sigma='0.0Jy/beam') os.system('rm -rf '+polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.A.manual.image') immath(outfile=polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.A.manual.image', mode='pola', imagename = polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.IQUV.manual.image', polithresh='%.8fJy/beam'%(5.0*prms)) # Stokes images of Phase calibrator (optional) 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 '+phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.IQUV.manual') tclean(vis = thevis, imagename = phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.IQUV.manual', field = phasecal, stokes = 'IQUV', spw = '', outframe = 'LSRK', specmode = 'cont', nterms = 1, imsize = [500, 500], cell = '0.05arcsec', deconvolver = 'clarkstokes', niter = 1000, weighting = 'briggs', robust = 0.5, mask = '', gridder = 'standard', pbcor = True, threshold = '1mJy', interactive = True) os.system('rm -rf '+phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.I.manual*') immath(imagename=phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.IQUV.manual.image',outfile=phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.I.manual.image',expr='IM0',stokes='I') os.system('rm -rf '+phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.Q.manual.image*') immath(imagename=phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.IQUV.manual.image',outfile=phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.Q.manual.image',expr='IM0',stokes='Q') os.system('rm -rf '+phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.U.manual.image*') immath(imagename=phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.IQUV.manual.image',outfile=phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.U.manual.image',expr='IM0',stokes='U') os.system('rm -rf '+phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.V.manual.image*') immath(imagename=phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.IQUV.manual.image',outfile=phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.V.manual.image',expr='IM0',stokes='V') resI=imfit(imagename = phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.I.manual.image', box = '220,220,280,280') resQ=imfit(imagename = phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.Q.manual.image', box = '220,220,280,280') resU=imfit(imagename = phasecalname+'_ph.spw'+originalSPWIDs+'.mfs.U.manual.image', box = '220,220,280,280') # 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) # IQUV continuum image of the target 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 '+target+'_sci.spw19.mfs.IQUV.manual*') # use original SPW IDs tclean(vis = thevis, imagename = target+'_sci.spw19.mfs.IQUV.manual', # use original SPW IDs field = target, stokes = 'IQUV', spw = '1,5,10,15', # depending on concatenation outcome, this may be more than one SPW outframe = 'LSRK', specmode = 'mfs', nterms = 1, imsize = [700,700], cell = '0.05arcsec', deconvolver = 'clarkstokes', niter = 10000, weighting = 'briggs', robust = 0.5, mask = '', gridder = 'standard', pbcor = True, threshold = '0.0mJy', interactive = True ) # IQUV image of target mystep = 3 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) os.system('rm -rf '+target+'_sci.spw21.cube.IQUV.manual*') # use original SPW IDs tclean(vis = thevis, imagename = target+'_sci.spw21.cube.IQUV.manual', # use original SPW IDs field = target, stokes = 'IQUV', spw = '2,7,12,17', # depending on concatenation outcome, this may be more than one SPW outframe = 'LSRK', specmode = 'cube', restfreq='230.54GHz', start='-30km/s', width='8km/s', nchan=8, nterms = 1, imsize = [700,700], cell = '0.05arcsec', deconvolver = 'clarkstokes', niter = 10000, weighting = 'briggs', robust = 0.5, mask = '', gridder = 'standard', pbcor = True, threshold = '0.0mJy', # '0.4mJy', last run interactive interactive = True ) # IQUV image of 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 '+target+'_sci.spw25.cube.IQUV.manual*') # use original SPW IDs tclean(vis = thevis, imagename = target+'_sci.spw25.cube.IQUV.manual', # use original SPW IDs field = target, stokes = 'IQUV', spw = '4,9,14,19', # depending on concatenation outcome, this may be more than one SPW outframe = 'LSRK', specmode = 'cube', restfreq='220.398GHz', start='-30km/s', width='8km/s', nchan=8, nterms = 1, imsize = [700,700], cell = '0.05arcsec', deconvolver = 'clarkstokes', niter = 10000, weighting = 'briggs', robust = 0.5, mask = '', gridder = 'standard', pbcor = True, threshold = '0.4mJy', interactive = True ) # Polarization images production mystep = 5 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) myimagebase=target+'_sci.spw19.mfs' calstat=imstat(imagename= myimagebase+'.IQUV.manual.residual', axes=[0,1]) rms=(calstat['rms']) prms = (rms[1]**2. + rms[2]**2.)**0.5 os.system('rm -rf '+myimagebase+'.P.manual.image.pbcor') immath(outfile=myimagebase+'.P.manual.image.pbcor', mode='poli', imagename = myimagebase+'.IQUV.manual.image.pbcor', sigma='0.0Jy/beam') os.system('rm -rf '+myimagebase+'.A.manual.image.pbcor') immath(outfile=myimagebase+'.A.manual.image.pbcor', mode='pola', imagename = myimagebase+'.IQUV.manual.image.pbcor', polithresh='%.8fJy/beam'%(5.0*prms)) # Export images to FITS format mystep = 6 if(mystep in thesteps): casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO') print('Step ', mystep, step_title[mystep]) myimages = [target+'_sci.spw19.mfs.IQUV.manual', target+'_sci.spw19.mfs.P.manual', target+'_sci.spw19.mfs.A.manual', target+'_sci.spw21.cube.IQUV.manual', target+'_sci.spw25.cube.IQUV.manual', polcalname+'_polleak.spw'+originalSPWIDs+'.mfs.IQUV.manual'] for myimagebase in myimages: exportfits(imagename = myimagebase+'.image.pbcor', fitsimage = myimagebase+'.pbcor.fits', overwrite = True) if os.path.exists(myimagebase+'.pb'): exportfits(imagename = myimagebase+'.pb', fitsimage = myimagebase+'.pb.fits', overwrite = True)