|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 00002 ############################################### 00003 # # 00004 # Regression/Benchmarking Script for H121 # 00005 # 5 MAY # 00006 ############################################### 00007 import time 00008 import os 00009 # Data must be pre-filled 00010 # l.ms, l2.ms 00011 00012 pathname=os.environ.get('CASAPATH').split()[0] 00013 datapath=pathname+'/data/regression/ATST1/H121/*3mmmay*.ms' 00014 00015 startTime=time.time() 00016 startProc=time.clock() 00017 00018 print '--Copy MS--' 00019 # Fill - really copy over pristine MS 00020 os.system('rm -rf *3mmmay*.ms h121*.c* h121*.3* 0224* 2200* nrao150*') 00021 copystring='cp -r '+datapath+' .' 00022 os.system(copystring) 00023 clearcal(vis='05-3mmmay-1998-h121.ms') 00024 clearcal(vis='09-3mmmay-1998-h121.ms') 00025 msfile1='05-3mmmay-1998-h121.ms' 00026 msfile2='09-3mmmay-1998-h121.ms' 00027 00028 copytime=time.time() 00029 00030 ## Set flux density scale for the flux density calibrater MWC349 00031 default('setjy') 00032 setjy(vis=msfile1,spw='7',field='1',scalebychan=False,fluxdensity=[1.12,0.,0.,0.]) 00033 00034 ## Get preliminary time-dependent phase solutions to improve 00035 ## coherent average for bandpass solution and derive bandpass 00036 ## calibration for 3mm LSB 00037 default('gaincal') 00038 gaincal(vis=msfile1,caltable='h121.3mm.ph.gcal0', 00039 field='4',spw='7:8~55', 00040 gaintype='GSPLINE',calmode='p',splinetime=10000., 00041 refant='1',phasewrap=260,gaincurve=False,opacity=0.0,preavg=120.) 00042 00043 print '--bandpass (3mm)--' 00044 #derive bandpass calibration 00045 # solint value 00046 default('bandpass') 00047 bandpass(vis=msfile1,caltable='h121.3mm.bpoly', 00048 field='4',spw='7:8~55', 00049 bandtype='BPOLY',degamp=6,degphase=12,solnorm=False, 00050 maskcenter=6,maskedge=0,refant='1', 00051 gaintable='h121.3mm.ph.gcal0',gaincurve=False,opacity=0.0) 00052 00053 # 3mm CO Band 00054 default('gaincal') 00055 gaincal(vis=msfile1,caltable='2200.3mm.ph.gcal0', 00056 field='0',spw='3:64~191', 00057 gaintype='GSPLINE',calmode='p',splinetime=10000., 00058 refant='1',preavg=10.,phasewrap=260,gaincurve=False,opacity=0.0) 00059 default('bandpass') 00060 bandpass(vis=msfile1,caltable='2200.3mm.bpoly', 00061 field='0',spw='3:8~247', 00062 bandtype='BPOLY',degamp=6,degphase=12,solnorm=False, 00063 maskcenter=6,maskedge=0,refant='1', 00064 gaintable='2200.3mm.ph.gcal0',gaincurve=False,opacity=0.0) 00065 00066 # Determine new and better phase solutions for 3mm LSB for all calibrators 00067 default('gaincal') 00068 gaincal(vis=msfile1,caltable='h121.3mm.ph.gcal', 00069 field='0,1,2,4',spw='7:8~55', 00070 gaintype='GSPLINE',calmode='p',splinetime=5000.,refant='1', 00071 phasewrap=260, 00072 gaincurve=False,opacity=0.0,gaintable='h121.3mm.bpoly',preavg=0.) 00073 00074 #Apply all solutions derived so far, determine calibrator's flux densities by solving for T 00075 #and use fluxscale 00076 default('gaincal') 00077 gaincal(vis=msfile1,caltable='h121.3mm.temp', 00078 field='0,1,2,4',spw='7:8~55', 00079 solint='600s',refant='1',gaintype='T', 00080 opacity=0.0,gaincurve=False, 00081 gaintable=['h121.3mm.ph.gcal','h121.3mm.bpoly']) 00082 00083 #fluxscale 00084 default('fluxscale') 00085 fluxscale(vis=msfile1,caltable='h121.3mm.temp',fluxtable='h121.3mm.flux', 00086 reference='MWC349*',transfer='NRAO150*,2200+420*,0224+671*') 00087 calphase3mmtime=time.time() 00088 #2200+420: 2.66 00089 #0224+671: 1.14 00090 #NRAO 150: 3.12 00091 00092 #use old values 00093 print '--Set fluxscale (setjy)--' 00094 default('setjy') 00095 setjy(vis=msfile1,field='0',spw='7',scalebychan=False,fluxdensity=[2.607,0.,0.,0.]) 00096 setjy(vis=msfile1,field='2',spw='7',scalebychan=False,fluxdensity=[1.107,0.,0.,0.]) 00097 setjy(vis=msfile1,field='4',spw='7',scalebychan=False,fluxdensity=[3.034,0.,0.,0.]) 00098 setjy3mmtime=time.time() 00099 00100 ## Amplitude calibration of 3mm LSB: 00101 ## 00102 ## phase solutions will be pre-applied as well as carried forward 00103 ## to the output solution table. 00104 print '--gaincal amp (3mm)--' 00105 default('gaincal') 00106 gaincal(vis=msfile1,caltable='h121.3mm.amp.gcal', 00107 field='0,1,2,4',spw='7:8~55', 00108 gaintype='GSPLINE',calmode='a',splinetime=20000.,refant='1', 00109 phasewrap=260,gaincurve=False,opacity=0.0, 00110 preavg=2500., 00111 gaintable=['h121.3mm.ph.gcal','h121.3mm.bpoly']) 00112 calamp3mmtime=time.time() 00113 00114 ## Correct the target source and all other 3mm LSB data 00115 ## 00116 ## note that only the 60 central channels will be calibrated 00117 ## since the BPOLY solution is only defined for these 00118 default('applycal') 00119 applycal(vis=msfile1, 00120 spw='7', 00121 gaincurve=False,opacity=0.0, 00122 gaintable=['h121.3mm.ph.gcal','h121.3mm.amp.gcal','h121.3mm.bpoly']) 00123 # Correct Target CO (1-0) 00124 default('applycal') 00125 applycal(vis=msfile1, 00126 spw='3', 00127 gaincurve=False,opacity=0.0, 00128 gaintable=['h121.3mm.ph.gcal','h121.3mm.amp.gcal','2200.3mm.bpoly']) 00129 00130 # Split calibrated target source data 00131 print '--split calibrater--' 00132 default('split') 00133 split(vis=msfile1,outputvis='h121.3mm.split.ms', 00134 # field=3,spw=7,nchan=48,start=8,step=1,datacolumn='corrected') 00135 field='3',spw='7:8~55',datacolumn='corrected') 00136 default('split') 00137 split(vis=msfile1,outputvis='nrao150.3mm.split.ms', 00138 # field=4,spw=7,nchan=48,start=8,step=1,datacolumn='corrected') 00139 field='4',spw='7:8~55',datacolumn='corrected') 00140 default('split') 00141 split(vis=msfile1,outputvis='h121.co10.split.ms', 00142 # field=3,spw=3,nchan=240,start=8,step=1,datacolumn='corrected') 00143 field='3',spw='3:8~247',datacolumn='corrected') 00144 splitsrctime=time.time() 00145 00146 ## Get a first image the target source in 3 mm continuum emission: 00147 default('clean') 00148 clean(vis='nrao150.3mm.split.ms',imagename='nrao150.3mm', 00149 psfmode='clark',imagermode='',niter=500,gain=0.1,nchan=1,start=0,width=48, 00150 spw='0',field='0',stokes='I',interpolation='nearest', 00151 weighting='briggs',robust=2.,cell=[1.,1.], 00152 imsize=[64,64],mode='mfs') 00153 default('clean') 00154 clean(vis='h121.co10.split.ms',imagename='h121.co10', 00155 psfmode='clark',imagermode='',niter=1000,gain=0.1,threshold=20., 00156 nchan=80,start=1,width=3,spw='0',interpolation='nearest', 00157 field='0',stokes='I', 00158 weighting='briggs',robust=2., 00159 cell=[0.5,0.5],imsize=[128,128],mode='channel') 00160 image1time=time.time() 00161 00162 ############################################### 00163 # # 00164 # Regression/Benchmarking Script for H121 # 00165 # 9 MAY # 00166 ############################################### 00167 00168 ## Set flux density scale for the flux density calibrater MWC349 00169 setjy(vis=msfile2,spw='7',field='4',scalebychan=False,fluxdensity=[1.12,0.,0.,0.]) 00170 00171 ## Get preliminary time-dependent phase solutions to improve 00172 ## coherent average for bandpass solution and derive bandpass 00173 ## calibration for 3mm LSB 00174 default('gaincal') 00175 gaincal(vis=msfile2,caltable='h121b.3mm.ph.gcal0', 00176 field='1',spw='7:8~55', 00177 gaintype='GSPLINE',calmode='p',splinetime=10000., 00178 refant='1',phasewrap=260, 00179 gaincurve=False,opacity=0.0,preavg=120.) 00180 00181 print '--bandpass (3mm)--' 00182 #derive bandpass calibration 00183 # solint value 00184 default('bandpass') 00185 bandpass(vis=msfile2,caltable='h121b.3mm.bpoly', 00186 field='1',spw='7:8~55', 00187 bandtype='BPOLY',degamp=6,degphase=12,solnorm=False, 00188 maskcenter=6,maskedge=0,refant='1', 00189 gaintable='h121.3mm.ph.gcal0',gaincurve=False,opacity=0.0) 00190 00191 # 3mm CO Band 00192 default('gaincal') 00193 gaincal(vis=msfile2,caltable='2200b.3mm.ph.gcal0', 00194 field='0',spw='3:64~191', 00195 gaintype='GSPLINE',calmode='p',splinetime=10000., 00196 refant='1',preavg=10.,phasewrap=260, 00197 gaincurve=False,opacity=0.0) 00198 default('bandpass') 00199 bandpass(vis=msfile2,caltable='2200b.3mm.bpoly', 00200 field='0',spw='3:8~247', 00201 bandtype='BPOLY',degamp=6,degphase=12,solnorm=False, 00202 maskcenter=6,maskedge=0,refant='1', 00203 gaintable='2200b.3mm.ph.gcal0',gaincurve=False,opacity=0.0) 00204 00205 # Determine new and better phase solutions for 3mm LSB for all calibrators 00206 default('gaincal') 00207 gaincal(vis=msfile2,caltable='h121b.3mm.ph.gcal', 00208 field='0,1,3,4',spw='7:8~55', 00209 gaintype='GSPLINE',calmode='p',splinetime=5000.,refant='1', 00210 phasewrap=260, 00211 gaincurve=False,opacity=0.0,gaintable='h121b.3mm.bpoly',preavg=0.) 00212 00213 #Apply all solutions derived so far, determine calibrator's flux densities by solving for T 00214 #and use fluxscale 00215 default('gaincal') 00216 gaincal(vis=msfile2,caltable='h121b.3mm.temp', 00217 field='0,1,3,4',spw='7:8~55', 00218 solint='600s',refant='1',gaintype='T', 00219 opacity=0.0,gaincurve=False, 00220 gaintable=['h121b.3mm.ph.gcal','h121b.3mm.bpoly']) 00221 00222 #fluxscale 00223 default('fluxscale') 00224 fluxscale(vis=msfile2,caltable='h121b.3mm.temp',fluxtable='h121b.3mm.flux', 00225 reference='MWC349*',transfer='3C454*,2200+420*,0224+671*') 00226 calphase3mmtime=time.time() 00227 #2200+420: 3.29 00228 #0224+671: 1.53 00229 #C454.3 : 6.127 00230 00231 #use old values 00232 print '--Set fluxscale (setjy)--' 00233 default('setjy') 00234 setjy(vis=msfile2,field='0',spw='7',scalebychan=False,fluxdensity=[3.235,0.,0.,0.]) 00235 setjy(vis=msfile2,field='1',spw='7',scalebychan=False,fluxdensity=[1.507,0.,0.,0.]) 00236 setjy(vis=msfile2,field='3',spw='7',scalebychan=False,fluxdensity=[5.982,0.,0.,0.]) 00237 setjy3mmtime=time.time() 00238 00239 ## Amplitude calibration of 3mm LSB: 00240 ## 00241 ## phase solutions will be pre-applied as well as carried forward 00242 ## to the output solution table. 00243 print '--gaincal amp (3mm)--' 00244 default('gaincal') 00245 gaincal(vis=msfile2,caltable='h121b.3mm.amp.gcal', 00246 field='0,1,3,4',spw='7:8~55', 00247 gaintype='GSPLINE',calmode='a',splinetime=20000.,refant='1', 00248 phasewrap=260, 00249 gaincurve=False,opacity=0.0, 00250 preavg=2500., 00251 gaintable=['h121b.3mm.ph.gcal','h121b.3mm.bpoly']) 00252 calamp3mmtime=time.time() 00253 00254 ## Correct the target source and all other 3mm LSB data 00255 ## 00256 ## note that only the 60 central channels will be calibrated 00257 ## since the BPOLY solution is only defined for these 00258 default('applycal') 00259 applycal(vis=msfile2, 00260 spw='7', 00261 gaincurve=False,opacity=0.0, 00262 gaintable=['h121b.3mm.ph.gcal','h121b.3mm.amp.gcal','h121b.3mm.bpoly']) 00263 # Correct Target CO (1-0) 00264 default('applycal') 00265 applycal(vis=msfile2, 00266 spw='3', 00267 gaincurve=False,opacity=0.0, 00268 gaintable=['h121b.3mm.ph.gcal','h121b.3mm.amp.gcal','2200b.3mm.bpoly']) 00269 00270 # Split calibrated target source data 00271 print '--split calibrater--' 00272 default('split') 00273 split(vis=msfile2,outputvis='h121b.3mm.split.ms', 00274 # field=2,spw=7,nchan=48,start=8,step=1,datacolumn='corrected') 00275 field='2',spw='7:8~55',datacolumn='corrected') 00276 default('split') 00277 split(vis=msfile2,outputvis='0224b.3mm.split.ms', 00278 # field=1,spw=7,nchan=48,start=8,step=1,datacolumn='corrected') 00279 field='1',spw='7:8~55',datacolumn='corrected') 00280 default('split') 00281 split(vis=msfile2,outputvis='h121b.co10.split.ms', 00282 # field=2,spw=3,nchan=240,start=8,step=1,datacolumn='corrected') 00283 field='2',spw='3:8~247',datacolumn='corrected') 00284 splitsrctime=time.time() 00285 00286 ## Get a first image the target source in 3 mm continuum emission: 00287 default('clean') 00288 clean(vis='0224b.3mm.split.ms',imagename='0224b.3mm', 00289 psfmode='clark',imagermode='',niter=500,gain=0.1,nchan=1,start=0,width=48, 00290 spw='0',field='0',stokes='I',interpolation='nearest', 00291 weighting='briggs',robust=2., 00292 cell=[1.,1.],imsize=[64,64],mode='mfs') 00293 default('clean') 00294 clean(vis='h121b.co10.split.ms',imagename='h121b.co10', 00295 psfmode='clark',imagermode='',niter=1000,gain=0.1,threshold=20., 00296 nchan=80,start=1,width=3,interpolation='nearest', 00297 spw='0',field='0',stokes='I', 00298 weighting='briggs',robust=2., 00299 cell=[1.,1.],imsize=[64,64],mode='channel') 00300 image1time=time.time() 00301 00302 #Concatenate datasets 00303 #print '--concatenate datasets--' 00304 #os.system('cp -r h121.3mm.split.ms h121.3mm.concat.ms') 00305 #os.system('cp -r h121.co10.split.ms h121.co10.concat.ms') 00306 #default('concat') 00307 #concat(concatvis='h121.3mm.concat.ms',vis='h121b.3mm.split.ms', 00308 # freqtol='0.1MHz',dirtol='0.1arcsec', timesort=True) 00309 #default('concat') 00310 #concat(concatvis='h121.co10.concat.ms',vis='h121b.co10.split.ms', 00311 # freqtol='0.1MHz',dirtol='0.1arcsec', timesort=True) 00312 #concattime=time.time() 00313 00314 # skip concat 00315 # use multims input in clean (to test CAS-1625) 00316 00317 #Image target in 3mm continuum emission 00318 default('clean') 00319 #clean(vis='h121.3mm.concat.ms',imagename='h121all.3mm', 00320 # psfmode='clark',niter=2000,gain=0.1,nchan=1,start=0,width=48, 00321 # spw='0',field='0',stokes='I',interpolation='nearest', 00322 # weighting='briggs',robust=2., 00323 # cell=[1.,1.],imsize=[64,64],mode='mfs') 00324 clean(vis=['h121.3mm.split.ms','h121b.3mm.split.ms'],imagename='h121all.3mm', 00325 psfmode='clark',imagermode='',niter=2000,gain=0.1,nchan=1,start=0,width=48, 00326 spw='0',field='0',stokes='I',interpolation='nearest', 00327 weighting='briggs',robust=2., 00328 cell=[1.,1.],imsize=[64,64],mode='mfs') 00329 default('clean') 00330 #clean(vis='h121.co10.concat.ms',imagename='h121c.co10', 00331 # psfmode='clark',niter=100,gain=0.1,nchan=80,start=1,width=3, 00332 # spw='0',field='0',stokes='I',interpolation='nearest', 00333 # weighting='briggs',robust=.5, 00334 # cell=[0.2,0.2],imsize=[256,256],mode='channel') 00335 clean(vis=['h121.co10.split.ms','h121b.co10.split.ms'],imagename='h121c.co10', 00336 psfmode='clark',imagermode='',niter=100,gain=0.1,nchan=80,start=1,width=3, 00337 spw='0',field='0',stokes='I',interpolation='nearest', 00338 weighting='briggs',robust=.5, 00339 cell=[0.2,0.2],imsize=[256,256],mode='channel') 00340 00341 endProc=time.clock() 00342 endTime=time.time() 00343 00344 # Regression 00345 00346 ms.open('h121.3mm.split.ms') 00347 thistest_3mm=max(ms.range(['amplitude']).get('amplitude')) 00348 ms.close() 00349 ms.open('h121b.3mm.split.ms') 00350 thistest_3mmb=max(ms.range(['amplitude']).get('amplitude')) 00351 ms.close() 00352 ia.open('h121.co10.image') 00353 statistics=ia.statistics() 00354 co10a=statistics['max'][0] 00355 co10arms=statistics['rms'][0] 00356 ia.close() 00357 ia.open('h121b.co10.image') 00358 statistics=ia.statistics() 00359 co10b=statistics['max'][0] 00360 co10brms=statistics['rms'][0] 00361 ia.close() 00362 ia.open('h121c.co10.image') 00363 statistics=ia.statistics() 00364 co10c=statistics['max'][0] 00365 co10crms=statistics['rms'][0] 00366 ia.close() 00367 ia.open('h121all.3mm.image') 00368 statistics=ia.statistics() 00369 cont3mmmax=statistics['max'][0] 00370 cont3mmrms=statistics['rms'][0] 00371 ia.close() 00372 00373 srca=3.235 00374 srcb=3.260 00375 #coa=0.4755 00376 #cob=0.4043 00377 #coc=0.3300 00378 # new values for CO images after cvel-clean param 00379 # unification (last relevant change r.14424) 00380 # The values are from the run on 64b RHEL(rishi), 00381 # with active # r.14436. Lower max values but 00382 # slightly better rms # as compared to ones get with casa3.1 00383 #coa=0.4172 00384 cob=0.3712 00385 coc=0.3362 00386 00387 # 2011-09-01: co10a is back up. 00388 coa = 0.4819 00389 00390 contmax=0.00963 00391 00392 diff_3mma=abs((srca-thistest_3mm)/srca) 00393 diff_3mmb=abs((srcb-thistest_3mmb)/srcb) 00394 diff_co10a=abs((coa-co10a)/coa) 00395 diff_co10b=abs((cob-co10b)/cob) 00396 diff_co10c=abs((coc-co10c)/coc) 00397 diff_cont=abs((contmax-cont3mmmax)/contmax) 00398 00399 import datetime 00400 datestring=datetime.datetime.isoformat(datetime.datetime.today()) 00401 outfile='h121.'+datestring+'.log' 00402 logfile=open(outfile,'w') 00403 00404 passfail = {True: '* Passed', False: '* FAILED'} 00405 print >>logfile,'********** Regression ***********' 00406 print >>logfile,'* *' 00407 print >>logfile, passfail[diff_3mma < 0.05], '3mm (5may) max amplitude test ' 00408 print >>logfile,'--3mm max (5may) amp '+str(thistest_3mm) 00409 print >>logfile, passfail[diff_3mmb < 0.05], '3mm (9may) max amplitude test ' 00410 print >>logfile,'--3mm max (9may) amp '+str(thistest_3mmb) 00411 print >>logfile, passfail[diff_co10a < 0.05], 'CO 1-0 image max (5may) test ' 00412 print >>logfile,'--CO 1-0 (5may) image max '+str(co10a) 00413 print >>logfile, passfail[diff_co10b < 0.05], 'CO 1-0 image max (9may) test ' 00414 print >>logfile,'--CO 1-0 (9may) image max '+str(co10b) 00415 print >>logfile, passfail[diff_co10c < 0.05], 'CO 1-0 image max (combined) test ' 00416 print >>logfile,'--CO 1-0 (combine) image max '+str(co10c) 00417 print >>logfile, passfail[diff_cont < 0.05], '3mm image max test ' 00418 print >>logfile,'--3mm image max '+str(cont3mmmax) 00419 00420 00421 if ((diff_3mma<0.05) & (diff_3mmb<0.05) & (diff_co10a<0.05) & (diff_co10b<0.05) & (diff_co10c<0.05) & (diff_cont<0.05)): 00422 regstate=True 00423 print >>logfile,'---' 00424 print >>logfile,'Passed Regression test for H121' 00425 print >>logfile,'---' 00426 print '' 00427 print 'Regression PASSED' 00428 print '' 00429 else: 00430 regstate=False 00431 print '' 00432 print 'Regression FAILED' 00433 print '' 00434 print >>logfile,'----FAILED Regression test for H121' 00435 print >>logfile,'*********************************' 00436 00437 print >>logfile,'' 00438 print >>logfile,'' 00439 print >>logfile,'********* Benchmarking *****************' 00440 print >>logfile,'* *' 00441 print >>logfile,'Total wall clock time was: '+str(endTime - startTime) 00442 print >>logfile,'Total CPU time was: '+str(endProc - startProc) 00443 print >>logfile,'Processing rate MB/s was: '+str(278./(endTime - startTime)) 00444 #print >>logfile,'* Breakdown: *' 00445 #print >>logfile,'* copy time was: '+str(copytime-startTime) 00446 #print >>logfile,'* flag time was: '+str(flagtime-copytime) 00447 #print >>logfile,'* cal ph 3mm time was: '+str(calphase3mmtime-flagtime) 00448 #print >>logfile,'* setjy 3mm time was: '+str(setjy3mmtime-calphase3mmtime) 00449 #print >>logfile,'* cal amp 3mm time was: '+str(calamp3mmtime-setjy3mmtime) 00450 #print >>logfile,'* correct 3mm time was: '+str(correct3mmtime-calamp3mmtime) 00451 #print >>logfile,'* split cal time was: '+str(splitcaltime-correct3mmtime) 00452 #print >>logfile,'* split src time was: '+str(splitsrctime-splitcaltime) 00453 #print >>logfile,'* image src time was: '+str(image3mmtime-splitsrctime) 00454 #print >>logfile,'* cal ph 1mm time was: '+str(calphase1mmtime-image3mmtime) 00455 #print >>logfile,'* setjy 1mm time was: '+str(setjy1mmtime-calphase1mmtime) 00456 #print >>logfile,'* cal amp 1mm time was: '+str(calamp1mmtime-setjy1mmtime) 00457 #print >>logfile,'* correct 1mm time was: '+str(correct1mmtime-calamp1mmtime) 00458 #print >>logfile,'* splitsrc 1mm time was: '+str(splitsrc1mmtime-correct1mmtime) 00459 #print >>logfile,'* image 1mm time was: '+str(image1mmtime-splitsrc1mmtime) 00460 #print >>logfile,'* HCO cal time was: '+str(hcocaltime-image1mmtime) 00461 #print >>logfile,'* image HCO time was: '+str(imagehcotime-hcocaltime) 00462 #print >>logfile,'* CO cal time was: '+str(cocaltime-imagehcotime) 00463 #print >>logfile,'* image CO time was: '+str(imagecotime-cocaltime) 00464 00465 logfile.close()
1.8.0
