|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################################################# 00002 # $Id:$ 00003 # Test Name: # 00004 # ngc3256-analysis-regression.py # 00005 # An ALMA Science Verification Data Analysis Regression # 00006 # using observations of NGC3256 from April 2011 # 00007 # # 00008 # Rationale for Inclusion: # 00009 # Need tests of ALMA analysis chain with selfcal # 00010 # # 00011 # Input data: # 00012 # six MSs # 00013 # # 00014 ############################################################################# 00015 00016 00017 step_title = { 0 : 'Data preparation', 00018 1 : 'Generate caltables', 00019 2 : 'Apriori flagging', 00020 3 : 'Delay cal for antenna DV07', 00021 4 : 'Apply wvr and delay cal tables, split out corrected data', 00022 5 : 'Apply tsys correction', 00023 6 : 'Concatenate', 00024 7 : 'More flagging and fixplanets on Titan', 00025 8 : 'Fast phase-only gaincal on the bandpass calibrator', 00026 9 : 'Bandpass calibration', 00027 10 : 'Set fluxscale on Titan', 00028 11 : 'Slow amplitude and phase gaincal', 00029 12 : 'Adjust absolute flux scale', 00030 13 : 'Apply bandpass and gain calibration', 00031 14 : 'Final flagging', 00032 15 : 'Final calibration', 00033 16 : 'Image the phase calibrator', 00034 17 : 'Apply calibration to flux calibrator', 00035 18 : 'Image flux calibrator', 00036 19 : 'Split out the calibrated target data', 00037 20 : 'Image the NGC3256 continuum', 00038 21 : 'Slow phase selfcal', 00039 22 : 'Image the NGC3256 continuum again after selfcal', 00040 23 : 'Determine and subtract continuum', 00041 24 : 'Clean the NGC3256 CO line cube', 00042 25 : 'Evaluate the NGC3256 CO line cube', 00043 26 : 'Clean the NGC3256 CNhi line cube', 00044 27 : 'Evaluate the NGC3256 CNhi line cube', 00045 28 : 'Clean the NGC3256 CNlo line cube', 00046 29 : 'Evaluate the NGC3256 CNlo line cube', 00047 30 : 'Verify regression results' 00048 } 00049 00050 # global defs 00051 basename=['uid___A002_X1d54a1_X5','uid___A002_X1d54a1_X174','uid___A002_X1d54a1_X2e3', 00052 'uid___A002_X1d5a20_X5','uid___A002_X1d5a20_X174','uid___A002_X1d5a20_X330'] 00053 makeplots=True 00054 therefant = 'DV04' 00055 00056 # regression results storage 00057 rmspcal = 0 00058 peakpcal = 0 00059 rmstitan = 0 00060 peaktitan = 0 00061 rmscont = 0 00062 peakcont = 0 00063 rmscontsc = 0 00064 peakcontsc = 0 00065 00066 ############################# 00067 00068 # Some infrastructure to make repeating individual parts 00069 # of this workflow more convenient. 00070 00071 thesteps = [] 00072 00073 try: 00074 print 'List of steps to be executed ...', mysteps 00075 thesteps = mysteps 00076 except: 00077 print 'global variable mysteps not set.' 00078 if (thesteps==[]): 00079 thesteps = range(0,len(step_title)) 00080 print 'mysteps empty. Executing all steps: ', thesteps 00081 00082 # The Python variable 'mysteps' will control which steps 00083 # are executed when you start the script using 00084 # execfile('alma-m100-analysis-regression.py') 00085 # e.g. setting 00086 # mysteps = [2,3,4] 00087 # before starting the script will make the script execute 00088 # only steps 2, 3, and 4 00089 # Setting mysteps = [] will make it execute all steps. 00090 00091 totaltime = 0 00092 inittime = time.time() 00093 ttime = inittime 00094 steptime = [] 00095 00096 def timing(): 00097 global totaltime 00098 global inittime 00099 global ttime 00100 global steptime 00101 global step_title 00102 global mystep 00103 global thesteps 00104 thetime = time.time() 00105 dtime = thetime - ttime 00106 steptime.append(dtime) 00107 totaltime += dtime 00108 ttime = thetime 00109 casalog.origin('TIMING') 00110 casalog.post( 'Step '+str(mystep)+': '+step_title[mystep], 'WARN') 00111 casalog.post( 'Time now: '+str(ttime), 'WARN') 00112 casalog.post( 'Time used this step: '+str(dtime), 'WARN') 00113 casalog.post( 'Total time used so far: ' + str(totaltime), 'WARN') 00114 casalog.post( 'Step Time used (s) Fraction of total time (percent) [description]', 'WARN') 00115 for i in range(0, len(steptime)): 00116 casalog.post( ' '+str(thesteps[i])+' '+str(steptime[i])+' '+str(steptime[i]/totaltime*100.) 00117 +' ['+step_title[thesteps[i]]+']', 'WARN') 00118 00119 # default tarfile name containing the input MSs 00120 mytarfile_name = 'NGC3256_Band3_UnCalibratedMSandTablesForReduction.tgz' 00121 # get the dataset name from the wrapper if possible 00122 mydict = locals() 00123 if mydict.has_key("tarfile_name"): 00124 if(mytarfile_name != mydict["tarfile_name"]): 00125 raise Exception, 'Wrong input file 1' 00126 00127 # data preparation 00128 mystep = 0 00129 if(mystep in thesteps): 00130 print 'Step ', mystep, step_title[mystep] 00131 00132 if not os.path.exists(mytarfile_name): 00133 raise Exception, 'Cannot find input file '+mytarfile_name 00134 for name in basename: 00135 os.system('rm -rf '+name+'.ms') 00136 00137 os.system('tar xvzf '+mytarfile_name) 00138 for name in basename: 00139 os.system('mv NGC3256_Band3_UnCalibratedMSandTablesForReduction/'+name+'.ms .') 00140 00141 timing() 00142 00143 # generate caltables 00144 mystep = 1 00145 if(mystep in thesteps): 00146 print 'Step ', mystep, step_title[mystep] 00147 00148 for name in basename: 00149 flagmanager(vis = name+'.ms', mode = 'save', versionname = 'Original') 00150 os.system('rm -rf cal-tsys_'+name+'.calnew cal-*'+name+'.Wnew') 00151 00152 gencal(vis=name+'.ms', 00153 caltype='tsys', 00154 caltable='cal-tsys_'+name+'.calnew') 00155 00156 wvrgcal(vis=name+'.ms', caltable='cal-'+name+'.Wnew', 00157 toffset=-1, segsource=True, 00158 tie=["Titan,1037-295,NGC3256"], statsource="1037-295", 00159 smooth=3) 00160 00161 if makeplots: 00162 for spw in ['1','3','5','7']: 00163 for name in basename: 00164 plotcal(caltable='cal-tsys_'+name+'.calnew', xaxis='freq', yaxis='amp', 00165 spw=spw, subplot=421, overplot=False, 00166 iteration='antenna', plotrange=[0, 0, 40, 180], plotsymbol='.', 00167 figfile='cal-tsys_per_spw_'+spw+'_'+name+'.png') 00168 00169 timing() 00170 00171 # Apriori flagging 00172 mystep = 2 00173 if(mystep in thesteps): 00174 print 'Step ', mystep, step_title[mystep] 00175 00176 for name in basename: 00177 flagmanager(vis = name+'.ms', mode = 'restore', versionname = 'Original') 00178 00179 flagdata(vis=name+'.ms', flagbackup = F, mode = 'shadow') 00180 flagdata(vis=name+'.ms',mode='manual', autocorr=True) 00181 flagdata(vis=name+'.ms', mode='manual', flagbackup = F, intent='*POINTING*') 00182 flagdata(vis=name+'.ms', mode='manual', flagbackup = F, intent='*ATMOSPHERE*') 00183 00184 flagmanager(vis = name+'.ms', mode = 'save', versionname = 'Apriori') 00185 00186 timing() 00187 00188 # Delay Correction for Antenna DV07 00189 mystep = 3 00190 if(mystep in thesteps): 00191 print 'Step ', mystep, step_title[mystep] 00192 00193 for i in range(3): # loop over the first three ms's 00194 name=basename[i] 00195 00196 os.system('rm -rf cal-'+name+'_del.Knew') 00197 gencal(vis=name+'.ms', caltable='cal-'+name+'_del.Knew', 00198 caltype='sbd', antenna='DV07', pol='X,Y', spw='1,3,5,7', 00199 parameter=[1.00, 1.10, -3.0, -3.0, -3.05, -3.05, -3.05, -3.05]) 00200 00201 timing() 00202 00203 # Applycal of delay and wvr corrections 00204 mystep = 4 00205 if(mystep in thesteps): 00206 print 'Step ', mystep, step_title[mystep] 00207 for i in range(3): # loop over the first three data sets 00208 name=basename[i] 00209 applycal(vis=name+'.ms', flagbackup=F, spw='1,3,5,7', 00210 interp=['nearest','nearest'], gaintable=['cal-'+name+'_del.Knew', 'cal-'+name+'.Wnew']) 00211 00212 for i in range(3,6): # loop over the last three data sets 00213 name=basename[i] 00214 applycal(vis=name+'.ms', flagbackup=F, spw='1,3,5,7', 00215 interp='nearest', gaintable='cal-'+name+'.Wnew') 00216 00217 for name in basename: 00218 os.system('rm -rf '+name+'_K_WVR.ms*') 00219 split(vis=name+'.ms', outputvis=name+'_K_WVR.ms', 00220 datacolumn='corrected', spw='0~7') 00221 00222 flagmanager(vis = name+'_K_WVR.ms', mode = 'save', versionname = 'Original') 00223 00224 timing() 00225 00226 # Applycal of tsys corrections, split out corrected data 00227 mystep = 5 00228 if(mystep in thesteps): 00229 print 'Step ', mystep, step_title[mystep] 00230 00231 flagmanager(vis ='uid___A002_X1d54a1_X174_K_WVR.ms', mode = 'restore', versionname = 'Original') 00232 00233 flagdata(vis='uid___A002_X1d54a1_X174_K_WVR.ms', mode='manual', 00234 antenna='DV04', flagbackup = T, scan='4,5,9', spw='7') 00235 00236 for name in basename: 00237 00238 for field in ['Titan','1037*','NGC*']: 00239 applycal(vis=name+'_K_WVR.ms', spw='1,3,5,7', 00240 flagbackup=F, field=field, gainfield=field, 00241 interp='linear,spline', gaintable=['cal-tsys_'+name+'.calnew']) 00242 00243 os.system('rm -rf '+name+'_line.ms*') 00244 split(vis=name+'_K_WVR.ms', outputvis=name+'_line.ms', 00245 datacolumn='corrected', spw='1,3,5,7') 00246 00247 00248 timing() 00249 00250 # concatenate the MSs 00251 mystep = 6 00252 if(mystep in thesteps): 00253 print 'Step ', mystep, step_title[mystep] 00254 00255 comvis=[] 00256 for name in basename: 00257 comvis.append(name+'_line.ms') 00258 00259 os.system('rm -rf ngc3256_line.ms*') 00260 concat(vis=comvis, concatvis='ngc3256_line.ms') 00261 00262 flagmanager(vis ='ngc3256_line.ms', mode = 'save', versionname = 'Original') 00263 00264 timing() 00265 00266 # fixplanets and flag the concatenated data 00267 mystep = 7 00268 if(mystep in thesteps): 00269 print 'Step ', mystep, step_title[mystep] 00270 00271 flagmanager(vis ='ngc3256_line.ms', mode = 'restore', versionname = 'Original') 00272 00273 flagdata(vis='ngc3256_line.ms', flagbackup=T, spw='*:0~16,*:125~127') 00274 00275 flagdata(vis = 'ngc3256_line.ms', flagbackup = T, 00276 timerange='>2011/04/16/12:00:00', field='Titan') 00277 00278 fixplanets(vis='ngc3256_line.ms', field='Titan', fixuvw=True) 00279 00280 flagdata(vis='ngc3256_line.ms', flagbackup=T, spw='3', 00281 correlation='YY', mode='manual', 00282 antenna='DV07', timerange='') 00283 00284 flagdata(vis='ngc3256_line.ms', flagbackup=T, spw='3', 00285 correlation='YY', mode='manual', 00286 antenna='DV08', timerange='>2011/04/17/03:00:00') 00287 00288 flagdata(vis='ngc3256_line.ms', flagbackup=T, spw='0', 00289 correlation='', mode='manual', 00290 antenna='PM03', timerange='2011/04/17/02:15:00~02:15:50') 00291 00292 flagdata(vis='ngc3256_line.ms', flagbackup=T, spw='2,3', 00293 correlation='', mode='manual', 00294 antenna='PM03', timerange='2011/04/16/04:13:50~04:18:00') 00295 00296 flagdata(vis='ngc3256_line.ms', flagbackup=T, spw='', 00297 correlation='', mode='manual', 00298 antenna='PM03&DV10', timerange='>2011/04/16/15:00:00') 00299 00300 timing() 00301 00302 # Fast phase-only gaincal on the bandpass calibrator 00303 mystep = 8 00304 if(mystep in thesteps): 00305 print 'Step ', mystep, step_title[mystep] 00306 00307 os.system('rm -rf cal-ngc3256.G1') 00308 gaincal(vis='ngc3256_line.ms', caltable='cal-ngc3256.G1', spw='*:40~80', field='1037*', 00309 selectdata=T, solint='int', refant=therefant, calmode='p') 00310 00311 if makeplots: 00312 plotcal(caltable = 'cal-ngc3256.G1', xaxis = 'time', yaxis = 'phase', 00313 poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw', 00314 figfile='cal-phase_vs_time_XX.G1.png', subplot = 221) 00315 plotcal(caltable = 'cal-ngc3256.G1', xaxis = 'time', yaxis = 'phase', 00316 poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration = 'spw', 00317 figfile='cal-phase_vs_time_YY.G1.png', subplot = 221) 00318 timing() 00319 00320 # Bandpass calibration 00321 mystep = 9 00322 if(mystep in thesteps): 00323 print 'Step ', mystep, step_title[mystep] 00324 00325 os.system('rm -rf cal-ngc3256.B1') 00326 bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1', 00327 gaintable = 'cal-ngc3256.G1', timerange='<2011/04/16/15:00:00', 00328 field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan', 00329 bandtype='B', fillgaps=1, refant = therefant, solnorm = T) 00330 00331 bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1', 00332 gaintable = 'cal-ngc3256.G1', timerange='>2011/04/16/15:00:00', 00333 field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan', 00334 bandtype='B', fillgaps=1, refant = therefant, solnorm = T, append=True) 00335 00336 if makeplots: 00337 plotcal(caltable = 'cal-ngc3256.B1', xaxis='freq', yaxis='amp', spw='', 00338 subplot=211, overplot=False, 00339 figfile='bandpass.B1.png', plotsymbol='.', timerange='') 00340 00341 timing() 00342 00343 # Set fluxscale on Titan 00344 mystep = 10 00345 if(mystep in thesteps): 00346 print 'Step ', mystep, step_title[mystep] 00347 00348 setjy(vis='ngc3256_line.ms', field='Titan', standard='Butler-JPL-Horizons 2010', 00349 spw='0,1,2,3', scalebychan=False, usescratch=False) 00350 00351 timing() 00352 00353 # Amplitude and Phase gaincal 00354 mystep = 11 00355 if(mystep in thesteps): 00356 print 'Step ', mystep, step_title[mystep] 00357 00358 gaincal(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.G2', spw = 00359 '*:16~112', field = '1037*,Titan', minsnr=1.0, 00360 solint= 'inf', selectdata=T, solnorm=False, refant = therefant, 00361 gaintable = 'cal-ngc3256.B1', calmode = 'ap') 00362 00363 if makeplots: 00364 plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'phase', 00365 poln='X', plotsymbol='o', plotrange = [0,0,-180,180], iteration= 'spw', 00366 figfile='cal-phase_vs_time_XX.G2.png', subplot = 221) 00367 plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'phase', 00368 poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], iteration= 'spw', 00369 figfile='cal-phase_vs_time_YY.G2.png', subplot = 221) 00370 plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'amp', 00371 poln='X', plotsymbol='o', plotrange = [], iteration = 'spw', 00372 figfile='cal-amp_vs_time_XX.G2.png', subplot = 221) 00373 plotcal(caltable = 'cal-ngc3256.G2', xaxis = 'time', yaxis = 'amp', 00374 poln='Y', plotsymbol='o', plotrange = [], iteration = 'spw', 00375 figfile='cal-amp_vs_time_YY.G2.png', subplot = 221) 00376 00377 timing() 00378 00379 # Adjust absolute fluxscale 00380 mystep = 12 00381 if(mystep in thesteps): 00382 print 'Step ', mystep, step_title[mystep] 00383 00384 fluxscale(vis="ngc3256_line.ms", caltable="cal-ngc3256.G2", 00385 fluxtable="cal-ngc3256.G2.flux", reference="Titan", 00386 transfer="1037*", refspwmap=[0,1,1,1]) 00387 00388 timing() 00389 00390 # Apply bandpass and gain calibration 00391 mystep = 13 00392 if(mystep in thesteps): 00393 print 'Step ', mystep, step_title[mystep] 00394 00395 applycal(vis='ngc3256_line.ms', flagbackup=F, field='NGC*,1037*', 00396 interp=['nearest','nearest'], gainfield = ['1037*', '1037*'], 00397 gaintable=['cal-ngc3256.G2.flux', 'cal-ngc3256.B1']) 00398 00399 flagmanager(vis = name+'.ms', mode = 'save', versionname = 'step12') 00400 00401 timing() 00402 00403 # Final flagging 00404 mystep = 14 00405 if(mystep in thesteps): 00406 print 'Step ', mystep, step_title[mystep] 00407 00408 flagmanager(vis = name+'.ms', mode = 'restore', versionname = 'step12') 00409 00410 flagdata(vis='ngc3256_line.ms', mode='manual', 00411 timerange='2011/04/16/04:13:35~04:13:45', flagbackup = T) 00412 00413 flagdata(vis='ngc3256_line.ms', mode='manual', 00414 timerange='2011/04/16/05:21:13~05:21:19', flagbackup = T) 00415 00416 flagdata(vis='ngc3256_line.ms', mode='manual', 00417 timerange='2011/04/16/04:16:40~04:16:49', flagbackup = T) 00418 00419 flagdata(vis='ngc3256_line.ms', mode='manual', 00420 timerange='2011/04/16/04:14:00~04:17:10', antenna='PM03', flagbackup = T) 00421 00422 flagdata(vis='ngc3256_line.ms', mode='manual', 00423 timerange='2011/04/17/00:35:30~01:20:20', antenna='DV04', spw='3', flagbackup = T) 00424 00425 flagmanager(vis = name+'.ms', mode = 'save', versionname = 'step13') 00426 00427 timing() 00428 00429 # Final calibration 00430 mystep = 15 00431 if(mystep in thesteps): 00432 print 'Step ', mystep, step_title[mystep] 00433 00434 delmod('ngc3256_line.ms') 00435 00436 gaincal(vis='ngc3256_line.ms', caltable='cal-ngc3256.G1n', spw='*:40~80', field='1037*', 00437 selectdata=T, solint='int', refant=therefant, calmode='p') 00438 00439 bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1n', 00440 gaintable = 'cal-ngc3256.G1n', timerange='<2011/04/16/15:00:00', 00441 field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan', 00442 bandtype='B', fillgaps=1, refant = therefant, solnorm = T) 00443 00444 bandpass(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.B1n', 00445 gaintable = 'cal-ngc3256.G1n', timerange='>2011/04/16/15:00:00', 00446 field = '1037*', minblperant=3, minsnr=2, solint='inf', combine='scan', 00447 bandtype='B', fillgaps=1, refant = therefant, solnorm = T, append=True) 00448 00449 setjy(vis='ngc3256_line.ms', field='Titan', standard='Butler-JPL-Horizons 2010', 00450 spw='0,1,2,3', usescratch=False) 00451 00452 gaincal(vis = 'ngc3256_line.ms', caltable = 'cal-ngc3256.G2n', spw = 00453 '*:16~112', field = '1037*,Titan', minsnr=1.0, 00454 solint= 'inf', selectdata=T, solnorm=False, refant = therefant, 00455 gaintable = 'cal-ngc3256.B1n', calmode = 'ap') 00456 00457 fluxscale( vis="ngc3256_line.ms", caltable="cal-ngc3256.G2n", 00458 fluxtable="cal-ngc3256.G2n.flux", reference="Titan", 00459 transfer="1037*", refspwmap=[0,1,1,1]) 00460 00461 applycal( vis='ngc3256_line.ms', flagbackup=F, field='NGC*,1037*', 00462 interp=['nearest','nearest'], gainfield = ['1037*', '1037*'], 00463 gaintable=['cal-ngc3256.G2n.flux', 'cal-ngc3256.B1n']) 00464 00465 timing() 00466 00467 # Image the phase calibrator 00468 mystep = 16 00469 if(mystep in thesteps): 00470 print 'Step ', mystep, step_title[mystep] 00471 00472 os.system('rm -rf result-phasecal_cont*') 00473 clean(vis='ngc3256_line.ms', imagename='result-phasecal_cont', field='1037*', 00474 spw='*:20~120', selectdata=T, mode='mfs', niter=500, 00475 gain=0.1, threshold='0.75mJy', psfmode='hogbom', 00476 interactive=False, mask=[62, 62, 67, 67], imsize=128, 00477 cell='1arcsec', weighting='briggs', robust=0.0, nterms=2, 00478 usescratch=False) 00479 00480 calstat=imstat(imagename='result-phasecal_cont.image.tt0', region='', box='85,8,120,120') 00481 rmspcal=(calstat['rms'][0]) 00482 print '>> rms in phase calibrator image: '+str(rmspcal) 00483 calstat=imstat(imagename='result-phasecal_cont.image.tt0', region='') 00484 peakpcal=(calstat['max'][0]) 00485 print '>> Peak in phase calibrator image: '+str(peakpcal) 00486 print '>> Dynamic range in phase calibrator image: '+str(peakpcal/rmspcal) 00487 00488 if makeplots: 00489 imview(raster={'file': 'result-phasecal_cont.image.tt0', 'colorwedge':T, 00490 'range':[-0.004, 0.250], 'scaling':-2.5, 'colormap':'Rainbow 2'}, 00491 out='result-phasecal_map.png', zoom=1) 00492 00493 imview(raster={'file': 'result-phasecal_cont.image.tt0', 'colorwedge':T, 00494 'colormap':'Rainbow 2'}, 00495 out='result-phasecal_map-lin.png', zoom=1) 00496 00497 timing() 00498 00499 # Apply calibration to flux calibrator 00500 mystep = 17 00501 if(mystep in thesteps): 00502 print 'Step ', mystep, step_title[mystep] 00503 00504 applycal(vis='ngc3256_line.ms', flagbackup=F, field='Titan', 00505 interp=['nearest', 'nearest'], gainfield = ['Titan', '1037*'], 00506 gaintable=['cal-ngc3256.G2.flux', 'cal-ngc3256.B1']) 00507 00508 timing() 00509 00510 # Image flux calibrator 00511 mystep = 18 00512 if(mystep in thesteps): 00513 print 'Step ', mystep, step_title[mystep] 00514 00515 os.system('rm -rf result-ampcal_cont*') 00516 clean(vis='ngc3256_line.ms', imagename='result-ampcal_cont', 00517 field='Titan', spw='0:20~120,1:20~120', mode='mfs', niter=200, 00518 threshold='5mJy', psfmode='hogbom', mask=[62, 62, 67, 67], imsize=128, 00519 cell='1arcsec', weighting='briggs', robust=0.0, 00520 usescratch=False) 00521 00522 calstat=imstat(imagename="result-ampcal_cont.image",region="",box="85,8,120,120") 00523 rmstitan=(calstat['rms'][0]) 00524 print ">> rms in amp calibrator image: "+str(rmstitan) 00525 calstat=imstat(imagename="result-ampcal_cont.image",region="") 00526 peaktitan=(calstat['max'][0]) 00527 print ">> Peak in amp calibrator image: "+str(peaktitan) 00528 print ">> Dynamic range in amp calibrator image: "+str(peaktitan/rmstitan) 00529 00530 if makeplots: 00531 imview(raster={'file': 'result-ampcal_cont.image', 'colorwedge':T, 00532 'range':[-0.02, 0.250], 'scaling':-1.5, 'colormap':'Rainbow 2'}, 00533 out='result-ampcal_map.png', zoom=1) 00534 00535 imview(raster={'file': 'result-ampcal_cont.image', 'colorwedge':T, 00536 'colormap':'Rainbow 2'}, 00537 out='result-ampcal_map-lin.png', zoom=1) 00538 00539 timing() 00540 00541 # split out the calibrated target data 00542 mystep = 19 00543 if(mystep in thesteps): 00544 print 'Step ', mystep, step_title[mystep] 00545 00546 os.system('rm -rf ngc3256_line_target.ms*') 00547 split(vis='ngc3256_line.ms', outputvis='ngc3256_line_target.ms', 00548 field='NGC*') 00549 00550 # Image the NGC3256 continuum 00551 mystep = 20 00552 if(mystep in thesteps): 00553 print 'Step ', mystep, step_title[mystep] 00554 00555 os.system('rm -rf result-ngc3256_cont*') 00556 clean( vis='ngc3256_line_target.ms', imagename='result-ngc3256_cont', 00557 spw='0:20~53;71~120,1:70~120,2:20~120,3:20~120', psfmode='hogbom', 00558 mode='mfs', niter=500, threshold='0.3mJy', mask=[42,43,59,60], 00559 imsize=100, cell='1arcsec', weighting='briggs', robust=0.0, 00560 interactive=False, usescratch=False) 00561 00562 calstat=imstat(imagename='result-ngc3256_cont.image', region='', box='10,10,90,35') 00563 rmscont=(calstat['rms'][0]) 00564 print '>> rms in continuum image: '+str(rmscont) 00565 calstat=imstat(imagename='result-ngc3256_cont.image', region='') 00566 peakcont=(calstat['max'][0]) 00567 print '>> Peak in continuum image: '+str(peakcont) 00568 print '>> Dynamic range in continuum image: '+str(peakcont/rmscont) 00569 00570 if makeplots: 00571 imview(raster={'file': 'result-ngc3256_cont.image', 'colorwedge':T, 00572 'range':[-0.001, 0.009], 'scaling':0, 'colormap':'Rainbow 2'}, 00573 out='result-ngc3256_cont.png', zoom=2) 00574 00575 timing() 00576 00577 # Slow phase selfcal 00578 mystep = 21 00579 if(mystep in thesteps): 00580 print 'Step ', mystep, step_title[mystep] 00581 00582 os.system('rm -rf cal-ngc3256_cont_30m.Gp') 00583 gaincal(vis='ngc3256_line_target.ms', field='NGC*', 00584 caltable='cal-ngc3256_cont_30m.Gp', 00585 spw='0:20~53;71~120,1:70~120,2:20~120,3:20~120', 00586 solint='1800s', refant='DV07', calmode='p', 00587 minblperant=3) 00588 00589 if makeplots: 00590 plotcal(caltable = 'cal-ngc3256_cont_30m.Gp', 00591 xaxis = 'time', yaxis = 'phase', 00592 poln='X', plotsymbol='o', plotrange = [0,0,-180,180], 00593 iteration = 'spw', figfile='cal-phase_vs_time_XX_30_Gp.png', 00594 subplot = 221) 00595 00596 plotcal(caltable = 'cal-ngc3256_cont_30m.Gp', 00597 xaxis = 'time', yaxis = 'phase', 00598 poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], 00599 iteration = 'spw', figfile='cal-phase_vs_time_YY_30_Gp.png', 00600 subplot = 221) 00601 00602 applycal(vis='ngc3256_line_target.ms', interp='linear', 00603 gaintable='cal-ngc3256_cont_30m.Gp') 00604 00605 timing() 00606 00607 # Image the NGC3256 continuum again after selfcal 00608 mystep = 22 00609 if(mystep in thesteps): 00610 print 'Step ', mystep, step_title[mystep] 00611 00612 os.system('rm -rf result-ngc3256_cont_sc1*') 00613 clean(vis='ngc3256_line_target.ms', imagename='result-ngc3256_cont_sc1', 00614 spw='0:20~53;71~120,1:70~120,2:20~120,3:20~120', psfmode='hogbom', 00615 mode='mfs', niter=500, threshold='0.13mJy', mask=[42,43,59,60], 00616 imsize=100, cell='1arcsec', weighting='briggs', robust=0.0, 00617 interactive=False, usescratch=False) 00618 00619 calstat=imstat(imagename='result-ngc3256_cont_sc1.image', region='', 00620 box='10,10,90,35') 00621 rmscontsc=(calstat['rms'][0]) 00622 print '>> rms in continuum image: '+str(rmscontsc) 00623 calstat=imstat(imagename='result-ngc3256_cont_sc1.image', region='') 00624 peakcontsc=(calstat['max'][0]) 00625 print '>> Peak in continuum image: '+str(peakcontsc) 00626 print '>> Dynamic range in continuum image: '+str(peakcontsc/rmscontsc) 00627 00628 if makeplots: 00629 imview(raster={'file': 'result-ngc3256_cont_sc1.image', 'colorwedge':T, 00630 'range':[-0.001, 0.009], 'scaling':0, 'colormap':'Rainbow 2'}, 00631 out='result-ngc3256_cont_sc1.png', zoom=2) 00632 00633 timing() 00634 00635 # Determine and subtract continuum 00636 mystep = 23 00637 if(mystep in thesteps): 00638 print 'Step ', mystep, step_title[mystep] 00639 00640 uvcontsub(vis = 'ngc3256_line_target.ms', 00641 fitspw='0:20~53;71~120,1:70~120,2:20~120,3:20~120', solint ='int', 00642 fitorder = 1, 00643 combine='spw') 00644 00645 timing() 00646 00647 # Clean the NGC3256 line cube 00648 mystep = 24 00649 if(mystep in thesteps): 00650 print 'Step ', mystep, step_title[mystep] 00651 00652 os.system('rm -rf result-ngc3256_line_CO.*') 00653 clean(vis='ngc3256_line_target.ms.contsub', imagename='result-ngc3256_line_CO', 00654 spw='0:38~87', mode='channel', start='', nchan=50, width='', 00655 psfmode='hogbom', outframe='LSRK', restfreq='115.271201800GHz', 00656 mask=[53,50,87,83], niter=500, interactive=False, imsize=128, cell='1arcsec', 00657 weighting='briggs', robust=0.0, threshold='5mJy', 00658 usescratch=False) 00659 00660 timing() 00661 00662 # Evaluate the NGC3256 line cube 00663 mystep = 25 00664 if(mystep in thesteps): 00665 print 'Step ', mystep, step_title[mystep] 00666 00667 os.system('rm -rf myresults.tbl') 00668 slsearch(outfile='myresults.tbl', freqrange = [84,116], species=['COv=0']) 00669 sl.open('myresults.tbl') 00670 sl.list() 00671 00672 tb.open('myresults.tbl') 00673 restfreq=tb.getcol('FREQUENCY')[0] 00674 tb.close() 00675 00676 os.system('rm -rf result-ngc3256_CO1-0.mom0*') 00677 immoments(imagename='result-ngc3256_line_CO.image', moments=[0], 00678 chans='15~34', box='38,38,90,90', axis='spectral', 00679 includepix=[0.02, 10000], outfile='result-ngc3256_CO1-0.mom0') 00680 00681 os.system('rm -rf result-ngc3256_CO1-0.mom1*') 00682 immoments(imagename='result-ngc3256_line_CO.image', moments=[1], 00683 chans='15~34', box='38,38,90,90', axis='spectral', 00684 includepix=[0.045, 10000], outfile='result-ngc3256_CO1-0.mom1') 00685 00686 os.system('rm -rf result-ngc3256_CO1-0.mom2*') 00687 immoments(imagename='result-ngc3256_line_CO.image', moments=[2], 00688 chans='5~44', box='38,38,90,90', axis='spectral', 00689 includepix=[0.035, 10000], outfile='result-ngc3256_CO1-0.mom2') 00690 00691 if makeplots: 00692 imview(contour={'file': 'result-ngc3256_CO1-0.mom0','levels': 00693 [5,10,20,40,80,160],'base':0,'unit':1}, 00694 raster={'file': 'result-ngc3256_CO1-0.mom1','range': [2630,2920], 00695 'colorwedge':T, 'colormap': 'Rainbow 2'}, out='result-CO_velfield.png') 00696 00697 imview(contour={'file': 'result-ngc3256_CO1-0.mom1','levels': 00698 [2650,2700,2750,2800,2850,2900],'base':0,'unit':1}, 00699 raster={'file': 'result-ngc3256_CO1-0.mom0', 'colorwedge':T, 00700 'colormap': 'Rainbow 2','scaling':-1.0,'range': [0.8,250]}, 00701 out='result-CO_map.png') 00702 00703 imview(contour={'file': 'result-ngc3256_CO1-0.mom2','levels': 00704 [20,30,40,50,60],'base':0,'unit':1}, 00705 raster={'file': 'result-ngc3256_CO1-0.mom2', 'colorwedge':T, 00706 'colormap': 'Greyscale 1','scaling':-1.0,'range': [0,74]}, 00707 out='result-CO_dispersion.png') 00708 00709 timing() 00710 00711 # Clean the NGC3256 CNhi line cube 00712 mystep = 26 00713 if(mystep in thesteps): 00714 print 'Step ', mystep, step_title[mystep] 00715 00716 os.system('rm -rf result-ngc3256_line_CNhi.*') 00717 clean(vis='ngc3256_line_target.ms', imagename='result-ngc3256_line_CNhi', 00718 outframe='LSRK', spw='1:50~76', start='', nchan=27, width='', 00719 restfreq='113.48812GHz', selectdata=T, mode='channel', 00720 niter=500, gain=0.1, psfmode='hogbom', mask=[53,50,87,83], 00721 interactive=False, imsize=128, cell='1arcsec', 00722 weighting='briggs', robust=0.0, threshold='2mJy', 00723 usescratch=False) 00724 00725 timing() 00726 00727 # Evaluate the NGC3256 CNhi line cube 00728 mystep = 27 00729 if(mystep in thesteps): 00730 print 'Step ', mystep, step_title[mystep] 00731 00732 os.system('rm -rf result-ngc3256_CNhi.mom.*') 00733 immoments( imagename='result-ngc3256_line_CNhi.image', moments=[0,1], 00734 chans='5~18', axis='spectral', box='38,38,90,90', 00735 includepix=[0.005, 10000], outfile='result-ngc3256_CNhi.mom') 00736 if makeplots: 00737 imview(contour={'file': 'result-ngc3256_CNhi.mom.integrated','range': []}, 00738 raster={'file': 'result-ngc3256_CNhi.mom.weighted_coord', 00739 'range': [2630,2920],'colorwedge':T, 00740 'colormap': 'Rainbow 2'}, out='result-CNhi_velfield.png') 00741 imview(contour={'file': 'result-ngc3256_CNhi.mom.weighted_coord','levels': 00742 [2650,2700,2750,2800,2850,2900],'base':0,'unit':1}, 00743 raster={'file': 'result-ngc3256_CNhi.mom.integrated','colorwedge':T, 00744 'colormap': 'Rainbow 2'}, out='result-CNhi_map.png') 00745 00746 timing() 00747 00748 # Clean the NGC3256 CNlo line cube 00749 mystep = 28 00750 if(mystep in thesteps): 00751 print 'Step ', mystep, step_title[mystep] 00752 00753 os.system('rm -rf result-ngc3256_line_CNlo.*') 00754 clean( vis='ngc3256_line_target.ms', imagename='result-ngc3256_line_CNlo', 00755 outframe='LSRK', spw='1:29~54', start='', nchan=26, width='', 00756 restfreq='113.17049GHz', selectdata=T, mode='channel', 00757 niter=300, gain=0.1, psfmode='hogbom', mask=[53,50,87,83], 00758 interactive=False, imsize=128, cell='1arcsec', 00759 weighting='briggs', robust=0.0, threshold='2mJy', 00760 usescratch=False) 00761 00762 timing() 00763 00764 # Evaluate the NGC3256 CNlo line cube 00765 mystep = 29 00766 if(mystep in thesteps): 00767 print 'Step ', mystep, step_title[mystep] 00768 00769 os.system('rm -rf result-ngc3256_CNlo.mom.*') 00770 immoments( imagename='result-ngc3256_line_CNlo.image', moments=[0,1], 00771 chans='8~18', axis='spectral', box='38,38,90,90', 00772 includepix=[0.003, 10000], outfile='result-ngc3256_CNlo.mom') 00773 00774 if makeplots: 00775 imview(contour={'file': 'result-ngc3256_CNlo.mom.integrated','range': []}, 00776 raster={'file': 'result-ngc3256_CNlo.mom.weighted_coord', 00777 'range': [2630,2920],'colorwedge':T, 00778 'colormap': 'Rainbow 2'}, out='result-CNlo_velfield.png') 00779 00780 imview(contour={'file': 'result-ngc3256_CNlo.mom.weighted_coord','levels': 00781 [2650,2700,2750,2800,2850,2900],'base':0,'unit':1}, 00782 raster={'file': 'result-ngc3256_CNlo.mom.integrated','colorwedge':T, 00783 'colormap': 'Rainbow 2'}, out='result-CNlo_map.png') 00784 00785 timing() 00786 00787 # Verify regression results 00788 mystep = 30 00789 if(mystep in thesteps): 00790 print 'Step ', mystep, step_title[mystep] 00791 00792 # reference values obtained with the NGC3256 Band 3 CASA guide for CASA 3.3 (9 May 2012, DP) 00793 refrmspcal33 = 0.000556594866794 00794 refpeakpcal33 = 1.86264276505 00795 refrmstitan33 = 0.00490083405748 00796 refpeaktitan33 = 0.361426800489 00797 refrmscont33 = 0.000324592343532 00798 refpeakcont33 = 0.00705300131813 00799 refrmscontsc33 = 8.03611983429e-05 00800 refpeakcontsc33 = 0.00976239796728 00801 00802 # reference values obtained with this script using CASA stable r19569 (10 May 2012, DP) 00803 #refrmspcal = 0.000612945703324 00804 #refpeakpcal = 1.86513638496 00805 #refrmstitan = 0.00491443555802 00806 #refpeaktitan = 0.361274629831 00807 #refrmscont = 0.00032463966636 00808 #refpeakcont = 0.00699360202998 00809 #refrmscontsc = 7.81473427196e-05 00810 #refpeakcontsc = 0.00983608141541 00811 00812 # reference values obtained with this script using CASA stable r21621 (15 Oct 2012, DP) 00813 refrmspcal = 0.000963084050454 00814 refpeakpcal = 1.87274956703 00815 refrmstitan = 0.00491269072518 00816 refpeaktitan = 0.361274719238 00817 refrmscont = 0.000305994151859 00818 refpeakcont = 0.00709502119571 00819 refrmscontsc = 7.5960662798e-05 00820 refpeakcontsc = 0.00963686406612 00821 00822 devrmspcal = abs(rmspcal-refrmspcal)/refrmspcal*100. 00823 devpeakpcal = abs(peakpcal-refpeakpcal)/refpeakpcal*100. 00824 devrmstitan = abs(rmstitan-refrmstitan)/refrmstitan*100. 00825 devpeaktitan = abs(peaktitan-refpeaktitan)/refpeaktitan*100. 00826 devrmscont = abs(rmscont-refrmscont)/refrmscont*100. 00827 devpeakcont = abs(peakcont-refpeakcont)/refpeakcont*100. 00828 devrmscontsc = abs(rmscontsc-refrmscontsc)/refrmscontsc*100. 00829 devpeakcontsc = abs(peakcontsc-refpeakcontsc)/refpeakcontsc*100. 00830 00831 casalog.origin('SUMMARY') 00832 casalog.post("\n***** Peak and RMS of the images of the primary phase calibrator, Titan, and the Target *****") 00833 casalog.post( "Field, Peak (expectation, expectation CASA 3.3), RMS (expectation, expectation CASA 3.3)") 00834 casalog.post( "------------------------------------------------------------------------------------------") 00835 casalog.post("Phasecal, "+str(peakpcal)+" ("+str(refpeakpcal)+", "+str(refpeakpcal33)+"), " 00836 +str(rmspcal)+" ("+str(refrmspcal)+", "+str(refrmspcal33)+")") 00837 casalog.post( "------------------------------------------------------------------------------------------") 00838 casalog.post("Titan, "+str(peaktitan)+" ("+str(refpeaktitan)+", "+str(refpeaktitan33)+"), " 00839 +str(rmstitan)+" ("+str(refrmstitan)+", "+str(refrmstitan33)+")") 00840 casalog.post( "------------------------------------------------------------------------------------------") 00841 casalog.post("NGC3256 continuum, "+str(peakcont)+" ("+str(refpeakcont)+", "+str(refpeakcont33)+"), " 00842 +str(rmscont)+" ("+str(refrmscont)+", "+str(refrmscont33)+")") 00843 casalog.post( "------------------------------------------------------------------------------------------") 00844 casalog.post("NGC3256 selfcaled continuum, "+str(peakcontsc)+" ("+str(refpeakcontsc)+", "+str(refpeakcontsc33)+"), " 00845 +str(rmscontsc)+" ("+str(refrmscontsc)+", "+str(refrmscontsc33)+")") 00846 casalog.post( "------------------------------------------------------------------------------------------") 00847 00848 print rmspcal, refrmspcal 00849 print peakpcal, refpeakpcal 00850 print rmstitan, refrmstitan 00851 print peaktitan, refpeaktitan 00852 print rmscont, refrmscont 00853 print peakcont, refpeakcont 00854 print rmscontsc, refrmscontsc 00855 print peakcontsc, refpeakcontsc 00856 00857 passed = True 00858 00859 if devpeakpcal>0.5: 00860 casalog.post( 'ERROR: Peak in primary phase calibrator image deviates from expectation by '+str(devpeakpcal)+' percent.', 'WARN') 00861 passed = False 00862 if devrmspcal>0.5: 00863 casalog.post( 'ERROR: RMS in primary phase calibrator image deviates from expectation by '+str(devrmspcal)+' percent.', 'WARN') 00864 passed = False 00865 if devpeaktitan>0.5: 00866 casalog.post( 'ERROR: Peak in image of Titan deviates from expectation by '+str(devpeaktitan)+' percent.', 'WARN') 00867 passed = False 00868 if devrmstitan>0.5: 00869 casalog.post( 'ERROR: RMS in image of Titan deviates from expectation by '+str(devrmstitan)+' percent.', 'WARN') 00870 passed = False 00871 if devpeakcont>0.5: 00872 casalog.post( 'ERROR: Peak in continuum image deviates from expectation by '+str(devpeakcont)+' percent.', 'WARN') 00873 passed = False 00874 if devrmscont>0.5: 00875 casalog.post( 'ERROR: RMS in continuum image deviates from expectation by '+str(devrmscont)+' percent.', 'WARN') 00876 passed = False 00877 if devpeakcontsc>0.5: 00878 casalog.post( 'ERROR: Peak in selfcal continuum image deviates from expectation by '+str(devpeakcontsc)+' percent.', 'WARN') 00879 passed = False 00880 if devrmscontsc>0.5: 00881 casalog.post( 'ERROR: RMS in selfcal continuum image deviates from expectation by '+str(devrmscontsc)+' percent.', 'WARN') 00882 passed = False 00883 00884 if not passed: 00885 casalog.post('Results are different from expectations by more than 0.5 percent.', 'WARN') 00886 else: 00887 casalog.post( "\nAll peak and RMS values within 0.5 percent of the expectation.") 00888 00889 pngfiles = ['cal-tsys_per_spw_1_uid___A002_X1d54a1_X174.png', 00890 'cal-tsys_per_spw_7_uid___A002_X1d5a20_X5.png', 00891 'cal-phase_vs_time_XX.G1.png', 00892 'cal-phase_vs_time_YY.G1.png', 00893 'bandpass.B1.png', 00894 'cal-phase_vs_time_XX.G2.png', 00895 'cal-phase_vs_time_YY.G2.png', 00896 'cal-amp_vs_time_XX.G2.png', 00897 'cal-amp_vs_time_YY.G2.png', 00898 'result-phasecal_map.png', 00899 'result-phasecal_map-lin.png', 00900 'result-ampcal_map.png', 00901 'result-ampcal_map-lin.png', 00902 'result-ngc3256_cont.png', 00903 'cal-phase_vs_time_XX_30_Gp.png', 00904 'cal-phase_vs_time_YY_30_Gp.png', 00905 'result-ngc3256_cont_sc1.png', 00906 'result-CO_velfield.png', 00907 'result-CO_map.png', 00908 'result-CO_dispersion.png', 00909 'result-CNhi_velfield.png', 00910 'result-CNhi_map.png', 00911 'result-CNlo_velfield.png', 00912 'result-CNlo_map.png'] 00913 00914 for pngfile in pngfiles: 00915 if not os.path.exists(pngfile): 00916 casalog.post('Plot file '+pngfile+' was not created!', 'WARN') 00917 passed = False 00918 00919 if not passed: 00920 raise Exception, 'Regression failed.' 00921 00922 print "Passed." 00923 00924 timing()
1.8.0
