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