- Timestamp:
- 07/20/11 15:20:39 (13 years ago)
- Location:
- trunk/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/MSFiller.cpp
r2234 r2237 95 95 { 96 96 os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ; 97 //double startSec = gettimeofday_sec() ;98 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;97 //double startSec = gettimeofday_sec() ; 98 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ; 99 99 //os_ << " filename = " << filename << endl ; 100 100 … … 147 147 isData_ = mstable_.tableDesc().isColumn( "DATA" ) ; 148 148 149 //double endSec = gettimeofday_sec() ;150 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;149 //double endSec = gettimeofday_sec() ; 150 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 151 151 return true ; 152 152 } … … 155 155 { 156 156 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ; 157 //double startSec = gettimeofday_sec() ;158 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;159 160 //double time0 = gettimeofday_sec() ;161 //os_ << "start init fill: " << time0 << LogIO::POST ;157 //double startSec = gettimeofday_sec() ; 158 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ; 159 160 //double time0 = gettimeofday_sec() ; 161 //os_ << "start init fill: " << time0 << LogIO::POST ; 162 162 163 163 // Initialize header 164 164 STHeader sdh ; 165 sdh.nchan = 0 ; 166 sdh.npol = 0 ; 167 sdh.nif = 0 ; 168 sdh.nbeam = 0 ; 169 sdh.observer = "" ; 170 sdh.project = "" ; 171 sdh.obstype = "" ; 172 sdh.antennaname = "" ; 173 sdh.antennaposition.resize( 0 ) ; 174 sdh.equinox = 0.0 ; 175 sdh.freqref = "" ; 176 sdh.reffreq = -1.0 ; 177 sdh.bandwidth = 0.0 ; 178 sdh.utc = 0.0 ; 179 sdh.fluxunit = "" ; 180 sdh.epoch = "" ; 181 sdh.poltype = "" ; 165 initHeader( sdh ) ; 182 166 183 167 // check if optional table exists … … 269 253 //fillHistory() ; 270 254 271 // shared pointers272 ROTableColumn *tcolr ;273 274 255 // MAIN 275 256 // Iterate over several ids … … 280 261 MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ; 281 262 if ( getPt_ ) { 282 //pointtab = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ;283 263 pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ; 284 264 } 285 tcolr = tpoolr->construct( anttab, "STATION" ) ; 286 String stationName = tcolr->asString( antenna_ ) ; 287 tpoolr->destroy( tcolr ) ; 288 tcolr = tpoolr->construct( anttab, "NAME" ) ; 289 String antennaName = tcolr->asString( antenna_ ) ; 290 tpoolr->destroy( tcolr ) ; 265 String stationName = asString( "STATION", antenna_, anttab, tpoolr ) ; 266 String antennaName = asString( "NAME", antenna_, anttab, tpoolr ) ; 291 267 sdh.antennaposition.resize( 3 ) ; 292 268 for ( int i = 0 ; i < 3 ; i++ ) … … 294 270 String telescopeName = "" ; 295 271 296 //double time1 = gettimeofday_sec() ;297 //os_ << "end fill init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;272 //double time1 = gettimeofday_sec() ; 273 //os_ << "end fill init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 298 274 299 275 // row based … … 334 310 TableIterator iter0( mstable_, "OBSERVATION_ID" ) ; 335 311 while( !iter0.pastEnd() ) { 336 //time0 = gettimeofday_sec() ;337 //os_ << "start 0th iteration: " << time0 << LogIO::POST ;312 //time0 = gettimeofday_sec() ; 313 //os_ << "start 0th iteration: " << time0 << LogIO::POST ; 338 314 Table t0 = iter0.table() ; 339 tcolr = tpoolr->construct( t0, "OBSERVATION_ID" ) ; 340 Int obsId = tcolr->asInt( 0 ) ; 341 tpoolr->destroy( tcolr ) ; 315 Int obsId = asInt( "OBSERVATION_ID", 0, t0, tpoolr ) ; 342 316 if ( sdh.observer == "" ) { 343 tcolr = tpoolr->construct( obstab, "OBSERVER" ) ; 344 sdh.observer = tcolr->asString( obsId ) ; 345 tpoolr->destroy( tcolr ) ; 317 sdh.observer = asString( "OBSERVER", obsId, obstab, tpoolr ) ; 346 318 } 347 319 if ( sdh.project == "" ) { 348 tcolr = tpoolr->construct( obstab, "PROJECT" ) ; 349 sdh.observer = tcolr->asString( obsId ) ; 350 tpoolr->destroy( tcolr ) ; 320 sdh.project = asString( "PROJECT", obsId, obstab, tpoolr ) ; 351 321 } 352 322 ROArrayMeasColumn<MEpoch> *tmpMeasCol = new ROArrayMeasColumn<MEpoch>( obstab, "TIME_RANGE" ) ; … … 357 327 } 358 328 if ( telescopeName == "" ) { 359 tcolr = tpoolr->construct( obstab, "TELESCOPE_NAME" ) ; 360 telescopeName = tcolr->asString( obsId ) ; 361 tpoolr->destroy( tcolr ) ; 329 telescopeName = asString( "TELESCOPE_NAME", obsId, obstab, tpoolr ) ; 362 330 } 363 331 Int nbeam = 0 ; 364 //time1 = gettimeofday_sec() ;365 //os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;332 //time1 = gettimeofday_sec() ; 333 //os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 366 334 // 367 335 // ITERATION: FEED1 … … 369 337 TableIterator iter1( t0, "FEED1" ) ; 370 338 while( !iter1.pastEnd() ) { 371 //time0 = gettimeofday_sec() ;372 //os_ << "start 1st iteration: " << time0 << LogIO::POST ;339 //time0 = gettimeofday_sec() ; 340 //os_ << "start 1st iteration: " << time0 << LogIO::POST ; 373 341 Table t1 = iter1.table() ; 374 342 // assume FEED1 == FEED2 375 tcolr = tpoolr->construct( t1, "FEED1" ) ; 376 Int feedId = tcolr->asInt( 0 ) ; 377 tpoolr->destroy( tcolr ) ; 343 Int feedId = asInt( "FEED1", 0, t1, tpoolr ) ; 378 344 nbeam++ ; 379 345 … … 386 352 *uintRF = 0 ; 387 353 388 //time1 = gettimeofday_sec() ;389 //os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;354 //time1 = gettimeofday_sec() ; 355 //os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 390 356 // 391 357 // ITERATION: FIELD_ID … … 393 359 TableIterator iter2( t1, "FIELD_ID" ) ; 394 360 while( !iter2.pastEnd() ) { 395 //time0 = gettimeofday_sec() ;396 //os_ << "start 2nd iteration: " << time0 << LogIO::POST ;361 //time0 = gettimeofday_sec() ; 362 //os_ << "start 2nd iteration: " << time0 << LogIO::POST ; 397 363 Table t2 = iter2.table() ; 398 tcolr = tpoolr->construct( t2, "FIELD_ID" ) ; 399 Int fieldId = tcolr->asInt( 0 ) ; 400 tpoolr->destroy( tcolr ) ; 401 tcolr = tpoolr->construct( fieldtab, "SOURCE_ID" ) ; 402 Int srcId = tcolr->asInt( fieldId ) ; 403 tpoolr->destroy( tcolr ) ; 404 tcolr = tpoolr->construct( fieldtab, "NAME" ) ; 405 String fieldName = tcolr->asString( fieldId ) + "__" + String::toString(fieldId) ; 406 tpoolr->destroy( tcolr ) ; 364 Int fieldId = asInt( "FIELD_ID", 0, t2, tpoolr ) ; 365 Int srcId = asInt( "SOURCE_ID", fieldId, fieldtab, tpoolr ) ; 366 String fieldName = asString( "NAME", fieldId, fieldtab, tpoolr ) ; 367 fieldName += "__" + String::toString(fieldId) ; 407 368 ROArrayMeasColumn<MDirection> *delayDirCol = new ROArrayMeasColumn<MDirection>( fieldtab, "DELAY_DIR" ) ; 408 369 Vector<MDirection> delayDir = (*delayDirCol)( fieldId ) ; 409 delete delayDirCol ; 410 Vector<Double> defaultScanrate( 2, 0.0 ) ; 411 Vector<Double> defaultDir = delayDir[0].getAngle( "rad" ).getValue() ; 412 if ( delayDir.nelements() > 1 ) 413 defaultScanrate = delayDir[1].getAngle( "rad" ).getValue() ; 414 370 delete delayDirCol ; 415 371 416 372 // FIELDNAME … … 419 375 420 376 421 //time1 = gettimeofday_sec() ;422 //os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;377 //time1 = gettimeofday_sec() ; 378 //os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 423 379 // 424 380 // ITERATION: DATA_DESC_ID … … 426 382 TableIterator iter3( t2, "DATA_DESC_ID" ) ; 427 383 while( !iter3.pastEnd() ) { 428 //time0 = gettimeofday_sec() ;429 //os_ << "start 3rd iteration: " << time0 << LogIO::POST ;384 //time0 = gettimeofday_sec() ; 385 //os_ << "start 3rd iteration: " << time0 << LogIO::POST ; 430 386 Table t3 = iter3.table() ; 431 tcolr = tpoolr->construct( t3, "DATA_DESC_ID" ) ; 432 Int ddId = tcolr->asInt( 0 ) ; 433 tpoolr->destroy( tcolr ) ; 434 tcolr = tpoolr->construct( ddtab, "POLARIZATION_ID" ) ; 435 Int polId = tcolr->asInt( ddId ) ; 436 tpoolr->destroy( tcolr ) ; 437 tcolr = tpoolr->construct( ddtab, "SPECTRAL_WINDOW_ID" ) ; 438 Int spwId = tcolr->asInt( ddId ) ; 439 tpoolr->destroy( tcolr ) ; 387 Int ddId = asInt( "DATA_DESC_ID", 0, t3, tpoolr ) ; 388 Int polId = asInt( "POLARIZATION_ID", ddId, ddtab, tpoolr ) ; 389 Int spwId = asInt( "SPECTRAL_WINDOW_ID", ddId, ddtab, tpoolr ) ; 440 390 441 391 // IFNO … … 444 394 445 395 // polarization information 446 tcolr = tpoolr->construct( poltab, "NUM_CORR" ) ; 447 Int npol = tcolr->asInt( polId ) ; 448 tpoolr->destroy( tcolr ) ; 396 Int npol = asInt( "NUM_CORR", polId, poltab, tpoolr ) ; 449 397 ROArrayColumn<Int> *roArrICol = new ROArrayColumn<Int>( poltab, "CORR_TYPE" ) ; 450 398 Vector<Int> corrtype = (*roArrICol)( polId ) ; 451 399 delete roArrICol ; 452 // os_ << "npol = " << npol << LogIO::POST ;453 // os_ << "corrtype = " << corrtype << LogIO::POST ;454 // source information455 // os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ;456 MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == srcId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ;457 if ( srctabSel.nrow() == 0 ) {458 srctabSel = srctab( srctab.col("SOURCE_ID") == srcId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ;459 }460 400 String srcName( "" ) ; 461 401 Vector<Double> srcPM( 2, 0.0 ) ; 462 402 Vector<Double> srcDir( 2, 0.0 ) ; 463 403 MDirection md ; 464 Int numLines = 0 ; 465 ROArrayColumn<Double> *roArrDCol = 0 ; 466 if ( srctabSel.nrow() > 0 ) { 467 // source name 468 tcolr = tpoolr->construct( srctabSel, "NAME" ) ; 469 srcName = tcolr->asString( 0 ) ; 470 tpoolr->destroy( tcolr ) ; 471 472 // source proper motion 473 roArrDCol = new ROArrayColumn<Double>( srctabSel, "PROPER_MOTION" ) ; 474 srcPM = (*roArrDCol)( 0 ) ; 475 delete roArrDCol ; 476 477 // source direction 478 roArrDCol = new ROArrayColumn<Double>( srctabSel, "DIRECTION" ) ; 479 srcDir = (*roArrDCol)( 0 ) ; 480 delete roArrDCol ; 481 482 // source direction as MDirection object 483 ROScalarMeasColumn<MDirection> *tmpMeasCol = new ROScalarMeasColumn<MDirection>( srctabSel, "DIRECTION" ) ; 484 md = (*tmpMeasCol)( 0 ) ; 485 delete tmpMeasCol ; 486 487 // number of lines 488 tcolr = tpoolr->construct( srctabSel, "NUM_LINES" ) ; 489 numLines = tcolr->asInt( 0 ) ; 490 tpoolr->destroy( tcolr ) ; 491 492 } 493 else { 494 md = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ; 495 } 404 Vector<Double> restFreqs ; 405 Vector<String> transitionName ; 406 Vector<Double> sysVels ; 407 // os_ << "npol = " << npol << LogIO::POST ; 408 // os_ << "corrtype = " << corrtype << LogIO::POST ; 409 410 // source and molecular transition 411 sourceInfo( srcId, spwId, srcName, md, srcPM, restFreqs, transitionName, sysVels, tpoolr ) ; 412 // os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ; 496 413 497 414 // SRCNAME … … 509 426 // SRCDIRECTION 510 427 darrRF.attachToRecord( trec, "SRCDIRECTION" ) ; 511 *darrRF = srcDir ; 428 // *darrRF = srcDir ; 429 *darrRF = md.getAngle().getValue( "rad" ) ; 512 430 513 431 //os_ << "srcDir = " << srcDir << LogIO::POST ; 514 432 515 // for MOLECULES subtable 516 // os_ << "numLines = " << numLines << LogIO::POST ; 517 518 Vector<Double> restFreqs( numLines, 0.0 ) ; 519 Vector<String> transitionName( numLines, "" ) ; 520 Vector<Double> sysVels ; 433 // SRCVELOCITY 521 434 Double sysVel = 0.0 ; 522 if ( numLines != 0 ) { 523 if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) { 524 sharedQDArrCol = new ROArrayQuantColumn<Double>( srctabSel, "REST_FREQUENCY" ) ; 525 Array< Quantum<Double> > qRestFreqs = (*sharedQDArrCol)( 0 ) ; 526 delete sharedQDArrCol ; 527 for ( int i = 0 ; i < numLines ; i++ ) { 528 restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ; 529 } 530 } 531 // os_ << "restFreqs = " << restFreqs << LogIO::POST ; 532 if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) { 533 ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ; 534 if ( transitionCol.isDefined( 0 ) ) 535 transitionName = transitionCol( 0 ) ; 536 //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ; 537 } 538 if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) { 539 roArrDCol = new ROArrayColumn<Double>( srctabSel, "SYSVEL" ) ; 540 sysVels = (*roArrDCol)( 0 ) ; 541 delete roArrDCol ; 542 } 543 if ( !sysVels.empty() ) { 544 //os_ << "sysVels.shape() = " << sysVels.shape() << LogIO::POST ; 545 // NB: assume all SYSVEL values are the same 546 sysVel = sysVels( IPosition(1,0) ) ; 547 } 548 } 549 550 // SRCVELOCITY 435 if ( !sysVels.empty() ) 436 sysVel = sysVels[0] ; 551 437 RecordFieldPtr<Double> doubleRF( trec, "SRCVELOCITY" ) ; 552 438 *doubleRF = sysVel ; … … 562 448 // spectral setup 563 449 uInt freqId ; 564 tcolr = tpoolr->construct( spwtab, "NUM_CHAN" ) ; 565 Int nchan = tcolr->asInt( spwId ) ; 450 Int nchan = asInt( "NUM_CHAN", spwId, spwtab, tpoolr ) ; 566 451 Bool iswvr = False ; 567 452 if ( nchan == 4 ) iswvr = True ; 568 tpoolr->destroy( tcolr) ;453 sdh.nchan = max( sdh.nchan, nchan ) ; 569 454 map<Int,uInt>::iterator iter = ifmap.find( spwId ) ; 570 455 if ( iter == ifmap.end() ) { 571 ROScalarQuantColumn<Double> *tmpQuantCol = new ROScalarQuantColumn<Double>( t3, "TIME" ) ; 572 me = MEpoch( (*tmpQuantCol)( 0 ), MEpoch::UTC ) ; 573 delete tmpQuantCol ; 574 MeasFrame mf( me, mp, md ) ; 575 tcolr = tpoolr->construct( spwtab, "MEAS_FREQ_REF" ) ; 576 MFrequency::Types freqRef = MFrequency::castType((uInt)(tcolr->asInt(spwId))) ; 577 tpoolr->destroy( tcolr ) ; 578 Bool even = False ; 579 if ( (nchan/2)*2 == nchan ) even = True ; 580 sdh.nchan = max( sdh.nchan, nchan ) ; 581 tmpQuantCol = new ROScalarQuantColumn<Double>( spwtab, "TOTAL_BANDWIDTH" ) ; 582 Double totbw = (*tmpQuantCol)( spwId ).getValue( "Hz" ) ; 583 delete tmpQuantCol ; 584 if ( nchan != 4 ) 585 sdh.bandwidth = max( sdh.bandwidth, totbw ) ; 586 if ( sdh.freqref == "" && nchan != 4) 587 //sdh.freqref = MFrequency::showType( freqRef ) ; 588 sdh.freqref = "LSRK" ; 589 if ( sdh.reffreq == -1.0 && nchan != 4 ) { 590 tmpQuantCol = new ROScalarQuantColumn<Double>( spwtab, "REF_FREQUENCY" ) ; 591 Quantum<Double> qreffreq = (*tmpQuantCol)( spwId ) ; 592 delete tmpQuantCol ; 593 // sdh.reffreq = qreffreq.getValue("Hz") ; 594 if ( freqRef == MFrequency::LSRK ) { 595 sdh.reffreq = qreffreq.getValue("Hz") ; 596 } 597 else { 598 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 599 sdh.reffreq = tolsr( qreffreq ).get("Hz").getValue() ; 600 } 601 } 602 Int refchan = nchan / 2 ; 603 IPosition refip( 1, refchan ) ; 604 Double refpix = 0.5*(nchan-1) ; 605 Double refval = 0.0 ; 606 sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_WIDTH" ) ; 607 Double increment = (*sharedQDArrCol)( spwId )( refip ).getValue( "Hz" ) ; 608 delete sharedQDArrCol ; 609 // os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ; 610 sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_FREQ" ) ; 611 Vector< Quantum<Double> > chanFreqs = (*sharedQDArrCol)( spwId ) ; 612 delete sharedQDArrCol ; 613 if ( nchan > 1 && chanFreqs[0].getValue("Hz") > chanFreqs[1].getValue("Hz") ) 614 increment *= -1.0 ; 615 // if ( even ) { 616 // IPosition refip0( 1, refchan-1 ) ; 617 // Double refval0 = chanFreqs(refip0).getValue("Hz") ; 618 // Double refval1 = chanFreqs(refip).getValue("Hz") ; 619 // refval = 0.5 * ( refval0 + refval1 ) ; 620 // } 621 // else { 622 // refval = chanFreqs(refip).getValue("Hz") ; 623 // } 624 if ( freqRef == MFrequency::LSRK ) { 625 if ( even ) { 626 IPosition refip0( 1, refchan-1 ) ; 627 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 628 Double refval1 = chanFreqs(refip).getValue("Hz") ; 629 refval = 0.5 * ( refval0 + refval1 ) ; 630 } 631 else { 632 refval = chanFreqs(refip).getValue("Hz") ; 633 } 634 } 635 else { 636 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 637 if ( even ) { 638 IPosition refip0( 1, refchan-1 ) ; 639 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 640 Double refval1 = chanFreqs(refip).getValue("Hz") ; 641 refval = 0.5 * ( refval0 + refval1 ) ; 642 refval = tolsr( refval ).get("Hz").getValue() ; 643 } 644 else { 645 refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ; 646 } 647 } 456 ROScalarMeasColumn<MEpoch> *tmpMeasCol = new ROScalarMeasColumn<MEpoch>( t3, "TIME" ) ; 457 me = (*tmpMeasCol)( 0 ) ; 458 delete tmpMeasCol ; 459 Double refpix ; 460 Double refval ; 461 Double increment ; 462 spectralSetup( spwId, 463 me, 464 mp, 465 md, 466 refpix, 467 refval, 468 increment, 469 nchan, 470 sdh.freqref, 471 sdh.reffreq, 472 sdh.bandwidth, 473 tpoolr ) ; 648 474 freqId = table_->frequencies().addEntry( refpix, refval, increment ) ; 649 //if ( ifmap.find( spwId ) == ifmap.end() ) {650 475 ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ; 651 476 //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ; … … 680 505 if ( !iswvr && sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ; 681 506 682 //time1 = gettimeofday_sec() ;683 //os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;507 //time1 = gettimeofday_sec() ; 508 //os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 684 509 // 685 510 // ITERATION: SCAN_NUMBER … … 687 512 TableIterator iter4( t3, "SCAN_NUMBER" ) ; 688 513 while( !iter4.pastEnd() ) { 689 //time0 = gettimeofday_sec() ;690 //os_ << "start 4th iteration: " << time0 << LogIO::POST ;514 //time0 = gettimeofday_sec() ; 515 //os_ << "start 4th iteration: " << time0 << LogIO::POST ; 691 516 Table t4 = iter4.table() ; 692 tcolr = tpoolr->construct( t4, "SCAN_NUMBER" ) ; 693 Int scanNum = tcolr->asInt( 0 ) ; 694 tpoolr->destroy( tcolr ) ; 517 Int scanNum = asInt( "SCAN_NUMBER", 0, t4, tpoolr ) ; 695 518 696 519 // SCANNO … … 698 521 *uintRF = scanNum - 1 ; 699 522 700 //time1 = gettimeofday_sec() ;701 //os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;523 //time1 = gettimeofday_sec() ; 524 //os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 702 525 // 703 526 // ITERATION: STATE_ID … … 705 528 TableIterator iter5( t4, "STATE_ID" ) ; 706 529 while( !iter5.pastEnd() ) { 707 //time0 = gettimeofday_sec() ;708 //os_ << "start 5th iteration: " << time0 << LogIO::POST ;530 //time0 = gettimeofday_sec() ; 531 //os_ << "start 5th iteration: " << time0 << LogIO::POST ; 709 532 Table t5 = iter5.table() ; 710 tcolr = tpoolr->construct( t5, "STATE_ID" ) ; 711 Int stateId = tcolr->asInt( 0 ) ; 712 tpoolr->destroy( tcolr ) ; 713 tcolr = tpoolr->construct( stattab, "OBS_MODE" ) ; 714 String obstype = tcolr->asString( stateId ) ; 715 tpoolr->destroy( tcolr ) ; 533 Int stateId = asInt( "STATE_ID", 0, t5, tpoolr ) ; 534 String obstype = asString( "OBS_MODE", 0, stattab, tpoolr ) ; 716 535 if ( sdh.obstype == "" ) sdh.obstype = obstype ; 717 536 718 537 Int nrow = t5.nrow() ; 719 //time1 = gettimeofday_sec() ;720 //os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;538 //time1 = gettimeofday_sec() ; 539 //os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 721 540 722 541 uInt cycle = 0 ; … … 724 543 Cube<Float> spArr ; 725 544 Cube<Bool> flArr ; 726 if ( isFloatData_ ) { 727 ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ; 728 ROArrayColumn<Float> mFloatDataCol( t5, "FLOAT_DATA" ) ; 729 spArr = mFloatDataCol.getColumn() ; 730 flArr = mFlagCol.getColumn() ; 731 if ( sdh.fluxunit == "" ) { 732 const TableRecord &dataColKeys = mFloatDataCol.keywordSet() ; 733 if ( dataColKeys.isDefined( "UNIT" ) ) 734 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ; 735 else if ( dataColKeys.isDefined( "QuantumUnits" ) ) 736 sdh.fluxunit = dataColKeys.asString( "QuantumUnits" ) ; 737 } 545 reshapeSpectraAndFlagtra( spArr, 546 flArr, 547 t5, 548 npol, 549 nchan, 550 nrow, 551 corrtype ) ; 552 if ( sdh.fluxunit == "" ) { 553 String colName = "FLOAT_DATA" ; 554 if ( isData_ ) colName = "DATA" ; 555 ROTableColumn dataCol( t5, colName ) ; 556 const TableRecord &dataColKeys = dataCol.keywordSet() ; 557 if ( dataColKeys.isDefined( "UNIT" ) ) 558 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ; 559 else if ( dataColKeys.isDefined( "QuantumUnits" ) ) 560 sdh.fluxunit = dataColKeys.asString( "QuantumUnits" ) ; 738 561 } 739 else if ( isData_ ) { 740 spArr.resize( npol, nchan, nrow ) ; 741 flArr.resize( npol, nchan, nrow ) ; 742 ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ; 743 ROArrayColumn<Complex> mDataCol( t5, "DATA" ) ; 744 for ( Int irow = 0 ; irow < nrow ; irow++ ) { 745 Bool crossOK = False ; 746 Matrix<Complex> sp = mDataCol( irow ) ; 747 Matrix<Bool> fl = mFlagCol( irow ) ; 748 Matrix<Float> spxy = spArr.xyPlane( irow ) ; 749 Matrix<Bool> flxy = flArr.xyPlane( irow ) ; 750 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 751 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX 752 || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) { 753 if ( !crossOK ) { 754 spxy.row( ipol ) = real( sp.row( ipol ) ) ; 755 flxy.row( ipol ) = fl.row( ipol ) ; 756 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::RL ) { 757 spxy.row( ipol+1 ) = imag( sp.row( ipol ) ) ; 758 flxy.row( ipol+1 ) = fl.row( ipol ) ; 759 } 760 else { 761 spxy.row( ipol+1 ) = imag( conj( sp.row( ipol ) ) ) ; 762 flxy.row( ipol+1 ) = fl.row( ipol ) ; 763 } 764 crossOK = True ; 765 } 766 } 767 else { 768 spxy.row( ipol ) = real( sp.row( ipol ) ) ; 769 flxy.row( ipol ) = fl.row( ipol ) ; 770 } 771 } 772 } 773 if ( sdh.fluxunit == "" ) { 774 const TableRecord &dataColKeys = mDataCol.keywordSet() ; 775 if ( dataColKeys.isDefined( "UNIT" ) ) 776 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ; 777 else if ( dataColKeys.isDefined( "QuantumUnits" ) ) 778 sdh.fluxunit = dataColKeys.asString( "QuantumUnits" ) ; 779 } 780 } 562 781 563 ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ; 782 564 Block<MEpoch> mTimeB( nrow ) ; 783 565 for ( Int irow = 0 ; irow < nrow ; irow++ ) 784 566 mTimeB[irow] = (*mTimeCol)( irow ) ; 785 ROTableColumn *mIntervalCol = tpoolr->construct( t5, "INTERVAL" ) ;786 ROTableColumn *mFlagRowCol = tpoolr->construct( t5, "FLAG_ROW" ) ;567 // ROTableColumn *mIntervalCol = tpoolr->construct( t5, "INTERVAL" ) ; 568 // ROTableColumn *mFlagRowCol = tpoolr->construct( t5, "FLAG_ROW" ) ; 787 569 Block<Int> sysCalIdx( nrow, -1 ) ; 788 570 if ( isSysCal_ ) { … … 793 575 Int srcType = getSrcType( stateId, tpoolr ) ; 794 576 uInt diridx = 0 ; 795 MDirection::Types dirType ;796 577 uInt wid = 0 ; 797 578 Int pidx = 0 ; … … 825 606 826 607 // FLAGROW 827 *flrRF = (uInt)mFlagRowCol->asBool( irow ) ; 608 // *flrRF = (uInt)mFlagRowCol->asBool( irow ) ; 609 *flrRF = (uInt)asBool( "FLAG_ROW", irow, t5, tpoolr ) ; 828 610 829 611 // SPECTRA, FLAG … … 837 619 838 620 // INTERVAL 839 *intervalRF = (Double)(mIntervalCol->asdouble( irow )) ; 621 // *intervalRF = (Double)(mIntervalCol->asdouble( irow )) ; 622 *intervalRF = asDouble( "INTERVAL", irow, t5, tpoolr ) ; 840 623 841 624 // TSYS … … 863 646 864 647 // DIRECTION, AZEL, SCANRATE 865 if ( getPt_ ) { 866 Vector<Double> dir ; 867 Vector<Double> scanrate ; 868 String refString ; 869 diridx = getDirection( diridx, dir, scanrate, refString, pointtab, mTimeB[irow].get("s").getValue() ) ; 870 MDirection::getType( dirType, refString ) ; 871 MeasFrame mf( mTimeB[irow], mp ) ; 872 //mf.resetEpoch( mTimeB[irow] ) ; 873 //mf.resetDirection( MDirection( MVDirection(dir), dirType ) ) ; 874 if ( refString == "J2000" ) { 875 *dirRF = dir ; 876 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 877 Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ; 878 *azRF = (Float)azel(0) ; 879 *elRF = (Float)azel(1) ; 880 } 881 else if ( refString(0,4) == "AZEL" ) { 882 *azRF = (Float)dir(0) ; 883 *elRF = (Float)dir(1) ; 884 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 885 Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ; 886 *dirRF = newdir ; 887 } 888 else { 889 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 890 Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ; 891 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 892 Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ; 893 *dirRF = newdir ; 894 *azRF = (Float)azel(0) ; 895 *elRF = (Float)azel(1) ; 896 } 897 if ( scanrate.size() != 0 ) { 898 *scrRF = scanrate ; 899 } 900 else { 901 *scrRF = defaultScanrate ; 902 } 903 } 904 else { 905 String ref = md.getRefString() ; 906 //Vector<Double> defaultDir = srcDir ; 907 MDirection::getType( dirType, "J2000" ) ; 908 MeasFrame mf( mTimeB[irow], mp ) ; 909 if ( ref != "J2000" ) { 910 ROScalarMeasColumn<MEpoch> tmCol( pointtab, "TIME" ) ; 911 mf.resetEpoch( tmCol( 0 ) ) ; 912 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 913 defaultDir = toj2000( defaultDir ).getAngle("rad").getValue() ; 914 } 915 mf.resetEpoch( mTimeB[irow] ) ; 916 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 917 Vector<Double> azel = toazel( defaultDir ).getAngle("rad").getValue() ; 918 *azRF = (Float)azel(0) ; 919 *elRF = (Float)azel(1) ; 920 *dirRF = defaultDir ; 921 *scrRF = defaultScanrate ; 922 } 648 Vector<Double> dir ; 649 Vector<Double> azel ; 650 Vector<Double> scanrate ; 651 String refString ; 652 if ( getPt_ ) 653 diridx = getDirection( diridx, dir, azel, scanrate, pointtab, mTimeB[irow], mp ) ; 654 else 655 getSourceDirection( dir, azel, scanrate, mTimeB[irow], mp, delayDir ) ; 656 *dirRF = dir ; 657 *azRF = azel[0] ; 658 *elRF = azel[1] ; 659 *scrRF = scanrate ; 923 660 924 661 // Polarization dependent things … … 927 664 *polnoRF = polnos[ipol] ; 928 665 929 //*spRF = sp.row( ipol ) ;930 //*ucarrRF = fl.row( ipol ) ;931 //*tsysRF = tsys.row( ipol ) ;932 666 spRF.define( sp.row( ipol ) ) ; 933 667 ucarrRF.define( fl.row( ipol ) ) ; … … 942 676 cycle++ ; 943 677 } 944 tpoolr->destroy( mIntervalCol ) ; 945 tpoolr->destroy( mFlagRowCol ) ; 946 947 // time1 = gettimeofday_sec() ; 948 // os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 678 679 //time1 = gettimeofday_sec() ; 680 //os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 949 681 950 682 iter5.next() ; 951 683 } 952 953 954 684 iter4.next() ; 955 685 } 956 957 958 686 iter3.next() ; 959 687 } 960 961 962 688 iter2.next() ; 963 689 } 964 965 966 690 iter1.next() ; 967 691 } … … 1008 732 sdh.freqref = "SOURCE"; 1009 733 } 1010 //if ( sdh.fluxunit == "" ) 734 1011 735 if ( sdh.fluxunit == "" || sdh.fluxunit == "CNTS" ) 1012 736 sdh.fluxunit = "K" ; … … 1043 767 } 1044 768 1045 //double endSec = gettimeofday_sec() ;1046 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;769 //double endSec = gettimeofday_sec() ; 770 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1047 771 } 1048 772 … … 1057 781 Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool ) 1058 782 { 1059 //double startSec = gettimeofday_sec() ;1060 //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;783 //double startSec = gettimeofday_sec() ; 784 //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ; 1061 785 1062 786 MSState statetab = mstable_.state() ; 1063 ROTableColumn *sharedCol ; 1064 sharedCol = tpool->construct( statetab, "OBS_MODE" ) ; 1065 String obsMode = sharedCol->asString( stateId ) ; 1066 tpool->destroy( sharedCol ) ; 1067 sharedCol = tpool->construct( statetab, "SIG" ) ; 1068 Bool sig = sharedCol->asBool( stateId ) ; 1069 tpool->destroy( sharedCol ) ; 1070 sharedCol = tpool->construct( statetab, "REF" ) ; 1071 Bool ref = sharedCol->asBool( stateId ) ; 1072 tpool->destroy( sharedCol ) ; 1073 sharedCol = tpool->construct( statetab, "CAL" ) ; 1074 Double cal = (Double)(sharedCol->asdouble( stateId )) ; 1075 tpool->destroy( sharedCol ) ; 787 String obsMode = asString( "OBS_MODE", stateId, statetab, tpool ) ; 788 Bool sig = asBool( "SIG", stateId, statetab, tpool ) ; 789 Bool ref = asBool( "REF", stateId, statetab, tpool ) ; 790 Double cal = asDouble( "CAL", stateId, statetab, tpool ) ; 1076 791 //os_ << "OBS_MODE = " << obsMode << LogIO::POST ; 1077 792 … … 1222 937 1223 938 //os_ << "srcType = " << srcType << LogIO::POST ; 1224 //double endSec = gettimeofday_sec() ;1225 //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;939 //double endSec = gettimeofday_sec() ; 940 //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1226 941 return srcType ; 1227 942 } … … 1230 945 Block<uInt> MSFiller::getPolNo( Int corrType ) 1231 946 { 1232 //double startSec = gettimeofday_sec() ;1233 //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;947 //double startSec = gettimeofday_sec() ; 948 //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ; 1234 949 Block<uInt> polno( 1 ) ; 1235 950 … … 1261 976 } 1262 977 //os_ << "polno = " << polno << LogIO::POST ; 1263 //double endSec = gettimeofday_sec() ;1264 //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;978 //double endSec = gettimeofday_sec() ; 979 //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1265 980 1266 981 return polno ; … … 1269 984 String MSFiller::getPolType( Int corrType ) 1270 985 { 1271 //double startSec = gettimeofday_sec() ;1272 //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;986 //double startSec = gettimeofday_sec() ; 987 //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ; 1273 988 String poltype = "" ; 1274 989 … … 1282 997 poltype = "linpol" ; 1283 998 1284 //double endSec = gettimeofday_sec() ;1285 //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;999 //double endSec = gettimeofday_sec() ; 1000 //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1286 1001 return poltype ; 1287 1002 } … … 1289 1004 void MSFiller::fillWeather() 1290 1005 { 1291 //double startSec = gettimeofday_sec() ;1292 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;1006 //double startSec = gettimeofday_sec() ; 1007 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ; 1293 1008 1294 1009 if ( !isWeather_ ) { … … 1430 1145 } 1431 1146 //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ; 1432 //double endSec = gettimeofday_sec() ;1433 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1147 //double endSec = gettimeofday_sec() ; 1148 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1434 1149 } 1435 1150 1436 1151 void MSFiller::fillFocus() 1437 1152 { 1438 //double startSec = gettimeofday_sec() ;1439 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;1153 //double startSec = gettimeofday_sec() ; 1154 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ; 1440 1155 // tentative 1441 1156 table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ; 1442 //double endSec = gettimeofday_sec() ;1443 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1157 //double endSec = gettimeofday_sec() ; 1158 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1444 1159 } 1445 1160 1446 1161 void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr ) 1447 1162 { 1448 //double startSec = gettimeofday_sec() ;1449 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;1163 //double startSec = gettimeofday_sec() ; 1164 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ; 1450 1165 1451 1166 if ( !isSysCal_ ) { … … 1502 1217 Table tab = table_->tcal().table() ; 1503 1218 ArrayColumn<Float> tcalCol( tab, "TCAL" ) ; 1504 ROTableColumn *sharedCol ;1505 1219 uInt oldnr = 0 ; 1506 1220 uInt newnr = 0 ; … … 1513 1227 while( !iter0.pastEnd() ) { 1514 1228 Table t0 = iter0.table() ; 1515 sharedCol = tpoolr->construct( t0, "FEED_ID" ) ; 1516 Int feedId = sharedCol->asInt( 0 ) ; 1517 tpoolr->destroy( sharedCol ) ; 1229 Int feedId = asInt( "FEED_ID", 0, t0, tpoolr ) ; 1518 1230 TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ; 1519 1231 while( !iter1.pastEnd() ) { 1520 1232 Table t1 = iter1.table() ; 1521 sharedCol = tpoolr->construct( t1, "SPECTRAL_WINDOW_ID" ) ; 1522 Int spwId = sharedCol->asInt( 0 ) ; 1523 tpoolr->destroy( sharedCol ) ; 1233 Int spwId = asInt( "SPECTRAL_WINDOW_ID", 0, t1, tpoolr ) ; 1524 1234 tmpTcalCol = new ROArrayColumn<Float>( t1, colTcal_ ) ; 1525 1235 ROScalarQuantColumn<Double> scTimeCol( t1, "TIME" ) ; … … 1555 1265 1556 1266 //tcalrec_.print( std::cout ) ; 1557 //double endSec = gettimeofday_sec() ;1558 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1267 //double endSec = gettimeofday_sec() ; 1268 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1559 1269 } 1560 1270 1561 1271 uInt MSFiller::getWeatherId( uInt idx, Double wtime ) 1562 1272 { 1563 //double startSec = gettimeofday_sec() ;1564 //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;1273 //double startSec = gettimeofday_sec() ; 1274 //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ; 1565 1275 uInt nrow = mwTime_.size() ; 1566 1276 if ( nrow < 2 ) … … 1599 1309 //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ; 1600 1310 1601 //double endSec = gettimeofday_sec() ;1602 //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1311 //double endSec = gettimeofday_sec() ; 1312 //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1603 1313 return wid ; 1604 1314 } … … 1606 1316 void MSFiller::getSysCalTime( Vector<MEpoch> &scTime, Vector<Double> &scInterval, Block<MEpoch> &tcol, Block<Int> &tidx ) 1607 1317 { 1608 //double startSec = gettimeofday_sec() ;1609 //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;1318 //double startSec = gettimeofday_sec() ; 1319 //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ; 1610 1320 1611 1321 if ( !isSysCal_ ) … … 1653 1363 } 1654 1364 } 1655 //double endSec = gettimeofday_sec() ;1656 //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;1365 //double endSec = gettimeofday_sec() ; 1366 //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ; 1657 1367 return ; 1658 1368 } … … 1660 1370 Block<uInt> MSFiller::getTcalId( Int fid, Int spwid, MEpoch &t ) 1661 1371 { 1662 //double startSec = gettimeofday_sec() ;1663 //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;1372 //double startSec = gettimeofday_sec() ; 1373 //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ; 1664 1374 //if ( table_->tcal().table().nrow() == 0 ) { 1665 1375 if ( !isSysCal_ ) { … … 1684 1394 tcalids[ipol] = ids[0] + ipol - 1 ; 1685 1395 1686 //double endSec = gettimeofday_sec() ;1687 //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1396 //double endSec = gettimeofday_sec() ; 1397 //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1688 1398 return tcalids ; 1689 1399 } … … 1691 1401 uInt MSFiller::getDirection( uInt idx, Vector<Double> &dir, Vector<Double> &srate, String &ref, MSPointing &tab, Double t ) 1692 1402 { 1693 //double startSec = gettimeofday_sec() ;1694 // os_ << "start MSFiller::getDirection() startSec=" << startSec << LogIO::POST ;1403 //double startSec = gettimeofday_sec() ; 1404 //os_ << "start MSFiller::getDirection1() startSec=" << startSec << LogIO::POST ; 1695 1405 // assume that cols is sorted by TIME 1696 1406 Bool doInterp = False ; … … 1798 1508 } 1799 1509 1800 //double endSec = gettimeofday_sec() ;1801 // os_ << "end MSFiller::getDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1510 //double endSec = gettimeofday_sec() ; 1511 //os_ << "end MSFiller::getDirection1() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1802 1512 return idx ; 1803 1513 } … … 1852 1562 return MFrequency::showType( frame ) ; 1853 1563 } 1564 1565 void MSFiller::reshapeSpectraAndFlagtra( Cube<Float> &sp, 1566 Cube<Bool> &fl, 1567 Table &tab, 1568 Int &npol, 1569 Int &nchan, 1570 Int &nrow, 1571 Vector<Int> &corrtype ) 1572 { 1573 //double startSec = gettimeofday_sec() ; 1574 //os_ << "start MFiller::reshapeSpectraAndFlagtra() startSec=" << startSec << LogIO::POST ; 1575 if ( isFloatData_ ) { 1576 ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ; 1577 ROArrayColumn<Float> mFloatDataCol( tab, "FLOAT_DATA" ) ; 1578 sp = mFloatDataCol.getColumn() ; 1579 fl = mFlagCol.getColumn() ; 1580 } 1581 else if ( isData_ ) { 1582 sp.resize( npol, nchan, nrow ) ; 1583 fl.resize( npol, nchan, nrow ) ; 1584 ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ; 1585 ROArrayColumn<Complex> mDataCol( tab, "DATA" ) ; 1586 for ( Int irow = 0 ; irow < nrow ; irow++ ) { 1587 Bool crossOK = False ; 1588 Matrix<Complex> mSp = mDataCol( irow ) ; 1589 Matrix<Bool> mFl = mFlagCol( irow ) ; 1590 Matrix<Float> spxy = sp.xyPlane( irow ) ; 1591 Matrix<Bool> flxy = fl.xyPlane( irow ) ; 1592 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 1593 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX 1594 || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) { 1595 if ( !crossOK ) { 1596 spxy.row( ipol ) = real( mSp.row( ipol ) ) ; 1597 flxy.row( ipol ) = mFl.row( ipol ) ; 1598 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::RL ) { 1599 spxy.row( ipol+1 ) = imag( mSp.row( ipol ) ) ; 1600 flxy.row( ipol+1 ) = mFl.row( ipol ) ; 1601 } 1602 else { 1603 spxy.row( ipol+1 ) = imag( conj( mSp.row( ipol ) ) ) ; 1604 flxy.row( ipol+1 ) = mFl.row( ipol ) ; 1605 } 1606 crossOK = True ; 1607 } 1608 } 1609 else { 1610 spxy.row( ipol ) = real( mSp.row( ipol ) ) ; 1611 flxy.row( ipol ) = mFl.row( ipol ) ; 1612 } 1613 } 1614 } 1615 } 1616 //double endSec = gettimeofday_sec() ; 1617 //os_ << "end MSFiller::reshapeSpectraAndFlagtra() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1618 } 1619 1620 uInt MSFiller::getDirection( uInt idx, 1621 Vector<Double> &dir, 1622 Vector<Double> &azel, 1623 Vector<Double> &srate, 1624 MSPointing &ptab, 1625 MEpoch &t, 1626 MPosition &antpos ) 1627 { 1628 //double startSec = gettimeofday_sec() ; 1629 //os_ << "start MFiller::getDirection2() startSec=" << startSec << LogIO::POST ; 1630 String refString ; 1631 MDirection::Types dirType ; 1632 uInt diridx = getDirection( idx, dir, srate, refString, ptab, t.get("s").getValue() ) ; 1633 MDirection::getType( dirType, refString ) ; 1634 MeasFrame mf( t, antpos ) ; 1635 if ( refString == "J2000" ) { 1636 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 1637 azel = toazel( dir ).getAngle("rad").getValue() ; 1638 } 1639 else if ( refString(0,4) == "AZEL" ) { 1640 azel = dir.copy() ; 1641 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 1642 dir = toj2000( dir ).getAngle("rad").getValue() ; 1643 } 1644 else { 1645 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 1646 azel = toazel( dir ).getAngle("rad").getValue() ; 1647 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 1648 dir = toj2000( dir ).getAngle("rad").getValue() ; 1649 } 1650 if ( srate.size() == 0 ) 1651 srate = Vector<Double>( 2, 0.0 ) ; 1652 //double endSec = gettimeofday_sec() ; 1653 //os_ << "end MSFiller::getDirection2() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1654 return diridx ; 1655 } 1656 1657 void MSFiller::getSourceDirection( Vector<Double> &dir, 1658 Vector<Double> &azel, 1659 Vector<Double> &srate, 1660 MEpoch &t, 1661 MPosition &antpos, 1662 Vector<MDirection> &srcdir ) 1663 { 1664 //double startSec = gettimeofday_sec() ; 1665 //os_ << "start MFiller::getSourceDirection() startSec=" << startSec << LogIO::POST ; 1666 Vector<Double> defaultScanrate( 2, 0.0 ) ; 1667 Vector<Double> defaultDir = srcdir[0].getAngle( "rad" ).getValue() ; 1668 if ( srcdir.nelements() > 1 ) 1669 defaultScanrate = srcdir[1].getAngle( "rad" ).getValue() ; 1670 String ref = srcdir[0].getRefString() ; 1671 MDirection::Types dirType ; 1672 MDirection::getType( dirType, ref ) ; 1673 MeasFrame mf( t, antpos ) ; 1674 if ( ref != "J2000" ) { 1675 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 1676 dir = toj2000( defaultDir ).getAngle("rad").getValue() ; 1677 } 1678 else 1679 dir = defaultDir ; 1680 if ( ref != "AZELGEO" ) { 1681 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ; 1682 azel = toazel( defaultDir ).getAngle("rad").getValue() ; 1683 } 1684 else 1685 azel = defaultDir ; 1686 srate = defaultScanrate ; 1687 //double endSec = gettimeofday_sec() ; 1688 //os_ << "end MSFiller::getSourceDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1689 } 1690 1691 void MSFiller::initHeader( STHeader &header ) 1692 { 1693 header.nchan = 0 ; 1694 header.npol = 0 ; 1695 header.nif = 0 ; 1696 header.nbeam = 0 ; 1697 header.observer = "" ; 1698 header.project = "" ; 1699 header.obstype = "" ; 1700 header.antennaname = "" ; 1701 header.antennaposition.resize( 0 ) ; 1702 header.equinox = 0.0 ; 1703 header.freqref = "" ; 1704 header.reffreq = -1.0 ; 1705 header.bandwidth = 0.0 ; 1706 header.utc = 0.0 ; 1707 header.fluxunit = "" ; 1708 header.epoch = "" ; 1709 header.poltype = "" ; 1710 } 1711 1712 String MSFiller::asString( String name, 1713 uInt idx, 1714 Table tab, 1715 boost::object_pool<ROTableColumn> *pool ) 1716 { 1717 ROTableColumn *col = pool->construct( tab, name ) ; 1718 String v = col->asString( idx ) ; 1719 pool->destroy( col ) ; 1720 return v ; 1721 } 1722 1723 Bool MSFiller::asBool( String name, 1724 uInt idx, 1725 Table &tab, 1726 boost::object_pool<ROTableColumn> *pool ) 1727 { 1728 ROTableColumn *col = pool->construct( tab, name ) ; 1729 Bool v = col->asBool( idx ) ; 1730 pool->destroy( col ) ; 1731 return v ; 1732 } 1733 1734 uInt MSFiller::asuInt( String name, 1735 uInt idx, 1736 Table &tab, 1737 boost::object_pool<ROTableColumn> *pool ) 1738 { 1739 ROTableColumn *col = pool->construct( tab, name ) ; 1740 uInt v = col->asuInt( idx ) ; 1741 pool->destroy( col ) ; 1742 return v ; 1743 } 1744 1745 Int MSFiller::asInt( String name, 1746 uInt idx, 1747 Table &tab, 1748 boost::object_pool<ROTableColumn> *pool ) 1749 { 1750 ROTableColumn *col = pool->construct( tab, name ) ; 1751 Int v = col->asInt( idx ) ; 1752 pool->destroy( col ) ; 1753 return v ; 1754 } 1755 1756 Float MSFiller::asFloat( String name, 1757 uInt idx, 1758 Table &tab, 1759 boost::object_pool<ROTableColumn> *pool ) 1760 { 1761 ROTableColumn *col = pool->construct( tab, name ) ; 1762 Float v = col->asfloat( idx ) ; 1763 pool->destroy( col ) ; 1764 return v ; 1765 } 1766 1767 Double MSFiller::asDouble( String name, 1768 uInt idx, 1769 Table &tab, 1770 boost::object_pool<ROTableColumn> *pool ) 1771 { 1772 ROTableColumn *col = pool->construct( tab, name ) ; 1773 Double v = col->asdouble( idx ) ; 1774 pool->destroy( col ) ; 1775 return v ; 1776 } 1777 1778 void MSFiller::sourceInfo( Int sourceId, 1779 Int spwId, 1780 String &name, 1781 MDirection &direction, 1782 Vector<casa::Double> &properMotion, 1783 Vector<casa::Double> &restFreqs, 1784 Vector<casa::String> &transitions, 1785 Vector<casa::Double> &sysVels, 1786 boost::object_pool<ROTableColumn> *tpoolr ) 1787 { 1788 //double startSec = gettimeofday_sec() ; 1789 //os_ << "start MFiller::sourceInfo() startSec=" << startSec << LogIO::POST ; 1790 1791 MSSource srctab = mstable_.source() ; 1792 MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ; 1793 if ( srctabSel.nrow() == 0 ) { 1794 srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ; 1795 } 1796 Int numLines = 0 ; 1797 if ( srctabSel.nrow() > 0 ) { 1798 // source name 1799 name = asString( "NAME", 0, srctabSel, tpoolr ) ; 1800 1801 // source proper motion 1802 ROArrayColumn<Double> roArrDCol( srctabSel, "PROPER_MOTION" ) ; 1803 properMotion = roArrDCol( 0 ) ; 1804 1805 // source direction as MDirection object 1806 ROScalarMeasColumn<MDirection> tmpMeasCol( srctabSel, "DIRECTION" ) ; 1807 direction = tmpMeasCol( 0 ) ; 1808 1809 // number of lines 1810 numLines = asInt( "NUM_LINES", 0, srctabSel, tpoolr ) ; 1811 } 1812 else { 1813 name = "" ; 1814 properMotion = Vector<Double>( 2, 0.0 ) ; 1815 direction = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ; 1816 } 1817 1818 restFreqs.resize( numLines ) ; 1819 transitions.resize( numLines ) ; 1820 sysVels.resize( numLines ) ; 1821 if ( numLines > 0 ) { 1822 if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) { 1823 ROArrayQuantColumn<Double> quantArrCol( srctabSel, "REST_FREQUENCY" ) ; 1824 Array< Quantum<Double> > qRestFreqs = quantArrCol( 0 ) ; 1825 for ( int i = 0 ; i < numLines ; i++ ) { 1826 restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ; 1827 } 1828 } 1829 //os_ << "restFreqs = " << restFreqs << LogIO::POST ; 1830 if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) { 1831 ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ; 1832 if ( transitionCol.isDefined( 0 ) ) 1833 transitions = transitionCol( 0 ) ; 1834 //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ; 1835 } 1836 if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) { 1837 ROArrayColumn<Double> roArrDCol( srctabSel, "SYSVEL" ) ; 1838 sysVels = roArrDCol( 0 ) ; 1839 } 1840 } 1841 1842 //double endSec = gettimeofday_sec() ; 1843 //os_ << "end MSFiller::sourceInfo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1844 } 1845 1846 void MSFiller::spectralSetup( Int spwId, 1847 MEpoch &me, 1848 MPosition &mp, 1849 MDirection &md, 1850 Double &refpix, 1851 Double &refval, 1852 Double &increment, 1853 Int &nchan, 1854 String &freqref, 1855 Double &reffreq, 1856 Double &bandwidth, 1857 boost::object_pool<ROTableColumn> *tpoolr ) 1858 { 1859 //double startSec = gettimeofday_sec() ; 1860 //os_ << "start MFiller::spectralSetup() startSec=" << startSec << LogIO::POST ; 1861 1862 MSSpectralWindow spwtab = mstable_.spectralWindow() ; 1863 MeasFrame mf( me, mp, md ) ; 1864 MFrequency::Types freqRef = MFrequency::castType( (uInt)asInt( "MEAS_FREQ_REF", spwId, spwtab, tpoolr ) ) ; 1865 Bool even = False ; 1866 if ( (nchan/2)*2 == nchan ) even = True ; 1867 ROScalarQuantColumn<Double> tmpQuantCol( spwtab, "TOTAL_BANDWIDTH" ) ; 1868 Double totbw = tmpQuantCol( spwId ).getValue( "Hz" ) ; 1869 if ( nchan != 4 ) 1870 bandwidth = max( bandwidth, totbw ) ; 1871 if ( freqref == "" && nchan != 4) 1872 //sdh.freqref = MFrequency::showType( freqRef ) ; 1873 freqref = "LSRK" ; 1874 if ( reffreq == -1.0 && nchan != 4 ) { 1875 tmpQuantCol.attach( spwtab, "REF_FREQUENCY" ) ; 1876 Quantum<Double> qreffreq = tmpQuantCol( spwId ) ; 1877 if ( freqRef == MFrequency::LSRK ) { 1878 reffreq = qreffreq.getValue("Hz") ; 1879 } 1880 else { 1881 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 1882 reffreq = tolsr( qreffreq ).get("Hz").getValue() ; 1883 } 1884 } 1885 Int refchan = nchan / 2 ; 1886 IPosition refip( 1, refchan ) ; 1887 refpix = 0.5*(nchan-1) ; 1888 refval = 0.0 ; 1889 ROArrayQuantColumn<Double> sharedQDArrCol( spwtab, "CHAN_WIDTH" ) ; 1890 increment = sharedQDArrCol( spwId )( refip ).getValue( "Hz" ) ; 1891 // os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ; 1892 sharedQDArrCol.attach( spwtab, "CHAN_FREQ" ) ; 1893 Vector< Quantum<Double> > chanFreqs = sharedQDArrCol( spwId ) ; 1894 if ( nchan > 1 && chanFreqs[0].getValue("Hz") > chanFreqs[1].getValue("Hz") ) 1895 increment *= -1.0 ; 1896 if ( freqRef == MFrequency::LSRK ) { 1897 if ( even ) { 1898 IPosition refip0( 1, refchan-1 ) ; 1899 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 1900 Double refval1 = chanFreqs(refip).getValue("Hz") ; 1901 refval = 0.5 * ( refval0 + refval1 ) ; 1902 } 1903 else { 1904 refval = chanFreqs(refip).getValue("Hz") ; 1905 } 1906 } 1907 else { 1908 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 1909 if ( even ) { 1910 IPosition refip0( 1, refchan-1 ) ; 1911 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 1912 Double refval1 = chanFreqs(refip).getValue("Hz") ; 1913 refval = 0.5 * ( refval0 + refval1 ) ; 1914 refval = tolsr( refval ).get("Hz").getValue() ; 1915 } 1916 else { 1917 refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ; 1918 } 1919 } 1920 1921 //double endSec = gettimeofday_sec() ; 1922 //os_ << "end MSFiller::spectralSetup() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1923 } 1924 1854 1925 } ; 1855 1926 -
trunk/src/MSFiller.h
r2234 r2237 91 91 // get direction for DIRECTION, AZIMUTH, and ELEVATION columns 92 92 casa::uInt getDirection( casa::uInt idx, casa::Vector<casa::Double> &dir, casa::Vector<casa::Double> &srate, casa::String &ref, casa::MSPointing &tab, casa::Double t ) ; 93 casa::uInt getDirection( casa::uInt idx, casa::Vector<casa::Double> &dir, casa::Vector<casa::Double> &azel, casa::Vector<casa::Double> &srate, casa::MSPointing &tab, casa::MEpoch &t, casa::MPosition &antpos ) ; 94 void getSourceDirection( casa::Vector<casa::Double> &dir, casa::Vector<casa::Double> &azel, casa::Vector<casa::Double> &srate, casa::MEpoch &t, casa::MPosition &antpos, casa::Vector<casa::MDirection> &srcdir ) ; 93 95 94 96 // create key for TCAL table … … 103 105 // get frequency frame 104 106 std::string getFrame() ; 105 107 108 // reshape SPECTRA and FLAGTRA 109 void reshapeSpectraAndFlagtra( casa::Cube<casa::Float> &sp, 110 casa::Cube<casa::Bool> &fl, 111 casa::Table &tab, 112 casa::Int &npol, 113 casa::Int &nchan, 114 casa::Int &nrow, 115 casa::Vector<casa::Int> &corrtype ) ; 116 117 // initialize header 118 void initHeader( STHeader &header ) ; 119 120 // get value from table column using object pool 121 casa::String asString( casa::String name, 122 casa::uInt idx, 123 casa::Table tab, 124 boost::object_pool<casa::ROTableColumn> *pool ) ; 125 casa::Bool asBool( casa::String name, 126 casa::uInt idx, 127 casa::Table &tab, 128 boost::object_pool<casa::ROTableColumn> *pool ) ; 129 casa::Int asInt( casa::String name, 130 casa::uInt idx, 131 casa::Table &tab, 132 boost::object_pool<casa::ROTableColumn> *pool ) ; 133 casa::uInt asuInt( casa::String name, 134 casa::uInt idx, 135 casa::Table &tab, 136 boost::object_pool<casa::ROTableColumn> *pool ) ; 137 casa::Float asFloat( casa::String name, 138 casa::uInt idx, 139 casa::Table &tab, 140 boost::object_pool<casa::ROTableColumn> *pool ) ; 141 casa::Double asDouble( casa::String name, 142 casa::uInt idx, 143 casa::Table &tab, 144 boost::object_pool<casa::ROTableColumn> *pool ) ; 145 146 void sourceInfo( casa::Int sourceId, 147 casa::Int spwId, 148 casa::String &name, 149 casa::MDirection &direction, 150 casa::Vector<casa::Double> &properMotion, 151 casa::Vector<casa::Double> &restFreqs, 152 casa::Vector<casa::String> &transitions, 153 casa::Vector<casa::Double> &sysVels, 154 boost::object_pool<casa::ROTableColumn> *pool ) ; 155 156 void spectralSetup( casa::Int spwId, 157 casa::MEpoch &me, 158 casa::MPosition &mp, 159 casa::MDirection &md, 160 casa::Double &refpix, 161 casa::Double &refval, 162 casa::Double &increment, 163 casa::Int &nchan, 164 casa::String &freqref, 165 casa::Double &reffreq, 166 casa::Double &bandwidth, 167 boost::object_pool<casa::ROTableColumn> *pool ) ; 168 106 169 casa::CountedPtr<Scantable> table_ ; 107 170 casa::MeasurementSet mstable_ ;
Note:
See TracChangeset
for help on using the changeset viewer.