DG Notes on 20170904 Event

Purpose

There are various scripts for doing the imaging and self calibration, but they require explanation on how to use them. This page is meant to provide step-by-step explanation on analyzing one event, the M5.5 flare that occurred on 2017-Sep-04. I earlier did an initial study of this event, but in a hurried fashion in order to use the results for a proposal (which was successful, by the way!). I would like to do a more careful job with the initial calibration and analysis, and then continue with a second round of self calibration.

Initial Preparation

The data begin as raw data in Miriad format, recorded to disk by the pipeline system, in the form of IDB files. These files each contain 10-min of data, with several such files together representing a scan. The easiest way to identify the time period of interest is to examine the overview spectral plots on the RHESSI Browser and check "EOVSA Radio Data" on the upper left. Use the time you identified to find the corresponding IDB file(s) under /data1/eovsa/fits/IDB/yyyymmdd/ (also accessible from the EOVSA web site. The scan containing this flare is a quite long one, starting 18:53:33 UT and continuing to 21:30:00 UT or so.

To know the exact time of interest takes some iteration, especially for such a long flare with several peaks and smaller events. After some checking, I settled on the three 10-min files starting at 20:23:33, 20:33:33, and 20:43:33 UT. Undoubtedly an additional one or two 10-min intervals after the end time of 20:53:33 would be useful for following the decay of the event, but in the interests of time analysis of these will be postponed.

Importing EOVSA Data

Once the desired timerange is known, the next step is to convert the raw IDB files to calibrated CASA ms files. The script below does this from within CASA, and must be run on pipeline, which has access to the raw data and the SQL database. This script may require 30 minutes or more to run.

from suncasa.tasks import task_calibeovsa as calibeovsa
from suncasa.tasks import task_importeovsa as timporteovsa
import dump_tsys as dt
from util import Time
import numpy as np
import os
# import copy
get_tp = True
workpath = '/data1/dgary/solar/20170904_Mflare/'
if not os.path.exists(workpath):
    os.makedirs(workpath)
os.chdir(workpath)
# This time range will read files IDB20170904202333, IDB20170904203333, IDB20170904204333, and IDB20170904205333
trange = Time(['2017-09-04 20:20:00', '2017-09-04 21:00:00'])
info = dt.rd_fdb(Time('2017-09-04'))
sidx = np.where(
    np.logical_and(info['SOURCEID'] == 'Sun', info['PROJECTID'] == 'NormalObserving') & np.logical_and(
        info['ST_TS'].astype(np.float) >= trange[0].lv,
        info['ST_TS'].astype(np.float) <= trange[1].lv))
filelist = info['FILE'][sidx]
outpath = 'msdata/'
if not os.path.exists(outpath):
    os.makedirs(outpath)
inpath = '/data1/eovsa/fits/IDB/{}/'.format(trange[0].datetime.strftime("%Y%m%d"))
ncpu = 1
#
#
timporteovsa.importeovsa(idbfiles=[inpath + ll for ll in filelist], ncpu=ncpu, timebin="0s", width=1, visprefix=outpath,
                         nocreatms=False,
                         doconcat=False, modelms="", doscaling=False, keep_nsclms=False, udb_corr=True)
msfiles = [outpath + ll + '.ms' for ll in filelist]
# invis=copy.copy(msfiles)
concatvis = os.path.basename(msfiles[0])[:11] + '_2023-2103.ms'
vis = calibeovsa.calibeovsa(msfiles, caltype=['refpha', 'phacal'], interp='nearest', doflag=True, flagant='13~15',
                            doimage=False, doconcat=True,
                            msoutdir='msdata', concatvis=concatvis, keep_orig_ms=True)

This will result in the four measurement sets, IDB20170904202333.ms, IDB20170904203333.ms, IDB20170904204333.ms, and IDB20170904205333.ms being created and written to the current directory, and then concatenated into a single ms with the name IDB20170904_2023-2103.ms.

At this stage, the calibrated, concatenated ms can be transferred to baozi for further processing. For example, the command

tar -czf IDB20170904_2023-2103.ms.tgz IDB20170904_2023-2103.ms

creates a compressed tar file that can be copied directly to baozi if the appropriate tunnel exists, or it can be staged in two steps. The file is 1 GB in size, so it could take some time to transfer. Actually, it only takes 5 minutes when I open a tunnel to pipeline on local port 8666, and then use

scp -P 8666 user@localhost:/data1/dgary/solar/20170904_Mflare/msdata/IDB20170904_2023-2103.ms.tgz .

in MobaXterm.

Selfcal on Baozi

Once the file is transferred to Baozi (I untarred the file into /srg/dgary/eovsa/solar/20170904_Mflare/msdata), the first thing that is needed is to set up the environment:

source /srg/.setenv_baozi

The strategy for the first stage of selfcal is to make a clean image for a short duration of the flare, when the source structure is relatively simple and the signal-to-noise ratio is high, then use that model to correct the antenna-based phases to bring the resulting map closer to the model. Three iterations of phase-only selfcal are done in each of two stages.

The initial setup involves declaring the working directory and ms names, and splitting out a single polarization (this one is 'XX'): from suncasa.utils import helioimage2fits as hf import os import numpy as np

1. task handles

dofullsun=0 doslfcal=1 dosbd=0 dofinalclean=0 workdir='/srg/dgary/eovsa/solar/20170904_Mflare/msdata/' refms = workdir+'IDB20170904_2023-2103.ms' slfcalms = workdir+'IDB20170904_2023-2103.ms.XX.peak' slfcaledms = workdir+'slfcal/IDB20170904_2023_2103.ms.XX.peak.slfcaled' if not os.path.exists(slfcalms):

   split(vis=refms,outputvis=slfcalms,correlation='XX')
   listobs(slfcalms,listfile=slfcalms+'.listobs')

clearcal(slfcalms) delmod(slfcalms) spws=[str(s) for s in range(1,31)] antennas='0~12' pol='XX' calprefix=workdir+'slfcal/caltbs/slf_2023' imgprefix=workdir+'slfcal/images/slf_2023'

To make the clean image in CASA, one must first declare a phase center in RA, Dec coordinates, and a clean box (mask) within which to apply the clean algorithm. One way to do this is to generate a full Sun image for a relatively wide timerange, bandwidth, and largish cell size during the event, and then measure the coordinates from the image. This is done with the following clean command: if dofullsun:

   trange='2017/09/04/20:33:33~2017/09/04/20:43:33'
   img_init=imgprefix+'_init'
   os.system('rm -rf '+img_init+'*')
   clean(vis=slfcalms,
           antenna='0~10,12',
           imagename=img_init,
           spw='1~15',
           mode='mfs',
           timerange=trange,
           imagermode='csclean',
           psfmode='clark',
           imsize=[512,512],
           cell=['5arcsec'],
           niter=50,
           gain=0.05,
           stokes=pol,
           restoringbeam=['15arcsec'],
           interactive=False,
           usescratch=True)

   viewer(img_init+'.image')

The last line displays the image in the CASA image viewer, from which the coordinates can be read with the mouse. From this image, the coordinate center of the emission was measured to be 'J2000 10h53m57 06d56m47'. This will be used as the phase center of the subsequent maps used as the selfcal model.

First Stage

There is a trade-off in doing selfcal in various bands, so the first stage calculates the model in overlapping CASA band ranges 1~5, 4~10, 8~15, 10~20, and 15~30. The EOVSA band numbers are 4 greater than these. Although the model image is made using grouped bands, that model is applied separately to individual bands within each range, i.e. the 1~5 model is applied to bands 1~4, the 4~10 model is applied to 5~9, the 8~15 model is applied to 10~14, the 15~20 model is applied to 15~19, and the 15~30 model is applied to 20~30. A script was used as below, where the setting of trange (trange='20:45:00~20:45:10'), phase center (phasecenter='J2000 10h53m57 06d56m47') and the clean mask (mask='circle [[128pix, 128pix], 60pix]') are the main changes needed from case to case.

if doslfcal:

   os.system('rm -rf '+calprefix+'*')
   os.system('rm -rf '+imgprefix+'*')
   trange='20:45:00~20:45:10'
   tb.open(slfcalms+'/SPECTRAL_WINDOW')
   reffreqs=tb.getcol('REF_FREQUENCY')
   bdwds=tb.getcol('TOTAL_BANDWIDTH')
   cfreqs=reffreqs+bdwds/2.
   tb.close()
   sbeam=40.
   strtmp=[t.replace(':',) for t in trange.split('~')]
   timestr='t'+strtmp[0]+'-'+strtmp[1]
   refantenna='0'
   nround=3
   niters=[100,100,200]
   robusts=[1.0,1.0,1.0]
   doapplycal=[1,1,1]
   uvranges=[,,]
   slftbs=[]
   for n in range(nround):
       slfcal_tb_g= calprefix+'.G'+str(n)
       for sp in spws:
           cfreq=cfreqs[int(sp)]
           bm=max(sbeam*cfreqs[1]/cfreq,10.)
           slfcal_img = imgprefix+'.spw'+sp.zfill(2)+'.slfcal'+str(n)
           if int(sp) < 5:
               spwran='1~5'
           if int(sp) >= 5 and int(sp) < 10:
               spwran='4~10'
           if int(sp) >= 10 and int(sp) < 15:
               spwran='8~15'
           if int(sp) >= 15 and int(sp) < 20:
               spwran='10~20'
           if int(sp) >= 20:
               spwran='15~30'
           print 'using spw {0:s} as model'.format(spwran)
           try:
               clean(vis=slfcalms,
                       antenna=antennas,
                       imagename=slfcal_img,
                       uvrange=uvranges[n],
                       #spw=sp,
                       spw=spwran,
                       mode='mfs',
                       timerange=trange,
                       imagermode='csclean',
                       psfmode='clark',
                       imsize=[256,256],
                       cell=['3arcsec'],
                       niter=niters[n],
                       gain=0.05,
                       stokes=pol,
                       weighting='briggs',
                       robust=robusts[n],
                       phasecenter='J2000 10h53m57 06d56m47',
                       #mask='box [ [ 75pix , 90pix] , [205pix, 165pix ] ]',
                       mask='circle [ [ 128pix ,128pix] , 60pix ]',
                       restoringbeam=[str(bm)+'arcsec'],
                       interactive=False,
                       usescratch=True)

           except:
               print 'error in cleaning spw: '+sp
               print 'using nearby spws for initial model'
               sp_e=int(sp)+2
               sp_i=int(sp)-2
               if sp_i < 0:
                   sp_i = 0
               if sp_e > 30:
                   sp_e = 30
               sp_=str(sp_i)+'~'+str(sp_e)
               try:
                   tget(clean)
                   spw=sp_
                   clean()
               except:
                   print 'still not successful. abort...'
                   break

           print 'processing spw: '+sp
           #gain solution, phase only
           gaincal(vis=slfcalms, refant=refantenna,antenna=antennas,caltable=slfcal_tb_g,spw=sp, uvrange=,\
                   #gaintable=slftbs,selectdata=True,timerange=trange,solint='inf',gaintype='G',calmode='p',\
                   gaintable=[],selectdata=True,timerange=trange,solint='inf',gaintype='G',calmode='p',\
                   combine=,minblperant=2,minsnr=3,append=True)
           if not os.path.exists(slfcal_tb_g):
               #listcal(vis=slfcalms, caltable=slfcal_table)
               #print 'solutions found in spw: '+sp
               print 'No solution found in spw: '+sp

       if os.path.exists(slfcal_tb_g):
           #slftbs.append(slfcal_tb_g)
           slftbs=[slfcal_tb_g]
           plotcal(caltable=slfcal_tb_g,antenna=antennas,xaxis='antenna',yaxis='phase',\
                   subplot=421,plotrange=[0,12,-180,180],iteration='spw')

       if doapplycal[n]:
           clearcal(slfcalms)
           delmod(slfcalms)
           applycal(vis=slfcalms,gaintable=slftbs,spw=','.join(spws),selectdata=True,\
                    antenna=antennas,interp='nearest',flagbackup=False,applymode='calonly',calwt=False)
       if n < nround-1:
           prompt=raw_input('Continuing to selfcal?')
           if prompt.lower() == 'n':
               if os.path.exists(slfcaledms):
                   os.system('rm -rf '+slfcaledms)
               split(slfcalms,slfcaledms,datacolumn='corrected')
               print 'Final calibrated ms is {0:s}'.format(slfcaledms)
               break
           if prompt.lower() == 'y':
               slfcalms_=slfcalms+str(n)
               if os.path.exists(slfcalms_):
                   os.system('rm -rf '+slfcalms_)
               split(slfcalms,slfcalms_,datacolumn='corrected')
               slfcalms=slfcalms_
       else:
           if os.path.exists(slfcaledms):
               os.system('rm -rf '+slfcaledms)
           split(slfcalms,slfcaledms,datacolumn='corrected')
           print 'Final calibrated ms is {0:s}'.format(slfcaledms)