|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################### 00002 # # 00003 # Regression/Benchmarking Script for NGC 4826 # 00004 # # 00005 ############################################### 00006 import time 00007 import os 00008 00009 #Data must be pre-filled 00010 #n4826_16apr98.ms, n4826_22apr98.ms 00011 00012 #N4826 - BIMA SONG Data 00013 #16apr98 00014 # source=ngc4826 00015 # phasecal=1310+323 00016 # fluxcal=3c273, Flux = 23 Jy on 16apr98 00017 # passcal= none - apparently, this is considered optional... odd 00018 #22apr98 00019 # source=ngc4826 00020 # phasecal=1230+123 00021 # fluxcal=3c273, Flux = 23 Jy on 22apr98 00022 # passcal= none 00023 # 00024 ## miriad: source Vlsr = 408; delta is 20 km/s 00025 # 00026 #NOTE: This data has been filled into MIRIAD, line-length correction 00027 # done, and then exported as separate files for each source. 00028 # 3c273 was not line length corrected since it was observed 00029 # for such a short amount of time that it did not need it. 00030 # 00031 ########################################################### 00032 ## Fill Data: 00033 ## 00034 00035 # miriad files copied from /home/cluster3/jmcmulli/ALMATST2/Data/N4826/ 00036 # 98apr16/ has 1310+323.ll/ 3c273/ ngc4826.ll/ 00037 # 98apr22/ has 1230+123.ll/ 3c273/ ngc4826.ll/ 00038 00039 ########################################################### 00040 ## Mosaic field spacing looks like: 00041 ## 00042 ## F3 (fieldid 4) F2 (fieldid 3) 00043 ## 00044 ## F4 (fieldid 5) F0 (fieldid 1) F1 (fieldid 2) 00045 ## 00046 ## F5 (fieldid 6) F6 (fieldid 7) 00047 ## 00048 ## 4x64 channels = 256 channels 00049 ## 00050 ## Primary Beam should be about 1.6' FWHM (7m dishes, 2.7mm wavelength) 00051 ## Resolution should be about 5-8" 00052 ############################################################# 00053 00054 pathname=os.environ.get('CASAPATH').split()[0] 00055 datapath=pathname+'/data/regression/ATST2/NGC4826/' 00056 00057 os.system('rm -rf n4826_22apr.* n4826_16apr.* n4826_tboth.ms') 00058 os.system('rm -rf srca.split.ms src.split.ms') 00059 os.system('rm -rf tmosaic* gcal*.split.* tgcal* n4826_t*') 00060 00061 startTime=time.time() 00062 startProc=time.clock() 00063 00064 print '--Copy/reinitialize--' 00065 #22 APR Calibration 00066 copystring22apr='cp -r '+datapath+'n4826_22apr98.ms n4826_22apr.ms' 00067 os.system(copystring22apr) 00068 tb.open('n4826_22apr.ms/HISTORY', nomodify=False) 00069 tb.removerows(range(tb.nrows())) 00070 tb.done() 00071 ###reset the data as pre calibrated state 00072 clearcal(vis='n4826_22apr.ms') 00073 ## Set the flux density of 3C273 to 23 Jy 00074 copytime=time.time() 00075 00076 ####set flux scale fo calibrater 00077 print '--setjy - 22apr98--' 00078 default('setjy') 00079 setjy(vis='n4826_22apr.ms',field='0',fluxdensity=[23.0,0.,0.,0.],scalebychan=False) 00080 setjy1time=time.time() 00081 ## Flag bad data non-interactively 00082 print '--flag data - 22apr97--' 00083 default('tflagdata') 00084 #flagdata(vis="n4826_22apr.ms", 00085 # antennaid=-1,baseline=[-1],chans=[-1], 00086 # clipfunction="ABS",clipcorr="YY", 00087 # clipminmax=[0.0, 80.0], 00088 # fieldid=-1,field="",spwid=-1,timerange="",unflag=False) 00089 tflagdata(vis="n4826_22apr.ms", mode='clip', 00090 correlation='ABS_YY', 00091 clipminmax=[0.0,80.0], 00092 clipoutside=True) 00093 #setclip=["ABS YY",[0.0, 80.0],True]) 00094 flag1time=time.time() 00095 ## Derive gain calibration solutions, try VLA-like calibration: 00096 print '--gaincal - 22par98--' 00097 default('gaincal') 00098 gaincal(vis='n4826_22apr.ms',caltable='n4826_22apr.gcal', 00099 field='0,1',spw='0', gaintype='G', 00100 minsnr=2.0, 00101 refant='ANT5',gaincurve=False,opacity=0.0,solint='inf',combine='') 00102 gaincal1time=time.time() 00103 ## Transfer the flux density scale: 00104 print '--fluxscale - 22apr98--' 00105 default('fluxscale') 00106 fluxscale(vis='n4826_22apr.ms',caltable='n4826_22apr.gcal', 00107 fluxtable='n4826_22apr.fcal', 00108 reference='3C273-F0',transfer=['1230+123-F0']); 00109 fluxscale1time=time.time() 00110 # 00111 ## Correct the calibrater/target source data: 00112 ## Use spwmap to apply gain solutions derived from spwid1 to all other spwids... 00113 ## 00114 print '--correct 22apr98--' 00115 default('applycal') 00116 applycal(vis='n4826_22apr.ms', 00117 field='1',spw='0', 00118 gaincurve=False,opacity=0.0,gaintable='n4826_22apr.fcal') 00119 applycal(vis='n4826_22apr.ms', 00120 field='2~8',spw='1~4', 00121 gaincurve=False,opacity=0.0,gaintable='n4826_22apr.fcal',spwmap=[0]) 00122 correct1time=time.time() 00123 00124 ## Split out calibrated target source and calibrater data: 00125 print '--split - 22apr98--' 00126 default('split') 00127 split(vis='n4826_22apr.ms',outputvis='gcal.split.ms', 00128 # field=1,spw=0,nchan=1,start=0,step=1,datacolumn='corrected') 00129 field='1',spw='0:0~1',datacolumn='corrected') 00130 default('split') 00131 split(vis='n4826_22apr.ms',outputvis='src.split.ms', 00132 # field=[2,3,4,5,6,7,8],spw=[1,2,3,4],nchan=64,start=0,step=1, 00133 field='2~8',spw='1~4:0~63', 00134 datacolumn='corrected') 00135 split1time=time.time() 00136 00137 ## Image the calibrater data: 00138 ## 00139 print '--image cal - 22apr98--' 00140 default('clean') 00141 clean(vis='gcal.split.ms',imagename='tgcal', 00142 cell=[1.,1.],imsize=[256,256],interpolation='nearest', 00143 field='0',spw='0',threshold=10.,ftmachine='ft', 00144 psfmode='clark',niter=100,stokes='I') 00145 imagecal1time=time.time() 00146 ## Image the target source mosaic: 00147 ## 00148 print '--image src - 22apr98--' 00149 default('clean') 00150 clean(vis='src.split.ms',imagename='tmosaicb', 00151 nchan=30,start=47,width=4,interpolation='nearest', 00152 # spw=range(0,3),field=range(7), 00153 spw='0,1,2',field='0,1,2,3,4,5,6', 00154 cell=[1.,1.],imsize=[256,256], 00155 stokes='I',mode='channel', 00156 psfmode='clark',imagermode='mosaic', ftmachine='mosaic', 00157 niter=300, 00158 scaletype='SAULT', 00159 pbcor=False, 00160 minpb=0.01, 00161 cyclefactor=2.0, usescratch=False, mask=[66,72,195,185]) 00162 imagesrc1time=time.time() 00163 ## Write image to a fits file: 00164 ## 00165 print '--write fits image--' 00166 default('exportfits') 00167 exportfits(imagename='tmosaicb.image',fitsimage='tmosaicb.fits') 00168 writefits1time=time.time() 00169 00170 00171 ################## 00172 #16 APR Calibration 00173 ########################################################### 00174 ## Concatenate the separate sources 00175 ## 00176 print '--copy/initialize - 16apr98 --' 00177 copystring16apr='cp -r '+datapath+'n4826_16apr98.ms n4826_16apr.ms' 00178 os.system(copystring16apr) 00179 tb.open('n4826_16apr.ms/HISTORY', nomodify=False) 00180 tb.removerows(range(tb.nrows())) 00181 tb.done() 00182 clearcal(vis='n4826_16apr.ms') 00183 copy2time=time.time() 00184 ## Set the flux density of 3C273 to 23 Jy 00185 print '--setjy - 16apr98--' 00186 default('setjy') 00187 setjy(vis='n4826_16apr.ms',field='0',fluxdensity=[23.0,0.,0.,0.],scalebychan=False) 00188 setjy2time=time.time() 00189 ## Flag end channels 00190 print '--flagdata - 16apr98 --' 00191 default('tflagdata') 00192 #flagdata(vis='n4826_16apr.ms',chans=[0,1,62,63],spwid=[2,3,4,5],fieldid=-1) 00193 tflagdata(vis='n4826_16apr.ms',spw='2~5:0;1;62;63', mode='manual') 00194 flagdata2time=time.time() 00195 ## Derive gain calibration solutions, try VLA-like calibration: 00196 print '--gaincal - 16apr98 --' 00197 default('gaincal') 00198 gaincal(vis='n4826_16apr.ms',caltable='n4826_16apr.gcal', 00199 field='0,1',spw='0,1', gaintype='G', 00200 minsnr=2.0, 00201 refant='ANT5',gaincurve=False,opacity=0.0,solint='inf',combine='') 00202 gaincal2time=time.time() 00203 # Found 14 good G Jones solutions. 00204 ## Transfer the flux density scale: 00205 print '--fluxscale - 16apr98 --' 00206 default('fluxscale') 00207 fluxscale(vis='n4826_16apr.ms',caltable='n4826_16apr.gcal', 00208 fluxtable='n4826_16apr.fcal', 00209 reference='3C273-F0',transfer=['1310+323-F0'],refspwmap=[0,0]); 00210 fluxscale2time=time.time() 00211 # Found reference field(s): 3C273-F0 00212 # Found transfer field(s): 1310+323-F0 00213 # Spw=2 will be referenced to spw=1 00214 # Flux density for 1310+323-F0 in SpW=2 (ref SpW=1) is: 00215 #1.47184 +/- 0.00787043 (SNR = 187.009) 00216 ## Correct the calibrater/target source data: 00217 ## Use new parm spwmap to apply gain solutions derived from spwid1 00218 ## to all other spwids... 00219 print '--correct - 16apr98 --' 00220 default('applycal') 00221 applycal(vis='n4826_16apr.ms', 00222 field='1',spw='1', 00223 gaincurve=False,opacity=0.0,gaintable='n4826_16apr.fcal') 00224 applycal(vis='n4826_16apr.ms', 00225 field='2~8',spw='2~5', 00226 gaincurve=False,opacity=0.0,gaintable='n4826_16apr.fcal',spwmap=[1]) 00227 correct2time=time.time() 00228 ## Split out calibrated target source and calibrater data: 00229 print '--split - 16apr98 --' 00230 default('split') 00231 split(vis='n4826_16apr.ms',outputvis='gcala.split.ms', 00232 # field=1,spw=1,nchan=1,start=0,step=1,datacolumn='corrected') 00233 field='1',spw='1:0~1',datacolumn='corrected') 00234 default('split') 00235 split(vis='n4826_16apr.ms',outputvis='srca.split.ms', 00236 # field=[2,3,4,5,6,7,8],spw=[2,3,4,5],nchan=64,start=0,step=1, 00237 field='2~8',spw='2~5:0~63', 00238 datacolumn='corrected') 00239 split2time=time.time() 00240 00241 # Extra flagging 00242 #flagdata(vis="srca.split.ms", 00243 # antennaid=[5],baseline=[-1],chans=[-1], 00244 # clipfunction="ABS",clipcorr="I",clipminmax=[0.0, 0.0], 00245 # fieldid=[-1],field="",spwid=-1, 00246 # timerange=['16-APR-1998/09:42:39.0', '16-APR-1998/10:24:46.0'], 00247 # unflag=False) 00248 default('tflagdata') 00249 tflagdata(vis="srca.split.ms", mode='clip', 00250 antenna='5', 00251 correlation="ABS_I", 00252 clipminmax=[0.0,0.0], 00253 clipoutside=True, 00254 timerange='1998/04/16/09:42:39.0~1998/04/16/10:24:46.0') 00255 #setclip=["ABS I",[0.0,0.0],True], 00256 00257 ## Image the calibrater data: 00258 print '--image cal - 16apr98 --' 00259 default('clean') 00260 clean(vis='gcala.split.ms',imagename='tgcala', 00261 cell=[1.,1.],imsize=[256,256],interpolation='nearest', 00262 field='0',spw='0',threshold=10.,ftmachine='ft', 00263 psfmode='clark',niter=100,stokes='I') 00264 imagecal2time=time.time() 00265 ############################################################# 00266 ## Image the target source mosaic: 00267 ## 00268 print '--image src - 16apr98 --' 00269 default('clean') 00270 ### mosaic data ...Sault weighting implies a noise unform image 00271 ### hence flux scale across image needed to be corrected to get right value 00272 clean(vis='srca.split.ms',imagename='tmosaica', 00273 nchan=30,start=47,width=4,interpolation='nearest', 00274 spw='0~2',field='0~6', 00275 cell=[1.,1.],imsize=[256,256], 00276 stokes='I',mode='channel', 00277 psfmode='clark',niter=300,imagermode='mosaic', ftmachine='mosaic',scaletype='SAULT',cyclefactor=3, usescratch=False) 00278 ################################################### 00279 ## Write image to a fits file: 00280 ## 00281 exportfits('tmosaica.image','tmosaica.fits') 00282 00283 #Combine and Image 00284 print '-- Initialize src.split.ms --' 00285 os.system('cp -r srca.split.ms n4826_tboth.ms'); 00286 cb.open('src.split.ms') 00287 cb.initcalset() 00288 cb.close() 00289 print '-- Concatenate --' 00290 ms.open(thems='n4826_tboth.ms',nomodify=False); 00291 ms.concatenate(msfile='src.split.ms',freqtol='50MHz') 00292 ms.close() 00293 00294 00295 #miriad:source velocity is 408; delta is 20 km/s; 24 maps 00296 default('clean') 00297 clean(vis='n4826_tboth.ms',imagename='tmosaic', 00298 nchan=30,start=46,width=4, 00299 spw='0~2',field='0~6',interpolation='nearest', 00300 cell=[1.,1.],imsize=[256,256], 00301 stokes='I',mode='channel', 00302 psfmode='clark',niter=500,imagermode='mosaic',ftmachine='mosaic',scaletype='SAULT', 00303 cyclefactor=3, usescratch=False) 00304 00305 ia.open(infile='tmosaic.image'); 00306 ia.moments(outfile='n4826_tmom0.im', 00307 moments=0,axis=3,includepix=[0.070,1000.0], 00308 mask='indexin(3,[3:24])') 00309 ia.moments(outfile='n4826_tmom1.im', 00310 moments=1,axis=3,includepix=[0.007,1000.0], 00311 mask='indexin(3,[3:24])') 00312 ia.close() 00313 00314 ################################################### 00315 ## Write image to a fits file: 00316 ## 00317 exportfits(imagename='tmosaic.image',fitsimage='tmosaic.fits') 00318 00319 endProc=time.clock() 00320 endTime=time.time() 00321 00322 # Regression 00323 # 22 APR 00324 ms.open('gcal.split.ms') 00325 thistest_cal_22apr=pl.mean(ms.range(["amplitude"]).get("amplitude")) 00326 ms.close() 00327 ms.open('src.split.ms') 00328 thistest_src_22apr=pl.mean(ms.range(["amplitude"]).get("amplitude")) 00329 ms.close 00330 ia.open('tgcal.image') 00331 statistics=ia.statistics() 00332 im_calmax22=statistics['max'][0] 00333 ia.close() 00334 ia.open('tmosaicb.image') 00335 statistics=ia.statistics() 00336 im_srcmax22=statistics['max'][0] 00337 00338 # 16 APR 00339 ms.open('gcala.split.ms') 00340 thistest_cal_16apr=pl.mean(ms.range(["amplitude"]).get("amplitude")) 00341 ms.close() 00342 ms.open('srca.split.ms') 00343 thistest_src_16apr=pl.mean(ms.range(["amplitude"]).get("amplitude")) 00344 ms.close 00345 ia.open('tgcala.image') 00346 statistics=ia.statistics() 00347 ia.close() 00348 im_calmax16=statistics['max'][0] 00349 ia.open('tmosaica.image') 00350 statistics=ia.statistics() 00351 ia.close() 00352 im_srcmax16=statistics['max'][0] 00353 00354 # COMBINED 00355 ms.open('n4826_tboth.ms') 00356 thistest_src=pl.mean(ms.range(["amplitude"]).get("amplitude")) 00357 ms.close() 00358 ia.open('n4826_tmom0.im') 00359 statistics=ia.statistics() 00360 ia.close() 00361 thistest_immax=statistics['max'][0] 00362 thistest_imrms=statistics['rms'][0] 00363 00364 ##changed values due to flat noise clean 00365 calmean22=10.51 00366 srcmean22=122.21 # 111.52 00367 calmax22=3.0459 00368 #srcmax22=1.54 00369 srcmax22=1.632 00370 calmean16=4.3269 00371 srcmean16=147.43 00372 srcmax16=1.584 00373 calmax16=1.43 00374 #srcmean=156.9898 00375 srcmean=147.43 00376 #immax=151.4 00377 immax=162.22 00378 #imrms=14.96 00379 imrms=13.18 00380 00381 # old (pre-minsnr=2) 00382 #calmean22=17.99 00383 #srcmean22=233.11 00384 #calmax22=2.9977 00385 #srcmax22=1.63 00386 #calmean16=75.292 00387 #srcmean16=155.521 00388 #srcmax16=1.45 00389 #calmax16=1.43 00390 #srcmean=233.5284 00391 #immax=157.6 00392 #imrms=16.18 00393 00394 diff_cal22apr=abs((calmean22-thistest_cal_22apr)/calmean22) 00395 diff_src22apr=abs((srcmean22-thistest_src_22apr)/srcmean22) 00396 diff_calmax22=abs((calmax22-im_calmax22)/calmax22) 00397 diff_srcmax22=abs((srcmax22-im_srcmax22)/srcmax22) 00398 diff_cal16apr=abs((calmean16-thistest_cal_16apr)/calmean16) 00399 diff_src16apr=abs((srcmean16-thistest_src_16apr)/srcmean16) 00400 diff_calmax16=abs((calmax16-im_calmax16)/calmax16) 00401 diff_srcmax16=abs((srcmax16-im_srcmax16)/srcmax16) 00402 00403 diff_src=abs((srcmean-thistest_src)/srcmean) 00404 diff_immax=abs((immax-thistest_immax)/immax) 00405 diff_imrms=abs((imrms-thistest_imrms)/imrms) 00406 00407 import datetime 00408 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00409 outfile='n4826.'+datestring+'.log' 00410 logfile=open(outfile,'w') 00411 00412 print >>logfile,'********** Regression ***********' 00413 print >>logfile,'* (Values: got, then expected) *' 00414 statstrs = {True: 'Passed', False: 'FAILED'} 00415 print >>logfile, '*', statstrs[diff_cal22apr<0.08], 'cal mean amp (22apr)' 00416 print >>logfile,'--Cal mean amp (22apr) '+str(thistest_cal_22apr)+', '+str(calmean22) 00417 print >>logfile, '*', statstrs[diff_src22apr<0.08], 'src mean amp (22apr)' 00418 print >>logfile,'--Src mean amp (22apr) '+str(thistest_src_22apr)+', '+str(srcmean22) 00419 print >>logfile, '*', statstrs[diff_calmax22<0.08], 'cal image max (22apr)' 00420 print >>logfile,'--Image max (cal;22apr) '+str(im_calmax22)+', '+str(calmax22) 00421 print >>logfile, '*', statstrs[diff_srcmax22<0.08], 'src image max (22apr)' 00422 print >>logfile,'--Image max (src;22apr): '+str(im_srcmax22)+', '+str(srcmax22) 00423 print >>logfile, '*', statstrs[diff_cal16apr<0.08], 'cal mean amp (16apr)' 00424 print >>logfile,'--Cal mean amp (16apr) '+str(thistest_cal_16apr)+', '+str(calmean16) 00425 print >>logfile, '*', statstrs[diff_src16apr<0.08], 'src mean amp (16apr)' 00426 print >>logfile,'--Src mean amp (16apr) '+str(thistest_src_16apr)+', '+str(srcmean16) 00427 print >>logfile, '*', statstrs[diff_calmax16<0.08], 'cal image max (16apr)' 00428 print >>logfile,'--Image max (cal;16apr): '+str(im_calmax16)+', '+str(calmax16) 00429 print >>logfile, '*', statstrs[diff_srcmax16<0.08], 'src image max (16apr)' 00430 print >>logfile,'--Image max (src;16apr): '+str(im_srcmax16)+', '+str(srcmax16) 00431 print >>logfile, '*', statstrs[diff_src<0.08], 'src mean amplitude test' 00432 print >>logfile,'--Src mean amp '+str(thistest_src)+', '+str(srcmean) 00433 print >>logfile, '*', statstrs[diff_immax<0.08], 'image max test' 00434 print >>logfile,'--Image max '+str(thistest_immax)+', '+str(immax) 00435 print >>logfile, '*', statstrs[diff_imrms<0.08], 'image rms test' 00436 print >>logfile,'--Image rms '+str(thistest_imrms)+', '+str(imrms) 00437 00438 if ((diff_cal22apr<0.08) & (diff_src22apr<0.08) & (diff_cal16apr<0.08) & (diff_src16apr<0.08) &(diff_src<0.08) & (diff_immax<0.08) & (diff_imrms<0.08)): 00439 regstate=True 00440 print >>logfile,'---' 00441 print >>logfile,'Passed Regression test for NGC 4826 Mosaic' 00442 print >>logfile,'---' 00443 else: 00444 regstate=False 00445 print >>logfile,'----FAILED Regression test for NGC 4826 Mosaic' 00446 print >>logfile,'*********************************' 00447 print >>logfile,'' 00448 print >>logfile,'********* Benchmarking *****************' 00449 print >>logfile,'* *' 00450 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 00451 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 00452 print >>logfile,'Processing rate MB/s was: '+str(300./(endTime - startTime)) 00453 print >>logfile,'* Breakdown: *' 00454 00455 logfile.close()
1.8.0
