|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ############################################################################# 00002 # $Id:$ 00003 # Test Name: # 00004 # alma-m100-analysis-regression.py # 00005 # An ALMA Science Verification Data Analysis Regression # 00006 # using observation of M100 from September 2011 # 00007 # # 00008 # Rationale for Inclusion: # 00009 # Need test of complete ALMA analysis chain # 00010 # # 00011 # Input data: # 00012 # two ASDMs # 00013 # two tsys cal tables # 00014 # the clean masks # 00015 # # 00016 ############################################################################# 00017 00018 # alma-m100-analysis.py 00019 # An ALMA Science Verification Data Analysis Regression 00020 # using observation of M100 from September 2011 00021 00022 #import analysisUtils as aU 00023 # for casa before release 3.4.0, aU is needed for the Tsys cal table 00024 # generation which we ommit in this regression by using canned Tsys cal tables 00025 # for 3.4.0, the tsys table can be generated without the help of aU 00026 00027 00028 step_title = { 0 : 'Data import', 00029 1 : 'Generate antenna position cal tables', 00030 2 : 'Generate tsys cal tables (if possible) or use canned ones', 00031 3 : 'Correct the Titan position', 00032 4 : 'Apriori flagging', 00033 5 : 'Generate WVR cal tables', 00034 6 : 'Generate delay calibration tables', 00035 7 : 'Apply antpos, wvr, tsys, and delay cal tables', 00036 8 : 'Split off non-wvr spws and save flags', 00037 9 : 'Flagging', 00038 10: 'Rebin to a reduced resolution of approx. 10 km/s', 00039 11: 'Fast phase-only gaincal for bandpass', 00040 12: 'Bandpass', 00041 13: 'Setjy', 00042 14: 'Fast phase-only gaincal', 00043 15: 'Slow phase-only gaincal', 00044 16: 'Slow amp and phase gaincal', 00045 17: 'Fluxscale', 00046 18: 'Applycal', 00047 19: 'Test image of the secondary phase cal', 00048 20: 'Test image of the primary phase cal', 00049 21: 'Test image of Titan', 00050 22: 'Split off calibrated M100 data', 00051 23: 'Concatenate M100 data', 00052 24: 'Average concatenated M100 data in time', 00053 25: 'Continuum image of M100', 00054 26: 'Determine and subtract continuum', 00055 27: 'Test image of central field', 00056 28: 'Clean line cube mosaic', 00057 29: 'Make moment maps', 00058 30: 'Verification of the regression results' 00059 } 00060 00061 # global defs 00062 basename=['X54','X220'] 00063 makeplots=False 00064 therefant = 'DV01' 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 00120 # default ASDM dataset name 00121 myasdm_dataset_name = "uid___A002_X2a5c2f_X54" 00122 myasdm_dataset2_name = "uid___A002_X2a5c2f_X220" 00123 mytsystable = 'cal-tsys_X54.fdm' 00124 mytsystable2 = 'cal-tsys_X220.fdm' 00125 00126 # get the dataset name from the wrapper if possible 00127 mydict = locals() 00128 if mydict.has_key("asdm_dataset_name"): 00129 if(myasdm_dataset_name != mydict["asdm_dataset_name"]): 00130 raise("Wrong input file 1") 00131 if mydict.has_key("asdm_dataset2_name"): 00132 if(myasdm_dataset2_name != mydict["asdm_dataset2_name"]): 00133 raise("Wrong input file 2") 00134 if mydict.has_key("tsys_table"): 00135 if(mytsystable != mydict["tsys_table"]): 00136 raise("Wrong input file 3") 00137 if mydict.has_key("tsys_table2"): 00138 if(mytsystable2 != mydict["tsys_table2"]): 00139 raise("Wrong input file 4") 00140 00141 00142 # data import 00143 mystep = 0 00144 if(mystep in thesteps): 00145 print 'Step ', mystep, step_title[mystep] 00146 00147 importasdm(asdm='uid___A002_X2a5c2f_X54',vis='X54.ms',asis='Stati* Anten*', overwrite=True) 00148 importasdm(asdm='uid___A002_X2a5c2f_X220',vis='X220.ms',asis='Stati* Anten*', overwrite=True) 00149 00150 if(makeplots): 00151 for name in basename: 00152 os.system('rm -rf plotants-'+name+'.png') 00153 plotants(name+'.ms', figfile='plotants-'+name+'.png') 00154 00155 timing() 00156 00157 00158 # Antenna positions 00159 mystep = 1 00160 if(mystep in thesteps): 00161 print 'Step ', mystep, step_title[mystep] 00162 00163 for name in basename: 00164 os.system('rm -rf cal-antpos_'+name) 00165 00166 gencal(vis = 'X54.ms', 00167 caltable = 'cal-antpos_X54', 00168 caltype = 'antpos', 00169 antenna = 'CM01,PM02,PM04,DV08,DV06,DV13,DV03,DV14', 00170 parameter = [-0.000379055076246,0.000910912511392,-0.000226045671848,7.8790821135e-05,0.000363811850548,-0.000224065035582,-0.000315962965482,0.000397399838991,-0.000581170175089,-6.35427422822e-05,0.00129573699087,-0.000625061802566,-0.000167516991496,0.000174060463905,-0.000417742878199,-0.00010875146836,0.000319179147482,-0.000588130671531,0.000142965465784,0.000455257482827,0.000168651808053,-0.00150847900659,0.00357818510383,0.000811365433037] 00171 ) 00172 gencal(vis = 'X220.ms', 00173 caltable = 'cal-antpos_X220', 00174 caltype = 'antpos', 00175 antenna = 'CM01,PM02,PM04,DV08,DV06,DV13,DV03,DV14', 00176 parameter = [-0.000379055076246,0.000910912511392,-0.000226045671848,7.8790821135e-05,0.000363811850548,-0.000224065035582,-0.000315962965482,0.000397399838991,-0.000581170175089,-6.35427422822e-05,0.00129573699087,-0.000625061802566,-0.000167516991496,0.000174060463905,-0.000417742878199,-0.00010875146836,0.000319179147482,-0.000588130671531,0.000142965465784,0.000455257482827,0.000168651808053,-0.00150847900659,0.00357818510383,0.000811365433037] 00177 ) 00178 00179 timing() 00180 00181 00182 # Tsys table generation 00183 mystep = 2 00184 if(mystep in thesteps): 00185 print 'Step ', mystep, step_title[mystep] 00186 00187 00188 if(casadef.casa_version>='3.4.0'): 00189 00190 for name in basename: 00191 os.system('rm -rf cal-tsys_'+name+'.fdm') 00192 gencal( 00193 vis=name+'.ms', 00194 caltype='tsys', 00195 caltable='cal-tsys_'+name+'.fdm') 00196 else: 00197 print "Generation of Tsys tables omitted to avoid usage of analysisUtils" 00198 print "We use canned tsys tables instead." 00199 00200 if(makeplots): 00201 for spw in ['9','11','13','15']: 00202 for name in basename: 00203 plotcal(caltable='cal-tsys_'+name+'.fdm', xaxis='freq', yaxis='amp', 00204 spw=spw, subplot=721, overplot=False, 00205 iteration='antenna', plotrange=[0, 0, 40, 180], plotsymbol='.', 00206 figfile='cal-tsys_per_spw_'+spw+'_'+name+'.png') 00207 00208 timing() 00209 00210 00211 # correct Titan position 00212 mystep = 3 00213 if(mystep in thesteps): 00214 print 'Step ', mystep, step_title[mystep] 00215 00216 for name in basename: 00217 fixplanets(vis=name+'.ms', field='Titan', fixuvw=True, refant=therefant) 00218 00219 timing() 00220 00221 00222 # Apriori flagging 00223 mystep = 4 00224 if(mystep in thesteps): 00225 print 'Step ', mystep, step_title[mystep] 00226 00227 # Create flagcmd input list 00228 flagcmd = ["mode='manual' antenna='CM01'", 00229 "mode='manual' intent='*POINTING*'", 00230 "mode='manual' intent='*ATMOSPHERE*'", 00231 "mode='shadow'", 00232 "mode='manual' autocorr=True"] 00233 00234 for name in basename: 00235 00236 flagmanager(vis=name+'.ms', mode='restore', versionname='Original') 00237 00238 tflagdata(vis=name + '.ms', mode='list', inpfile=flagcmd, flagbackup=False) 00239 00240 if(makeplots): 00241 # Plot amplitude vs time 00242 plotms(vis=name+'.ms', xaxis='time', yaxis='amp', spw='1', 00243 averagedata=T, avgchannel='4000', coloraxis='field', 00244 iteraxis='spw') 00245 00246 tb.clearlocks() 00247 00248 timing() 00249 00250 00251 # WVR cal table generation 00252 mystep = 5 00253 if(mystep in thesteps): 00254 print 'Step ', mystep, step_title[mystep] 00255 00256 for name in basename: 00257 os.system('rm -rf cal-wvr_'+name) 00258 00259 # os.system('/home/arcproc/opsw/alma/wvrgcal-stable/bin/wvrgcal --ms X54.ms --output cal-wvr_X54 --toffset -1 --wvrflag CM01 --segsource --tie "1224+213 Phase","M100" --statsource "1224+213 Phase"') 00260 00261 # os.system('/home/arcproc/opsw/alma/wvrgcal-stable/bin/wvrgcal --ms X220.ms --output cal-wvr_X220 --toffset -1 --wvrflag CM01 --segsource --tie "1224+213 Phase","M100" --statsource "1224+213 Phase"') 00262 wvrgcal(vis="X54.ms", caltable='cal-wvr_X54', toffset=-1, wvrflag=['CM01'], segsource=True, tie=["1224+213 Phase,M100"], statsource="1224+213 Phase") 00263 wvrgcal(vis='X220.ms', caltable='cal-wvr_X220', toffset=-1, wvrflag=['CM01'], segsource=True, tie=["1224+213 Phase,M100"], statsource="1224+213 Phase") 00264 00265 timing() 00266 00267 00268 # delay calibration table generation 00269 mystep = 6 00270 if(mystep in thesteps): 00271 print 'Step ', mystep, step_title[mystep] 00272 00273 for name in basename: 00274 os.system('rm -rf cal-delay_'+name+'.K') 00275 gaincal(vis=name+'.ms',caltable='cal-delay_'+name+'.K', 00276 field='*Phase*',spw='1,3,5,7', 00277 solint='inf',combine='scan',refant=therefant, 00278 gaintable='cal-antpos_'+name, 00279 gaintype='K') 00280 00281 timing() 00282 00283 00284 # applycal 00285 mystep = 7 00286 if(mystep in thesteps): 00287 print 'Step ', mystep, step_title[mystep] 00288 00289 tsysinterp = 'nearest' 00290 tsysspwmap=[] 00291 if(casadef.casa_version>='3.4.0'): 00292 print "Using linear interpolation for Tsys in applycal ..." 00293 tsysinterp='linear' 00294 tsysspwmap=[0,9,0,11,0,13,0,15] 00295 00296 for myfield in ['3c273 - Bandpass','Titan','3c273 - Phase','1224+213 Phase','M100']: 00297 print "Field ", myfield 00298 for name in basename: 00299 applycal(vis=name+'.ms', 00300 spw='1,3,5,7', 00301 field=myfield, gainfield=[myfield,myfield,'',''], 00302 interp=['nearest',tsysinterp,'nearest','nearest'], 00303 spwmap=[[],tsysspwmap,[],[]], 00304 gaintable=['cal-wvr_'+name,'cal-tsys_'+name+'.fdm','cal-delay_'+name+'.K','cal-antpos_'+name], 00305 flagbackup=F) 00306 00307 timing() 00308 00309 00310 # Split and save flags 00311 mystep = 8 00312 if(mystep in thesteps): 00313 print 'Step ', mystep, step_title[mystep] 00314 00315 for name in basename: 00316 os.system('rm -rf '+name+'-line.ms') 00317 split(vis=name+'.ms', 00318 outputvis=name+'-line.ms', 00319 spw='1,3,5,7', 00320 datacolumn='corrected') 00321 flagmanager(vis = name+'-line.ms', 00322 mode = 'save', 00323 versionname = 'apriori') 00324 timing() 00325 00326 00327 # flagging 00328 mystep = 9 00329 if(mystep in thesteps): 00330 print 'Step ', mystep, step_title[mystep] 00331 00332 # Create flagcmd input list (could also call tflagdata twice alternatively) 00333 flagcmd = ["mode='manual' field='' spw='0~3:0~10;3800~3839'", 00334 "mode='manual' field='' spw='0~3:239;447~448;720~721;2847~2848'"] 00335 00336 for name in basename: 00337 00338 flagmanager(vis=name + '-line.ms', mode='restore', 00339 versionname='apriori') 00340 00341 tflagdata(vis=name + '-line.ms', mode='list', inpfile=flagcmd, flagbackup=False) 00342 00343 # some integrations are off 00344 tflagdata(vis='X220-line.ms', mode='manual', 00345 timerange='19:52:55~19:53:04', flagbackup=False) 00346 00347 tflagdata(vis='X54-line.ms', 00348 antenna='PM01', 00349 timerange='19:03:35~19:03:42', 00350 mode='manual', 00351 flagbackup=False) 00352 00353 tflagdata(vis='X54-line.ms', 00354 antenna='DV04', 00355 timerange='19:38:45~19:38:55', 00356 mode='manual', 00357 flagbackup=False) 00358 00359 timing() 00360 00361 # Bin it up to lower spectral resolution to about 10 km/s 00362 mystep = 10 00363 if(mystep in thesteps): 00364 print 'Step ', mystep, step_title[mystep] 00365 00366 for name in basename: 00367 os.system('rm -rf '+name+'-line-vs.ms') 00368 split(vis=name+'-line.ms', 00369 outputvis=name+'-line-vs.ms', 00370 datacolumn='data', 00371 width='8') 00372 00373 flagmanager(vis=name+'-line-vs.ms', mode='save', versionname='apriori') 00374 00375 timing() 00376 00377 00378 # Fast phase-only gaincal for bandpass 00379 mystep = 11 00380 if(mystep in thesteps): 00381 print 'Step ', mystep, step_title[mystep] 00382 00383 for name in basename: 00384 os.system('rm -rf cal-'+name+'-BPint.Gp') 00385 gaincal(vis=name+'-line-vs.ms', 00386 caltable='cal-'+name+'-BPint.Gp', 00387 spw='*:190~310', 00388 field='*Bandpass*', 00389 selectdata=F, solint='int', refant=therefant, calmode='p') 00390 00391 if(makeplots): 00392 plotcal(caltable='cal-'+name+'-BPint.Gp', 00393 xaxis = 'time', yaxis = 'phase', 00394 poln='X', plotsymbol='o', plotrange = [0,0,-180,180], 00395 iteration = 'spw', 00396 figfile='cal-'+name+'-phase_vs_time_XX.BPint.Gp.png', subplot = 221) 00397 00398 plotcal(caltable='cal-'+name+'-BPint.Gp', 00399 xaxis = 'time', yaxis = 'phase', 00400 poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], 00401 iteration = 'spw', 00402 figfile='cal-'+name+'-phase_vs_time_YY.BPint.Gp.png', subplot = 221) 00403 00404 timing() 00405 00406 # Bandpass 00407 mystep = 12 00408 if(mystep in thesteps): 00409 print 'Step ', mystep, step_title[mystep] 00410 00411 for name in basename: 00412 os.system('rm -rf cal-'+name+'.B1') 00413 00414 bandpass(vis=name+'-line-vs.ms', 00415 caltable='cal-'+name+'.B1', 00416 field='*Bandpass*', 00417 bandtype='B', fillgaps=10, solnorm = T, combine='', 00418 selectdata=F, 00419 solint='inf', 00420 refant=therefant, 00421 gaintable='cal-'+name+'-BPint.Gp') 00422 00423 if(makeplots): 00424 for spw in ['0','1','2','3']: 00425 plotcal(caltable = 'cal-'+name+'.B1', 00426 xaxis='freq', yaxis='phase', spw=spw, antenna='', 00427 iteration='antenna', 00428 subplot=431, overplot=False, plotrange = [0,0,-70,70], 00429 plotsymbol='.', timerange='', 00430 figfile='cal-'+name+'-phase.spw'+spw+'.B1.png') 00431 00432 plotcal(caltable = 'cal-'+name+'.B1', 00433 xaxis='freq', yaxis='amp', spw=spw, 00434 iteration='antenna', 00435 subplot=431, overplot=False, 00436 plotsymbol='.', timerange='', 00437 figfile='cal-'+name+'-amplitude.spw'+spw+'.B1.png') 00438 timing() 00439 00440 00441 # Setjy 00442 # Strong line for Titan is obvious 00443 # Noisy for uvdistances less than 40 klambda 00444 mystep = 13 00445 if(mystep in thesteps): 00446 print 'Step ', mystep, step_title[mystep] 00447 00448 for name in basename: 00449 flagmanager(vis=name+'-line-vs.ms', mode='restore', versionname='apriori') 00450 00451 tflagdata(vis=name + '-line-vs.ms', 00452 field='Titan', 00453 mode='manual', 00454 uvrange='0~40', 00455 spw='', 00456 flagbackup=False) 00457 00458 tflagdata(vis=name + '-line-vs.ms', 00459 field='Titan', 00460 mode='manual', 00461 uvrange='', 00462 spw='0:200~479', 00463 flagbackup=False) 00464 00465 setjy(vis=name+'-line-vs.ms', 00466 field='Titan', 00467 standard='Butler-JPL-Horizons 2010', 00468 scalebychan=False, 00469 spw='0,1,2,3') 00470 00471 timing() 00472 00473 # Fast phase-only gaincal 00474 mystep = 14 00475 if(mystep in thesteps): 00476 print 'Step ', mystep, step_title[mystep] 00477 00478 for name in basename: 00479 os.system('rm -rf cal-'+name+'-int.Gp') 00480 gaincal(vis=name+'-line-vs.ms', 00481 caltable='cal-'+name+'-int.Gp', 00482 spw='*:25~455', 00483 field='*Phase*,*Band*,Titan', 00484 gaintable='cal-'+name+'.B1', 00485 selectdata=F, solint='int', 00486 refant=therefant, calmode='p') 00487 00488 if(makeplots): 00489 plotcal(caltable='cal-'+name+'-int.Gp', 00490 xaxis = 'time', yaxis = 'phase', 00491 poln='X', plotsymbol='o', plotrange = [0,0,-180,180], 00492 iteration = 'spw', 00493 figfile='cal-'+name+'-phase_vs_time_XX.int.Gp.png', subplot = 221) 00494 00495 plotcal(caltable='cal-'+name+'-int.Gp', 00496 xaxis = 'time', yaxis = 'phase', 00497 poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], 00498 iteration = 'spw', 00499 figfile='cal-'+name+'-phase_vs_time_YY.int.Gp.png', subplot = 221) 00500 00501 timing() 00502 00503 00504 # Slow phase-only gaincal 00505 mystep = 15 00506 if(mystep in thesteps): 00507 print 'Step ', mystep, step_title[mystep] 00508 00509 for name in basename: 00510 os.system('rm -rf cal-'+name+'-scan.Gp') 00511 gaincal(vis=name+'-line-vs.ms', 00512 caltable='cal-'+name+'-scan.Gp', 00513 spw='*:25~455', 00514 field='*Phase*,*Band*,Titan', 00515 gaintable='cal-'+name+'.B1', 00516 selectdata=F, solint='inf', 00517 refant=therefant, calmode='p') 00518 00519 if(makeplots): 00520 plotcal(caltable='cal-'+name+'-scan.Gp', 00521 xaxis = 'time', yaxis = 'phase', 00522 poln='X', plotsymbol='o', plotrange = [0,0,-180,180], 00523 iteration = 'spw', 00524 figfile='cal-'+name+'-phase_vs_time_XX.scan.Gp.png', subplot = 221) 00525 00526 plotcal(caltable='cal-'+name+'-scan.Gp', 00527 xaxis = 'time', yaxis = 'phase', 00528 poln='Y', plotsymbol='o', plotrange = [0,0,-180,180], 00529 iteration = 'spw', 00530 figfile='cal-'+name+'-phase_vs_time_YY.scan.Gp.png', subplot = 221) 00531 00532 timing() 00533 00534 00535 # Slow amplitude and phase gaincal 00536 mystep = 16 00537 if(mystep in thesteps): 00538 print 'Step ', mystep, step_title[mystep] 00539 00540 for name in basename: 00541 os.system('rm -rf cal-'+name+'-scan.Gap') 00542 gaincal(vis=name+'-line-vs.ms', 00543 caltable='cal-'+name+'-scan.Gap', 00544 spw='*:25~455', 00545 field='*Phase*,*Band*,Titan', 00546 gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp'], 00547 selectdata=F, solint='inf', 00548 refant=therefant, calmode='ap') 00549 00550 if(makeplots): 00551 plotcal(caltable='cal-'+name+'-scan.Gap', 00552 xaxis = 'time', yaxis = 'amp', 00553 poln='X', plotsymbol='o', 00554 iteration = 'spw', 00555 figfile='cal-'+name+'-amp_vs_time_XX.scan.Gap.png', subplot = 221) 00556 00557 plotcal(caltable='cal-'+name+'-scan.Gap', 00558 xaxis = 'time', yaxis = 'amp', 00559 poln='Y', plotsymbol='o', 00560 iteration = 'spw', 00561 figfile='cal-'+name+'-amp_vs_time_YY.scan.Gap.png', subplot = 221) 00562 00563 timing() 00564 00565 # Fluxscale 00566 mystep = 17 00567 if(mystep in thesteps): 00568 print 'Step ', mystep, step_title[mystep] 00569 00570 for name in basename: 00571 fluxscale(vis=name+'-line-vs.ms',caltable='cal-'+name+'-scan.Gap', 00572 fluxtable='cal-'+name+'.flux',reference='Titan',transfer='*Phase*,*Band*') 00573 00574 # Restore flags on calibrators? 00575 #for name in basename: 00576 # flagmanager(vis=name+'-line-vs.ms',mode='restore',versionname='apriori') 00577 00578 timing() 00579 00580 00581 # Applycal 00582 mystep = 18 00583 if(mystep in thesteps): 00584 print 'Step ', mystep, step_title[mystep] 00585 00586 for name in basename: 00587 # to the bandpass cal 00588 applycal(vis=name+'-line-vs.ms',field='*Band*', 00589 gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'], 00590 interp=['nearest','nearest','nearest'], 00591 gainfield=['*Band*','*Band*','*Band*'], 00592 calwt=F, 00593 flagbackup=T) 00594 00595 # to the secondary phase cal 00596 applycal(vis=name+'-line-vs.ms',field='3c273 - Phase', 00597 gaintable=['cal-'+name+'.B1','cal-'+name+'-scan.Gp','cal-'+name+'.flux'], 00598 interp=['nearest','linear','linear'], 00599 gainfield=['*Band*','1224*','1224*'], 00600 calwt=F, 00601 flagbackup=T) 00602 00603 # to the primary phase cal 00604 applycal(vis=name+'-line-vs.ms',field='1224*', 00605 gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'], 00606 interp=['nearest','nearest','nearest'], 00607 gainfield=['*Band*','1224*','1224*'], 00608 calwt=F, 00609 flagbackup=T) 00610 00611 # to Titan 00612 applycal(vis=name+'-line-vs.ms',field='Titan', 00613 gaintable=['cal-'+name+'.B1','cal-'+name+'-int.Gp','cal-'+name+'.flux'], 00614 interp=['nearest','nearest','nearest'], 00615 gainfield=['*Band*','Titan','Titan'], 00616 calwt=F, 00617 flagbackup=T) 00618 00619 # to M100 00620 applycal(vis=name+'-line-vs.ms',field='M100', 00621 gaintable=['cal-'+name+'.B1','cal-'+name+'-scan.Gp','cal-'+name+'.flux'], 00622 interp=['nearest','linear','linear'], 00623 gainfield=['*Band*','1224*','1224*'], 00624 calwt=F, 00625 flagbackup=T) 00626 00627 00628 # For X146 the calibrated fluxes for the secondary phase cal are different for the different pols. What 00629 # is going on? 00630 00631 timing() 00632 00633 00634 # Test image of the secondary phase cal 00635 mystep = 19 00636 if(mystep in thesteps): 00637 print 'Step ', mystep, step_title[mystep] 00638 00639 for name in basename: 00640 os.system('rm -rf test-'+name+'-sec_phasecal*') 00641 clean(vis=name+'-line-vs.ms', 00642 imagename='test-'+name+'-sec_phasecal', 00643 field='3c*Ph*',spw='0~3', 00644 nterms=2, 00645 mode='mfs',niter=100, 00646 interactive=False, 00647 mask=[96, 96, 104, 104],imsize=200,cell='0.5arcsec') 00648 if makeplots: 00649 for name in basename: 00650 imview(raster={'file': 'test-'+name+'-sec_phasecal.image.tt0', 'colorwedge':T, 00651 'range':[-0.02, 8.0], 'scaling':-1.5, 'colormap':'Rainbow 2'}, 00652 out='test-'+name+'-sec_phasecal.png', zoom=1) 00653 00654 timing() 00655 00656 00657 # Test image of the primary phase cal 00658 mystep = 20 00659 if(mystep in thesteps): 00660 print 'Step ', mystep, step_title[mystep] 00661 00662 for name in basename: 00663 os.system('rm -rf test-'+name+'-prim_phasecal*') 00664 clean(vis=name+'-line-vs.ms', 00665 imagename='test-'+name+'-prim_phasecal', 00666 field='1224*',spw='0~3', 00667 nterms=2, 00668 mode='mfs',niter=100, 00669 interactive=False, 00670 mask=[96, 96, 104, 104],imsize=200,cell='0.5arcsec') 00671 00672 if makeplots: 00673 for name in basename: 00674 imview(raster={'file': 'test-'+name+'-prim_phasecal.image.tt0', 'colorwedge':T, 00675 'range':[-0.005, 0.9], 'scaling':-2.5, 'colormap':'Rainbow 2'}, 00676 out='test-'+name+'-prim_phasecal.png', zoom=1) 00677 calstat=imstat(imagename='test-'+name+'-prim_phasecal.image.tt0', region='', box='30,30,170,80') 00678 rms=(calstat['rms'][0]) 00679 print '>> rms in phase calibrator image: '+str(rms) 00680 calstat=imstat(imagename='test-'+name+'-prim_phasecal.image.tt0', region='') 00681 peak=(calstat['max'][0]) 00682 print '>> Peak in phase calibrator image: '+str(peak) 00683 print '>> Dynamic range in phase calibrator image: '+str(peak/rms) 00684 00685 timing() 00686 00687 # Test image of Titan 00688 mystep = 21 00689 if(mystep in thesteps): 00690 print 'Step ', mystep, step_title[mystep] 00691 00692 for name in basename: 00693 os.system('rm -rf test-'+name+'-Titan*') 00694 clean(vis=name+'-line-vs.ms', 00695 imagename='test-'+name+'-Titan', 00696 field='Titan',spw='0~3', 00697 mode='mfs',niter=100, 00698 interactive=False, 00699 mask=[96, 96, 104, 104],imsize=200,cell='0.5arcsec') 00700 00701 timing() 00702 00703 00704 # Split off the calibrated data on M100 00705 mystep = 22 00706 if(mystep in thesteps): 00707 print 'Step ', mystep, step_title[mystep] 00708 00709 for name in basename: 00710 os.system('rm -rf '+name+'-calibrated.ms*') 00711 split(vis=name+'-line-vs.ms',field='M100', 00712 outputvis=name+'-calibrated.ms', 00713 datacolumn = 'corrected', 00714 keepflags=F) 00715 00716 timing() 00717 00718 # Concat 00719 mystep = 23 00720 if(mystep in thesteps): 00721 print 'Step ', mystep, step_title[mystep] 00722 00723 os.system('rm -rf M100all.ms*') 00724 concat(vis=['X54-calibrated.ms', 'X220-calibrated.ms'], 00725 concatvis='M100all.ms', copypointing=False) 00726 00727 timing() 00728 00729 00730 # rebin the data 00731 mystep = 24 00732 if(mystep in thesteps): 00733 print 'Step ', mystep, step_title[mystep] 00734 00735 os.system('rm -rf M100all_lores.ms*') 00736 split(vis='M100all.ms', outputvis='M100all_lores.ms', 00737 datacolumn='data', 00738 timebin='60s') 00739 00740 timing() 00741 00742 00743 # Continuum image 00744 mystep = 25 00745 if(mystep in thesteps): 00746 print 'Step ', mystep, step_title[mystep] 00747 00748 os.system('rm -rf M100cont.*') 00749 clean(vis = 'M100all_lores.ms', 00750 imagename = 'M100cont', 00751 field='2~47', 00752 spw='0:10~210;256~440,1~3:10~460', 00753 mode = 'mfs', 00754 niter = 1000, 00755 mask='M100cont-orig.mask', 00756 imagermode = 'mosaic', 00757 interactive = F, 00758 imsize = 200, 00759 cell = '0.5arcsec', 00760 phasecenter='J2000 12h22m54.9 +15d49m15') 00761 00762 # Continuum peak is 0.5 mJy. Too weak for self-cal... 00763 00764 timing() 00765 00766 00767 # uvcontsub2 00768 mystep = 26 00769 if(mystep in thesteps): 00770 print 'Step ', mystep, step_title[mystep] 00771 00772 os.system('rm -rf M100all_lores.ms.c*') 00773 uvcontsub2(vis='M100all_lores.ms',field='',fitspw='0:10~205;260~440', 00774 combine='',solint='inf',fitorder=1,spw='0',want_cont=False) 00775 00776 timing() 00777 00778 00779 # Test image of central field 00780 mystep = 27 00781 if(mystep in thesteps): 00782 print 'Step ', mystep, step_title[mystep] 00783 00784 os.system('rm -rf test-M100line.*') 00785 clean( 00786 vis='M100all_lores.ms.contsub', 00787 imagename='test-M100line', 00788 field='26', 00789 spw='0:231~248', 00790 mode='mfs', 00791 niter=500,gain=0.1,threshold='0.0mJy', 00792 imagermode='csclean', 00793 interactive=False, 00794 mask='test-M100line-orig.mask', 00795 outframe='',veltype='radio', 00796 imsize=200,cell='0.5arcsec', 00797 phasecenter='', 00798 stokes='I', 00799 weighting='briggs',robust=0.5, 00800 npercycle=100,cyclefactor=1.5,cyclespeedup=-1) 00801 00802 timing() 00803 00804 00805 # Clean line cube mosaic 00806 mystep = 28 00807 if(mystep in thesteps): 00808 print 'Step ', mystep, step_title[mystep] 00809 00810 os.system('rm -rf M100line.*') 00811 clean(vis='M100all_lores.ms.contsub',imagename='M100line', 00812 field='2~47', 00813 spw='0:220~259', 00814 mode='channel', 00815 niter=1000,gain=0.1,threshold='0.0mJy',psfmode='clark', 00816 imagermode='mosaic',ftmachine='mosaic',mosweight=False, 00817 scaletype='SAULT', 00818 interactive=False, 00819 mask='M100line-orig.mask', 00820 nchan=40,start=220, 00821 width=1, 00822 outframe='',veltype='radio', 00823 imsize=600,cell='0.5arcsec', 00824 phasecenter='J2000 12h22m54.9 +15d49m10', 00825 restfreq='115.271201800GHz',stokes='I', 00826 weighting='briggs',robust=0.5, 00827 pbcor=False,minpb=0.2, 00828 npercycle=100,cyclefactor=1.5,cyclespeedup=-1) 00829 00830 timing() 00831 00832 # Moments 00833 mystep = 29 00834 if(mystep in thesteps): 00835 print 'Step ', mystep, step_title[mystep] 00836 00837 os.system('rm -rf M100-CO.mom?') 00838 immoments(imagename='M100line.image', 00839 moments=[0], 00840 axis='spectral', 00841 region='',box='100,110,515,500', 00842 chans='7~35', 00843 mask='', 00844 outfile='M100-CO.mom0', 00845 includepix=[0.03, 1000000]) 00846 00847 immoments(imagename='M100line.image', 00848 moments=[1], 00849 axis='spectral', 00850 region='',box='100,110,515,500', 00851 chans='7~35', 00852 mask='', 00853 outfile='M100-CO.mom1', 00854 includepix=[0.035, 1000000]) 00855 00856 if makeplots: 00857 os.system('rm -rf M100-CO_velfield.png') 00858 imview(contour={'file': 'M100-CO.mom0','levels': 00859 [1,2,5,10,20,40,80,160],'base':0,'unit':1}, 00860 raster={'file': 'M100-CO.mom1','range': [1440,1700], 00861 'colorwedge':T, 'colormap': 'Rainbow 2'}, out='M100-CO_velfield.png') 00862 os.system('rm -rf M100-CO_map.png') 00863 imview(contour={'file': 'M100-CO.mom1','levels': 00864 [1430,1460,1490,1520,1550,1580,1610,1640,1670,1700],'base':0,'unit':1}, 00865 raster={'file': 'M100-CO.mom0', 'colorwedge':T, 00866 'colormap': 'Rainbow 2','scaling':-1.8,'range': [0.5,20]}, 00867 out='M100-CO_map.png') 00868 00869 os.system('rm -rf M100-CO_contmap.png') 00870 imview(contour={'file': 'M100cont.image','levels': 00871 [0.00025,0.0004],'base':0,'unit':1}, 00872 zoom=3, 00873 raster={'file': 'M100-CO.mom0', 'colorwedge':T, 00874 'colormap': 'Rainbow 2','scaling':0,'range': [0.8,40]}, 00875 out='M100-CO_contmap.png') 00876 00877 timing() 00878 00879 mystep = 30 00880 if(mystep in thesteps): 00881 print 'Step ', mystep, step_title[mystep] 00882 00883 resrms = [] 00884 respeak = [] 00885 00886 # expectation values set 14 March 2012 based on analysis using CASA 3.3 00887 exppeak33 = [1.11061167717, 1.08436012268] 00888 exprms33 = [0.000449335755548, 0.000499602989294] 00889 #expectation values set 15 March 2012 based on analysis using CASA active r18746 00890 exppeakr18746 = [1.11075341702, 1.08440160751] 00891 exprmsr18746 = [0.000527012278326, 0.000579607207328] # note worse RMS 00892 #expectation values set 17 Sept 2012 based on analysis using CASA active r21038 00893 # (change was due to a modification of the fluxscale task in r21037) 00894 exppeak = [1.11501789093, 1.08849704266] 00895 exprms = [0.000528678297997, 0.000582209031563] 00896 00897 00898 for name in basename: 00899 00900 calstat=imstat(imagename='test-'+name+'-prim_phasecal.image.tt0', region='', box='30,30,170,80') 00901 rms=(calstat['rms'][0]) 00902 casalog.post( name+': rms in phase calibrator image: '+str(rms)) 00903 calstat=imstat(imagename='test-'+name+'-prim_phasecal.image.tt0', region='') 00904 peak=(calstat['max'][0]) 00905 casalog.post( name+': Peak in phase calibrator image: '+str(peak)) 00906 casalog.post( name+': Dynamic range in phase calibrator image: '+str(peak/rms)) 00907 00908 resrms.append(rms) 00909 respeak.append(peak) 00910 00911 resrmsm = [] 00912 respeakm = [] 00913 00914 # expectation values set 14 March 2012 based on analysis using CASA 3.3 00915 exppeakm33 = 0.164112448692 00916 exprmsm33 = 0.0083269206807 00917 #expectation values set 15 March 2012 based on analysis using CASA active r18746 00918 exppeakm = 0.163885176182 00919 exprmsm = 0.00828791689128 00920 #expectation values set 17 Sept 2012 based on analysis using CASA active r21038 00921 exppeakm = 0.164736732841 00922 exprmsm = 0.00830688327551 00923 00924 calstat=imstat(imagename='test-M100line.image', region='', box='42,115,65,134') 00925 resrmsm=(calstat['rms'][0]) 00926 calstat=imstat(imagename='test-M100line.image', region='') 00927 respeakm=(calstat['max'][0]) 00928 00929 casalog.post( ' rms in M100: '+str(resrmsm)) 00930 casalog.post( ' Peak in M100: '+str(respeakm)) 00931 casalog.post( ' Dynamic range in M100: '+str(respeakm/resrmsm)) 00932 00933 00934 timing() 00935 00936 # print results to logger 00937 casalog.origin('SUMMARY') 00938 casalog.post("\n***** Peak and RMS of the images of the primary phase calibrator *****") 00939 casalog.post( "Dataset, Peak (expectation, expectation CASA 3.3), RMS (expectation, expectation CASA 3.3)") 00940 casalog.post( "------------------------------------------------------------------------------------------") 00941 for i in [0,1]: 00942 casalog.post( basename[i]+","+ str(respeak[i])+ "("+str(exppeak[i])+","+ str(exppeak33[i])+")," 00943 + str(resrms[i])+ "("+str(exprms[i])+","+str(exprms33[i])+")") 00944 00945 casalog.post( "------------------------------------------------------------------------------------------") 00946 00947 casalog.post( "\n***** Peak and RMS of the image of the central field of the M100 mosaic *****") 00948 00949 casalog.post( "M100: Peak (expectation, expectation CASA 3.3), RMS (expectation, expectation CASA 3.3)") 00950 casalog.post( "------------------------------------------------------------------------------------------") 00951 casalog.post( str(respeakm)+ "("+str(exppeakm)+","+ str(exppeakm33)+"),"+ str(resrmsm)+ "("+str(exprmsm)+","+str(exprmsm33)+")") 00952 casalog.post( "------------------------------------------------------------------------------------------") 00953 00954 00955 passed = True 00956 00957 for i in [0,1]: 00958 peakdev = abs(respeak[i]-exppeak[i])/exppeak[i]*100. 00959 if (peakdev > 0.5): 00960 casalog.post( 'ERROR: Peak in primary phase calibrator image '+str(i)+' deviates from expectation by '+str(peakdev)+' percent.', 'WARN') 00961 passed = False 00962 00963 rmsdev = abs(resrms[i]-exprms[i])/exprms[i]*100. 00964 if (rmsdev > 0.5): 00965 casalog.post( 'ERROR: RMS in primary phase calibrator image '+str(i)+' deviates from expectation by '+str(rmsdev)+' percent.','WARN') 00966 passed = False 00967 00968 peakmdev = abs(respeakm-exppeakm)/exppeakm*100. 00969 if (peakmdev > 0.5): 00970 casalog.post( 'ERROR: Peak in M100 central field image '+str(i)+' deviates from expectation by '+str(peakmdev)+' percent.','WARN') 00971 passed = False 00972 00973 rmsmdev = abs(resrmsm-exprmsm)/exprmsm*100. 00974 if (rmsmdev > 0.5): 00975 casalog.post( 'ERROR: RMS in M100 central field image '+str(i)+' deviates from expectation by '+str(rmsmdev)+' percent.','WARN') 00976 passed = False 00977 00978 if not os.path.exists('M100-CO.mom0'): 00979 casalog.post( 'ERROR: M100 line cube moment 0 map was not created!','WARN') 00980 passed = False 00981 00982 if not os.path.exists('M100-CO.mom1'): 00983 casalog.post( 'ERROR: M100 line cube moment 1 map was not created!','WARN') 00984 passed = False 00985 00986 if not passed: 00987 raise Exception, 'Results are different from expectations by more than 0.5 percent.' 00988 00989 casalog.post( "\nAll peak and RMS values within 0.5 percent of the expectation.") 00990 00991 00992 print 'Done.'
1.8.0
