- Timestamp:
- 09/12/11 12:07:41 (13 years ago)
- Location:
- trunk/src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/MSFiller.cpp
r2289 r2291 1 1 2 // 2 3 // C++ Interface: MSFiller … … 4 5 // Description: 5 6 // 6 // This class is specific filler for MS format 7 // This class is specific filler for MS format 8 // New version that is implemented using TableVisitor instead of TableIterator 7 9 // 8 // Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 201 010 // Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011 9 11 // 10 12 // Copyright: See COPYING file that comes with this distribution … … 12 14 // 13 15 16 #include <assert.h> 14 17 #include <iostream> 15 18 #include <map> … … 51 54 #include "STHeader.h" 52 55 53 // #include <ctime>54 // #include <sys/time.h>55 56 56 #include "MathUtils.h" 57 57 … … 60 60 61 61 namespace asap { 62 // double MSFiller::gettimeofday_sec() 63 // { 64 // struct timeval tv ; 65 // gettimeofday( &tv, NULL ) ; 66 // return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ; 67 // } 62 63 class BaseMSFillerVisitor: public TableVisitor { 64 uInt lastRecordNo ; 65 Int lastObservationId ; 66 Int lastFeedId ; 67 Int lastFieldId ; 68 Int lastDataDescId ; 69 Int lastScanNo ; 70 Int lastStateId ; 71 Double lastTime ; 72 protected: 73 const Table &table; 74 uInt count; 75 public: 76 BaseMSFillerVisitor(const Table &table) 77 : table(table) 78 { 79 count = 0; 80 } 81 82 virtual void enterObservationId(const uInt recordNo, Int columnValue) { } 83 virtual void leaveObservationId(const uInt recordNo, Int columnValue) { } 84 virtual void enterFeedId(const uInt recordNo, Int columnValue) { } 85 virtual void leaveFeedId(const uInt recordNo, Int columnValue) { } 86 virtual void enterFieldId(const uInt recordNo, Int columnValue) { } 87 virtual void leaveFieldId(const uInt recordNo, Int columnValue) { } 88 virtual void enterDataDescId(const uInt recordNo, Int columnValue) { } 89 virtual void leaveDataDescId(const uInt recordNo, Int columnValue) { } 90 virtual void enterScanNo(const uInt recordNo, Int columnValue) { } 91 virtual void leaveScanNo(const uInt recordNo, Int columnValue) { } 92 virtual void enterStateId(const uInt recordNo, Int columnValue) { } 93 virtual void leaveStateId(const uInt recordNo, Int columnValue) { } 94 virtual void enterTime(const uInt recordNo, Double columnValue) { } 95 virtual void leaveTime(const uInt recordNo, Double columnValue) { } 96 97 virtual Bool visitRecord(const uInt recordNo, 98 const Int ObservationId, 99 const Int feedId, 100 const Int fieldId, 101 const Int dataDescId, 102 const Int scanNo, 103 const Int stateId, 104 const Double time) { return True ; } 105 106 virtual Bool visit(Bool isFirst, const uInt recordNo, 107 const uInt nCols, void const *const colValues[]) { 108 Int observationId, feedId, fieldId, dataDescId, scanNo, stateId; 109 Double time; 110 { // prologue 111 uInt i = 0; 112 { 113 const Int *col = (const Int *)colValues[i++]; 114 observationId = col[recordNo]; 115 } 116 { 117 const Int *col = (const Int *)colValues[i++]; 118 feedId = col[recordNo]; 119 } 120 { 121 const Int *col = (const Int *)colValues[i++]; 122 fieldId = col[recordNo]; 123 } 124 { 125 const Int *col = (const Int *)colValues[i++]; 126 dataDescId = col[recordNo]; 127 } 128 { 129 const Int *col = (const Int *)colValues[i++]; 130 scanNo = col[recordNo]; 131 } 132 { 133 const Int *col = (const Int *)colValues[i++]; 134 stateId = col[recordNo]; 135 } 136 { 137 const Double *col = (const Double *)colValues[i++]; 138 time = col[recordNo]; 139 } 140 assert(nCols == i); 141 } 142 143 if (isFirst) { 144 enterObservationId(recordNo, observationId); 145 enterFeedId(recordNo, feedId); 146 enterFieldId(recordNo, fieldId); 147 enterDataDescId(recordNo, dataDescId); 148 enterScanNo(recordNo, scanNo); 149 enterStateId(recordNo, stateId); 150 enterTime(recordNo, time); 151 } else { 152 if (lastObservationId != observationId) { 153 leaveTime(lastRecordNo, lastTime); 154 leaveStateId(lastRecordNo, lastStateId); 155 leaveScanNo(lastRecordNo, lastScanNo); 156 leaveDataDescId(lastRecordNo, lastDataDescId); 157 leaveFieldId(lastRecordNo, lastFieldId); 158 leaveFeedId(lastRecordNo, lastFeedId); 159 leaveObservationId(lastRecordNo, lastObservationId); 160 161 enterObservationId(recordNo, observationId); 162 enterFeedId(recordNo, feedId); 163 enterFieldId(recordNo, fieldId); 164 enterDataDescId(recordNo, dataDescId); 165 enterScanNo(recordNo, scanNo); 166 enterStateId(recordNo, stateId); 167 enterTime(recordNo, time); 168 } else if (lastFeedId != feedId) { 169 leaveTime(lastRecordNo, lastTime); 170 leaveStateId(lastRecordNo, lastStateId); 171 leaveScanNo(lastRecordNo, lastScanNo); 172 leaveDataDescId(lastRecordNo, lastDataDescId); 173 leaveFieldId(lastRecordNo, lastFieldId); 174 leaveFeedId(lastRecordNo, lastFeedId); 175 176 enterFeedId(recordNo, feedId); 177 enterFieldId(recordNo, fieldId); 178 enterDataDescId(recordNo, dataDescId); 179 enterScanNo(recordNo, scanNo); 180 enterStateId(recordNo, stateId); 181 enterTime(recordNo, time); 182 } else if (lastFieldId != fieldId) { 183 leaveTime(lastRecordNo, lastTime); 184 leaveStateId(lastRecordNo, lastStateId); 185 leaveScanNo(lastRecordNo, lastScanNo); 186 leaveDataDescId(lastRecordNo, lastDataDescId); 187 leaveFieldId(lastRecordNo, lastFieldId); 188 189 enterFieldId(recordNo, fieldId); 190 enterDataDescId(recordNo, dataDescId); 191 enterScanNo(recordNo, scanNo); 192 enterStateId(recordNo, stateId); 193 enterTime(recordNo, time); 194 } else if (lastDataDescId != dataDescId) { 195 leaveTime(lastRecordNo, lastTime); 196 leaveStateId(lastRecordNo, lastStateId); 197 leaveScanNo(lastRecordNo, lastScanNo); 198 leaveDataDescId(lastRecordNo, lastDataDescId); 199 200 enterDataDescId(recordNo, dataDescId); 201 enterScanNo(recordNo, scanNo); 202 enterStateId(recordNo, stateId); 203 enterTime(recordNo, time); 204 } else if (lastScanNo != scanNo) { 205 leaveTime(lastRecordNo, lastTime); 206 leaveStateId(lastRecordNo, lastStateId); 207 leaveScanNo(lastRecordNo, lastScanNo); 208 209 enterScanNo(recordNo, scanNo); 210 enterStateId(recordNo, stateId); 211 enterTime(recordNo, time); 212 } else if (lastStateId != stateId) { 213 leaveTime(lastRecordNo, lastTime); 214 leaveStateId(lastRecordNo, lastStateId); 215 216 enterStateId(recordNo, stateId); 217 enterTime(recordNo, time); 218 } else if (lastTime != time) { 219 leaveTime(lastRecordNo, lastTime); 220 enterTime(recordNo, time); 221 } 222 } 223 count++; 224 Bool result = visitRecord(recordNo, observationId, feedId, fieldId, dataDescId, 225 scanNo, stateId, time); 226 227 { // epilogue 228 lastRecordNo = recordNo; 229 230 lastObservationId = observationId; 231 lastFeedId = feedId; 232 lastFieldId = fieldId; 233 lastDataDescId = dataDescId; 234 lastScanNo = scanNo; 235 lastStateId = stateId; 236 lastTime = time; 237 } 238 return result ; 239 } 240 241 virtual void finish() { 242 if (count > 0) { 243 leaveTime(lastRecordNo, lastTime); 244 leaveStateId(lastRecordNo, lastStateId); 245 leaveScanNo(lastRecordNo, lastScanNo); 246 leaveDataDescId(lastRecordNo, lastDataDescId); 247 leaveFieldId(lastRecordNo, lastFieldId); 248 leaveFeedId(lastRecordNo, lastFeedId); 249 leaveObservationId(lastRecordNo, lastObservationId); 250 } 251 } 252 }; 253 254 class MSFillerVisitor: public BaseMSFillerVisitor, public MSFillerUtils { 255 public: 256 MSFillerVisitor(const Table &from, Scantable &to) 257 : BaseMSFillerVisitor(from), 258 scantable( to ) 259 { 260 antennaId = 0 ; 261 rowidx = 0 ; 262 tablerow = TableRow( scantable.table() ) ; 263 feedEntry = Vector<Int>( 64, -1 ) ; 264 nbeam = 0 ; 265 ifmap.clear() ; 266 const TableDesc &desc = table.tableDesc() ; 267 if ( desc.isColumn( "DATA" ) ) 268 dataColumnName = "DATA" ; 269 else if ( desc.isColumn( "FLOAT_DATA" ) ) 270 dataColumnName = "FLOAT_DATA" ; 271 getpt = False ; 272 isWeather = False ; 273 isSysCal = False ; 274 cycleNo = 0 ; 275 numSysCalRow = 0 ; 276 header = scantable.getHeader() ; 277 fluxUnit( header.fluxunit ) ; 278 279 // MS subtables 280 const TableRecord &hdr = table.keywordSet(); 281 obstab = hdr.asTable( "OBSERVATION" ) ; 282 sctab = hdr.asTable( "SYSCAL" ) ; 283 spwtab = hdr.asTable( "SPECTRAL_WINDOW" ) ; 284 statetab = hdr.asTable( "STATE" ) ; 285 ddtab = hdr.asTable( "DATA_DESCRIPTION" ) ; 286 poltab = hdr.asTable( "POLARIZATION" ) ; 287 fieldtab = hdr.asTable( "FIELD" ) ; 288 anttab = hdr.asTable( "ANTENNA" ) ; 289 srctab = hdr.asTable( "SOURCE" ) ; 290 291 // attach to columns 292 // MS MAIN 293 intervalCol.attach( table, "INTERVAL" ) ; 294 flagRowCol.attach( table, "FLAG_ROW" ) ; 295 flagCol.attach( table, "FLAG" ) ; 296 if ( dataColumnName.compare( "DATA" ) == 0 ) 297 dataCol.attach( table, dataColumnName ) ; 298 else 299 floatDataCol.attach( table, dataColumnName ) ; 300 301 // set dummy epoch 302 mf.set( currentTime ) ; 303 304 // initialize dirType 305 dirType = MDirection::N_Types ; 306 307 // 308 // add rows to scantable 309 // 310 // number of polarization is up to 4 311 uInt addrow = table.nrow() * maxNumPol() ; 312 scantable.table().addRow( addrow ) ; 313 314 // attach to columns 315 // Scantable MAIN 316 TableRecord &r = tablerow.record() ; 317 timeRF.attachToRecord( r, "TIME" ) ; 318 intervalRF.attachToRecord( r, "INTERVAL" ) ; 319 directionRF.attachToRecord( r, "DIRECTION" ) ; 320 azimuthRF.attachToRecord( r, "AZIMUTH" ) ; 321 elevationRF.attachToRecord( r, "ELEVATION" ) ; 322 scanRateRF.attachToRecord( r, "SCANRATE" ) ; 323 weatherIdRF.attachToRecord( r, "WEATHER_ID" ) ; 324 cycleNoRF.attachToRecord( r, "CYCLENO" ) ; 325 flagRowRF.attachToRecord( r, "FLAGROW" ) ; 326 polNoRF.attachToRecord( r, "POLNO" ) ; 327 tcalIdRF.attachToRecord( r, "TCAL_ID" ) ; 328 spectraRF.attachToRecord( r, "SPECTRA" ) ; 329 flagtraRF.attachToRecord( r, "FLAGTRA" ) ; 330 tsysRF.attachToRecord( r, "TSYS" ) ; 331 beamNoRF.attachToRecord( r, "BEAMNO" ) ; 332 ifNoRF.attachToRecord( r, "IFNO" ) ; 333 freqIdRF.attachToRecord( r, "FREQ_ID" ) ; 334 moleculeIdRF.attachToRecord( r, "MOLECULE_ID" ) ; 335 sourceNameRF.attachToRecord( r, "SRCNAME" ) ; 336 sourceProperMotionRF.attachToRecord( r, "SRCPROPERMOTION" ) ; 337 sourceDirectionRF.attachToRecord( r, "SRCDIRECTION" ) ; 338 sourceVelocityRF.attachToRecord( r, "SRCVELOCITY" ) ; 339 focusIdRF.attachToRecord( r, "FOCUS_ID" ) ; 340 fieldNameRF.attachToRecord( r, "FIELDNAME" ) ; 341 sourceTypeRF.attachToRecord( r, "SRCTYPE" ) ; 342 scanNoRF.attachToRecord( r, "SCANNO" ) ; 343 344 // put values 345 RecordFieldPtr<Int> refBeamNoRF( r, "REFBEAMNO" ) ; 346 *refBeamNoRF = -1 ; 347 RecordFieldPtr<Int> fitIdRF( r, "FIT_ID" ) ; 348 *fitIdRF = -1 ; 349 RecordFieldPtr<Float> opacityRF( r, "OPACITY" ) ; 350 *opacityRF = 0.0 ; 351 } 352 353 virtual void enterObservationId(const uInt recordNo, Int columnValue) { 354 //printf("%u: ObservationId: %d\n", recordNo, columnValue); 355 // update header 356 if ( header.observer.empty() ) 357 getScalar( String("OBSERVER"), (uInt)columnValue, obstab, header.observer ) ; 358 if ( header.project.empty() ) 359 getScalar( "PROJECT", (uInt)columnValue, obstab, header.project ) ; 360 if ( header.utc == 0.0 ) { 361 Vector<MEpoch> amp ; 362 getArrayMeas( "TIME_RANGE", (uInt)columnValue, obstab, amp ) ; 363 header.utc = amp[0].get( "d" ).getValue() ; 364 } 365 if ( header.antennaname.empty() ) 366 getScalar( "TELESCOPE_NAME", (uInt)columnValue, obstab, header.antennaname ) ; 367 } 368 virtual void leaveObservationId(const uInt recordNo, Int columnValue) { 369 // update header 370 header.nbeam = max( header.nbeam, (Int)nbeam ) ; 371 372 nbeam = 0 ; 373 feedEntry = -1 ; 374 } 375 virtual void enterFeedId(const uInt recordNo, Int columnValue) { 376 //printf("%u: FeedId: %d\n", recordNo, columnValue); 377 378 // update feed entry 379 if ( allNE( feedEntry, columnValue ) ) { 380 feedEntry[nbeam] = columnValue ; 381 nbeam++ ; 382 } 383 384 // put values 385 *beamNoRF = (uInt)columnValue ; 386 *focusIdRF = (uInt)0 ; 387 } 388 virtual void leaveFeedId(const uInt recordNo, Int columnValue) { 389 Int nelem = (Int)feedEntry.nelements() ; 390 if ( nbeam > nelem ) { 391 feedEntry.resize( nelem+64, True ) ; 392 Slicer slice( IPosition( 1, nelem ), IPosition( 1, feedEntry.nelements()-1 ) ) ; 393 feedEntry( slice ) = -1 ; 394 } 395 } 396 virtual void enterFieldId(const uInt recordNo, Int columnValue) { 397 //printf("%u: FieldId: %d\n", recordNo, columnValue); 398 // update sourceId and fieldName 399 getScalar( "SOURCE_ID", (uInt)columnValue, fieldtab, sourceId ) ; 400 String fieldName ; 401 getScalar( "NAME", (uInt)columnValue, fieldtab, fieldName ) ; 402 fieldName += "__" + String::toString( columnValue ) ; 403 404 // put values 405 *fieldNameRF = fieldName ; 406 } 407 virtual void leaveFieldId(const uInt recordNo, Int columnValue) { 408 sourceId = -1 ; 409 } 410 virtual void enterDataDescId(const uInt recordNo, Int columnValue) { 411 //printf("%u: DataDescId: %d\n", recordNo, columnValue); 412 // update polarization and spectral window ids 413 getScalar( "POLARIZATION_ID", (uInt)columnValue, ddtab, polId ) ; 414 getScalar( "SPECTRAL_WINDOW_ID", (uInt)columnValue, ddtab, spwId ) ; 415 416 // polarization setup 417 getScalar( "NUM_CORR", (uInt)polId, poltab, npol ) ; 418 Vector<Int> corrtype ; 419 getArray( "CORR_TYPE", (uInt)polId, poltab, corrtype ) ; 420 polnos = getPolNos( corrtype ) ; 421 422 // process SOURCE table 423 String sourceName ; 424 Vector<Double> sourcePM, restFreqs, sysVels ; 425 Vector<String> transition ; 426 processSource( sourceId, spwId, sourceName, sourceDir, sourcePM, 427 restFreqs, transition, sysVels ) ; 428 429 // spectral setup 430 uInt freqId ; 431 Double reffreq, bandwidth ; 432 String freqref ; 433 getScalar( "NUM_CHAN", (uInt)spwId, spwtab, nchan ) ; 434 Bool iswvr = (Bool)(nchan == 4) ; 435 map<Int,uInt>::iterator iter = ifmap.find( spwId ) ; 436 if ( iter == ifmap.end() ) { 437 MEpoch me ; 438 getScalarMeas( "TIME", recordNo, table, me ) ; 439 spectralSetup( spwId, me, antpos, sourceDir, 440 freqId, nchan, 441 freqref, reffreq, bandwidth ) ; 442 ifmap.insert( pair<Int,uInt>(spwId,freqId) ) ; 443 } 444 else { 445 freqId = iter->second ; 446 } 447 sp.resize( npol, nchan ) ; 448 fl.resize( npol, nchan ) ; 449 450 451 // molecular setup 452 STMolecules mtab = scantable.molecules() ; 453 uInt molId = mtab.addEntry( restFreqs, transition, transition ) ; 454 455 // process SYSCAL table 456 if ( isSysCal ) 457 processSysCal( spwId ) ; 458 459 // update header 460 if ( !iswvr ) { 461 header.nchan = max( header.nchan, nchan ) ; 462 header.bandwidth = max( header.bandwidth, bandwidth ) ; 463 if ( header.reffreq == -1.0 ) 464 header.reffreq = reffreq ; 465 header.npol = max( header.npol, npol ) ; 466 if ( header.poltype.empty() ) 467 header.poltype = getPolType( corrtype[0] ) ; 468 if ( header.freqref.empty() ) 469 header.freqref = freqref ; 470 } 471 472 // put values 473 *ifNoRF = (uInt)spwId ; 474 *freqIdRF = freqId ; 475 *moleculeIdRF = molId ; 476 *sourceNameRF = sourceName ; 477 sourceProperMotionRF.define( sourcePM ) ; 478 Vector<Double> srcD = sourceDir.getAngle().getValue( "rad" ) ; 479 sourceDirectionRF.define( srcD ) ; 480 if ( !sysVels.empty() ) 481 *sourceVelocityRF = sysVels[0] ; 482 else { 483 *sourceVelocityRF = (Double)0.0 ; 484 } 485 } 486 virtual void leaveDataDescId(const uInt recordNo, Int columnValue) { 487 npol = 0 ; 488 nchan = 0 ; 489 numSysCalRow = 0 ; 490 } 491 virtual void enterScanNo(const uInt recordNo, Int columnValue) { 492 //printf("%u: ScanNo: %d\n", recordNo, columnValue); 493 // put value 494 // scan number is 1-based in MS while 0-based in Scantable 495 *scanNoRF = (uInt)columnValue - 1 ; 496 } 497 virtual void leaveScanNo(const uInt recordNo, Int columnValue) { 498 cycleNo = 0 ; 499 } 500 virtual void enterStateId(const uInt recordNo, Int columnValue) { 501 //printf("%u: StateId: %d\n", recordNo, columnValue); 502 // SRCTYPE 503 Int srcType = getSrcType( columnValue ) ; 504 505 // update header 506 if ( header.obstype.empty() ) 507 getScalar( "OBS_MODE", (uInt)columnValue, statetab, header.obstype ) ; 508 509 // put value 510 *sourceTypeRF = srcType ; 511 } 512 virtual void leaveStateId(const uInt recordNo, Int columnValue) { } 513 virtual void enterTime(const uInt recordNo, Double columnValue) { 514 //printf("%u: Time: %f\n", recordNo, columnValue); 515 currentTime = MEpoch( Quantity( columnValue, "s" ), MEpoch::UTC ) ; 516 517 // DIRECTION, AZEL, and SCANRATE 518 Vector<Double> direction, azel ; 519 Vector<Double> scanrate( 2, 0.0 ) ; 520 if ( getpt ) 521 getDirection( direction, azel, scanrate ) ; 522 else 523 getSourceDirection( direction, azel, scanrate ) ; 524 525 // INTERVAL 526 Double interval = intervalCol.asdouble( recordNo ) ; 527 528 // WEATHER_ID 529 uInt wid = 0 ; 530 if ( isWeather ) 531 wid = getWeatherId() ; 532 533 // put value 534 Double t = currentTime.get( "d" ).getValue() ; 535 *timeRF = t ; 536 *intervalRF = interval ; 537 directionRF.define( direction ) ; 538 *azimuthRF = (Float)azel[0] ; 539 *elevationRF = (Float)azel[1] ; 540 scanRateRF.define( scanrate ) ; 541 *weatherIdRF = wid ; 542 } 543 virtual void leaveTime(const uInt recordNo, Double columnValue) { } 544 virtual Bool visitRecord(const uInt recordNo, 545 const Int observationId, 546 const Int feedId, 547 const Int fieldId, 548 const Int dataDescId, 549 const Int scanNo, 550 const Int stateId, 551 const Double time) 552 { 553 //printf("%u: %d, %d, %d, %d, %d, %d, %f\n", recordNo, 554 //observationId, feedId, fieldId, dataDescId, scanNo, stateId, time); 555 556 // SPECTRA and FLAGTRA 557 //Matrix<Float> sp; 558 //Matrix<uChar> fl; 559 spectraAndFlagtra( recordNo, sp, fl ) ; 560 561 // FLAGROW 562 Bool flr = flagRowCol.asBool( recordNo ) ; 563 564 // TSYS 565 Matrix<Float> tsys ; 566 uInt scIdx = getSysCalIndex() ; 567 if ( numSysCalRow > 0 ) { 568 tsys = sysCalTsysCol( syscalRow[scIdx] ) ; 569 } 570 else { 571 tsys.resize( npol, 1 ) ; 572 tsys = 1.0 ; 573 } 574 575 // TCAL_ID 576 Block<uInt> tcalids( npol, 0 ) ; 577 if ( numSysCalRow > 0 ) { 578 tcalids = getTcalId( syscalTime[scIdx] ) ; 579 } 580 581 // put value 582 *cycleNoRF = cycleNo ; 583 *flagRowRF = (uInt)flr ; 584 585 // for each polarization component 586 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 587 // put value depending on polarization component 588 *polNoRF = polnos[ipol] ; 589 *tcalIdRF = tcalids[ipol] ; 590 spectraRF.define( sp.row( ipol ) ) ; 591 flagtraRF.define( fl.row( ipol ) ) ; 592 tsysRF.define( tsys.row( ipol ) ) ; 593 594 // commit row 595 tablerow.put( rowidx ) ; 596 rowidx++ ; 597 } 598 599 // increment CYCLENO 600 cycleNo++ ; 601 602 return True ; 603 } 604 virtual void finish() 605 { 606 BaseMSFillerVisitor::finish(); 607 //printf("Total: %u\n", count); 608 // remove redundant rows 609 //cout << "filled " << rowidx << " rows out of " << scantable.nrow() << " rows" << endl ; 610 if ( scantable.nrow() > rowidx ) { 611 uInt numRemove = scantable.nrow() - rowidx ; 612 //cout << "numRemove = " << numRemove << endl ; 613 Vector<uInt> rows( numRemove ) ; 614 indgen( rows, rowidx ) ; 615 scantable.table().removeRow( rows ) ; 616 } 617 618 // antenna name and station name 619 String antennaName ; 620 getScalar( "NAME", (uInt)antennaId, anttab, antennaName ) ; 621 String stationName ; 622 getScalar( "STATION", (uInt)antennaId, anttab, stationName ) ; 623 624 // update header 625 header.nif = ifmap.size() ; 626 header.antennaposition = antpos.get( "m" ).getValue() ; 627 if ( header.antennaname.empty() || header.antennaname == antennaName ) 628 header.antennaname = antennaName ; 629 else 630 header.antennaname += "//" + antennaName ; 631 if ( !stationName.empty() && stationName != antennaName ) 632 header.antennaname += "@" + stationName ; 633 if ( header.fluxunit.empty() || header.fluxunit == "CNTS" ) 634 header.fluxunit = "K" ; 635 header.epoch = "UTC" ; 636 header.equinox = 2000.0 ; 637 if (header.freqref == "TOPO") { 638 header.freqref = "TOPOCENT"; 639 } else if (header.freqref == "GEO") { 640 header.freqref = "GEOCENTR"; 641 } else if (header.freqref == "BARY") { 642 header.freqref = "BARYCENT"; 643 } else if (header.freqref == "GALACTO") { 644 header.freqref = "GALACTOC"; 645 } else if (header.freqref == "LGROUP") { 646 header.freqref = "LOCALGRP"; 647 } else if (header.freqref == "CMB") { 648 header.freqref = "CMBDIPOL"; 649 } else if (header.freqref == "REST") { 650 header.freqref = "SOURCE"; 651 } 652 scantable.setHeader( header ) ; 653 } 654 void setAntenna( Int id ) 655 { 656 antennaId = id ; 657 658 Vector< Quantum<Double> > pos ; 659 getArrayQuant( "POSITION", (uInt)antennaId, anttab, pos ) ; 660 antpos = MPosition( MVPosition( pos ), MPosition::ITRF ) ; 661 mf.set( antpos ) ; 662 } 663 void setPointingTable( const Table &tab, String columnToUse="DIRECTION" ) 664 { 665 // input POINTING table must be 666 // 1) selected by antenna 667 // 2) sorted by TIME 668 ROScalarColumn<Double> tcol( tab, "TIME" ) ; 669 ROArrayColumn<Double> dcol( tab, columnToUse ) ; 670 tcol.getColumn( pointingTime ) ; 671 dcol.getColumn( pointingDirection ) ; 672 const TableRecord &rec = dcol.keywordSet() ; 673 String pointingRef = rec.asRecord( "MEASINFO" ).asString( "Ref" ) ; 674 MDirection::getType( dirType, pointingRef ) ; 675 getpt = True ; 676 } 677 void setWeatherTime( const Vector<Double> &t, const Vector<Double> &it ) 678 { 679 isWeather = True ; 680 weatherTime = t ; 681 weatherInterval = it ; 682 } 683 void setSysCalRecord( const Record &r ) 684 //void setSysCalRecord( const map< String,Vector<uInt> > &r ) 685 { 686 isSysCal = True ; 687 syscalRecord = r ; 688 689 const TableDesc &desc = sctab.tableDesc() ; 690 uInt nrow = sctab.nrow() ; 691 syscalRow.resize( nrow ) ; 692 syscalTime.resize( nrow ) ; 693 syscalInterval.resize( nrow ) ; 694 String tsysCol = "NONE" ; 695 Vector<String> tsysCols = stringToVector( "TSYS_SPECTRUM,TSYS" ) ; 696 for ( uInt i = 0 ; i < tsysCols.nelements() ; i++ ) { 697 if ( tsysCol == "NONE" && desc.isColumn( tsysCols[i] ) ) 698 tsysCol = tsysCols[i] ; 699 } 700 sysCalTsysCol.attach( sctab, tsysCol ) ; 701 } 702 STHeader getHeader() { return header ; } 703 uInt getNumBeam() { return nbeam ; } 704 uInt getFilledRowNum() { return rowidx ; } 705 private: 706 void fluxUnit( String &u ) 707 { 708 ROTableColumn col( table, dataColumnName ) ; 709 const TableRecord &rec = col.keywordSet() ; 710 if ( rec.isDefined( "UNIT" ) ) 711 u = rec.asString( "UNIT" ) ; 712 else if ( rec.isDefined( "QuantumUnits" ) ) 713 u = rec.asString( "QuantumUnits" ) ; 714 if ( u.empty() ) 715 u = "K" ; 716 } 717 void processSource( Int sourceId, Int spwId, 718 String &name, MDirection &dir, Vector<Double> &pm, 719 Vector<Double> &rf, Vector<String> &trans, Vector<Double> &vel ) 720 { 721 // find row 722 uInt nrow = srctab.nrow() ; 723 Int idx = -1 ; 724 ROTableRow row( srctab ) ; 725 for ( uInt irow = 0 ; irow < nrow ; irow++ ) { 726 const TableRecord &r = row.get( irow ) ; 727 if ( r.asInt( "SOURCE_ID" ) == sourceId ) { 728 Int tmpSpwId = r.asInt( "SPECTRAL_WINDOW_ID" ) ; 729 if ( tmpSpwId == spwId || tmpSpwId == -1 ) { 730 idx = (Int)irow ; 731 break ; 732 } 733 } 734 } 735 736 // fill 737 Int numLines = 0 ; 738 if ( idx != -1 ) { 739 const TableRecord &r = row.get( idx ) ; 740 name = r.asString( "NAME" ) ; 741 getScalarMeas( "DIRECTION", idx, srctab, dir ) ; 742 pm = r.toArrayDouble( "PROPER_MOTION" ) ; 743 numLines = r.asInt( "NUM_LINES" ) ; 744 } 745 else { 746 name = "" ; 747 pm = Vector<Double>( 2, 0.0 ) ; 748 dir = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ; 749 } 750 if ( !getpt ) { 751 String ref = dir.getRefString() ; 752 MDirection::getType( dirType, ref ) ; 753 } 754 755 rf.resize( numLines ) ; 756 trans.resize( numLines ) ; 757 vel.resize( numLines ) ; 758 if ( numLines > 0 ) { 759 Block<Bool> isDefined = row.getDefined() ; 760 Vector<String> colNames = row.columnNames() ; 761 Vector<Int> indexes( 3, -1 ) ; 762 Vector<String> cols = stringToVector( "REST_FREQUENCY,TRANSITION,SYSVEL" ) ; 763 for ( uInt icol = 0 ; icol < colNames.nelements() ; icol++ ) { 764 if ( anyEQ( indexes, -1 ) ) { 765 for ( uInt jcol = 0 ; jcol < cols.nelements() ; jcol++ ) { 766 if ( colNames[icol] == cols[jcol] ) 767 indexes[jcol] = icol ; 768 } 769 } 770 } 771 if ( indexes[0] != -1 && isDefined[indexes[0]] == True ) { 772 Vector< Quantum<Double> > qrf ; 773 getArrayQuant( "REST_FREQUENCY", idx, srctab, qrf ) ; 774 for ( int i = 0 ; i < numLines ; i++ ) 775 rf[i] = qrf[i].getValue( "Hz" ) ; 776 } 777 if ( indexes[1] != -1 && isDefined[indexes[1]] == True ) { 778 getArray( "TRANSITION", idx, srctab, trans ) ; 779 } 780 if ( indexes[2] != -1 && isDefined[indexes[2]] == True ) { 781 Vector< Quantum<Double> > qsv ; 782 getArrayQuant( "SYSVEL", idx, srctab, qsv ) ; 783 for ( int i = 0 ; i < numLines ; i++ ) 784 vel[i] = qsv[i].getValue( "m/s" ) ; 785 } 786 } 787 } 788 void spectralSetup( Int &spwId, MEpoch &me, MPosition &mp, MDirection &md, 789 uInt &freqId, Int &nchan, 790 String &freqref, Double &reffreq, Double &bandwidth ) 791 { 792 // fill 793 Int measFreqRef ; 794 getScalar( "MEAS_FREQ_REF", spwId, spwtab, measFreqRef ) ; 795 MFrequency::Types freqRef = MFrequency::castType( measFreqRef ) ; 796 //freqref = MFrequency::showType( freqRef ) ; 797 freqref = "LSRK" ; 798 Quantum<Double> q ; 799 getScalarQuant( "TOTAL_BANDWIDTH", spwId, spwtab, q ) ; 800 bandwidth = q.getValue( "Hz" ) ; 801 getScalarQuant( "REF_FREQUENCY", spwId, spwtab, q ) ; 802 reffreq = q.getValue( "Hz" ) ; 803 Double refpix = 0.5 * ( (Double)nchan-1.0 ) ; 804 Int refchan = ( nchan - 1 ) / 2 ; 805 Bool even = (Bool)( nchan % 2 == 0 ) ; 806 Vector< Quantum<Double> > qa ; 807 getArrayQuant( "CHAN_WIDTH", spwId, spwtab, qa ) ; 808 Double increment = qa[refchan].getValue( "Hz" ) ; 809 getArrayQuant( "CHAN_FREQ", spwId, spwtab, qa ) ; 810 if ( nchan == 1 ) { 811 Int netSideband ; 812 getScalar( "NET_SIDEBAND", spwId, spwtab, netSideband ) ; 813 if ( netSideband == 1 ) increment *= -1.0 ; 814 } 815 else { 816 if ( qa[0].getValue( "Hz" ) > qa[1].getValue( "Hz" ) ) 817 increment *= -1.0 ; 818 } 819 Double refval = qa[refchan].getValue( "Hz" ) ; 820 if ( even ) 821 refval = 0.5 * ( refval + qa[refchan+1].getValue( "Hz" ) ) ; 822 if ( freqRef != MFrequency::LSRK ) { 823 MeasFrame mframe( me, mp, md ) ; 824 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mframe ) ) ; 825 refval = tolsr( Quantum<Double>( refval, "Hz" ) ).get( "Hz" ).getValue() ; 826 } 827 828 // add new row to FREQUENCIES 829 Table ftab = scantable.frequencies().table() ; 830 freqId = ftab.nrow() ; 831 ftab.addRow() ; 832 TableRow row( ftab ) ; 833 TableRecord &r = row.record() ; 834 RecordFieldPtr<uInt> idRF( r, "ID" ) ; 835 *idRF = freqId ; 836 RecordFieldPtr<Double> refpixRF( r, "REFPIX" ) ; 837 RecordFieldPtr<Double> refvalRF( r, "REFVAL" ) ; 838 RecordFieldPtr<Double> incrRF( r, "INCREMENT" ) ; 839 *refpixRF = refpix ; 840 *refvalRF = refval ; 841 *incrRF = increment ; 842 row.put( freqId ) ; 843 } 844 void spectraAndFlagtra( uInt recordNo, Matrix<Float> &sp, Matrix<uChar> &fl ) 845 { 846 Matrix<Bool> b = flagCol( recordNo ) ; 847 if ( dataColumnName.compare( "FLOAT_DATA" ) == 0 ) { 848 sp = floatDataCol( recordNo ) ; 849 convertArray( fl, b ) ; 850 } 851 else { 852 Bool notyet = True ; 853 Matrix<Complex> c = dataCol( recordNo ) ; 854 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 855 if ( ( header.poltype == "linear" || header.poltype == "circular" ) 856 && ( polnos[ipol] == 2 || polnos[ipol] == 3 ) ) { 857 if ( notyet ) { 858 Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ; 859 IPosition start( 1, 0 ) ; 860 IPosition end( 1, 2*nchan-1 ) ; 861 IPosition inc( 1, 2 ) ; 862 if ( polnos[ipol] == 2 ) { 863 sp.row( ipol ) = tmp( start, end, inc ) ; 864 Vector<Bool> br = b.row( ipol ) ; 865 Vector<uChar> flr = fl.row( ipol ) ; 866 convertArray( flr, br ) ; 867 start = IPosition( 1, 1 ) ; 868 Int jpol = ipol+1 ; 869 while( polnos[jpol] != 3 && jpol < npol ) 870 jpol++ ; 871 sp.row( jpol ) = tmp( start, end, inc ) ; 872 flr.reference( fl.row( jpol ) ) ; 873 convertArray( flr, br ) ; 874 } 875 else if ( polnos[ipol] == 3 ) { 876 sp.row( ipol ) = sp.row( ipol ) * (Float)(-1.0) ; 877 Int jpol = ipol+1 ; 878 while( polnos[jpol] != 2 && jpol < npol ) 879 jpol++ ; 880 Vector<Bool> br = b.row( ipol ) ; 881 Vector<uChar> flr = fl.row( jpol ) ; 882 sp.row( jpol ) = tmp( start, end, inc ) ; 883 convertArray( flr, br ) ; 884 start = IPosition( 1, 1 ) ; 885 sp.row( ipol ) = tmp( start, end, inc ) * (Float)(-1.0) ; 886 flr.reference( fl.row( ipol ) ) ; 887 convertArray( flr, br ) ; 888 } 889 notyet = False ; 890 } 891 } 892 else { 893 Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ; 894 IPosition start( 1, 0 ) ; 895 IPosition end( 1, 2*nchan-1 ) ; 896 IPosition inc( 1, 2 ) ; 897 sp.row( ipol ) = tmp( start, end, inc ) ; 898 Vector<Bool> br = b.row( ipol ) ; 899 Vector<uChar> flr = fl.row( ipol ) ; 900 convertArray( flr, br ) ; 901 } 902 } 903 } 904 } 905 uInt binarySearch( Vector<Double> &timeList, Double target ) 906 { 907 Int low = 0 ; 908 Int high = timeList.nelements() ; 909 uInt idx = 0 ; 910 911 while ( low <= high ) { 912 idx = (Int)( 0.5 * ( low + high ) ) ; 913 Double t = timeList[idx] ; 914 if ( t < target ) 915 low = idx + 1 ; 916 else if ( t > target ) 917 high = idx - 1 ; 918 else { 919 return idx ; 920 } 921 } 922 923 idx = max( 0, min( low, high ) ) ; 924 return idx ; 925 } 926 void getDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate ) 927 { 928 // @todo At the moment, do binary search every time 929 // if this is bottleneck, frequency of binary search must be reduced 930 Double t = currentTime.get( "s" ).getValue() ; 931 uInt idx = min( binarySearch( pointingTime, t ), pointingTime.nelements()-1 ) ; 932 Matrix<Double> d ; 933 if ( pointingTime[idx] == t ) 934 d = pointingDirection.xyPlane( idx ) ; 935 else if ( pointingTime[idx] < t ) { 936 if ( idx == pointingTime.nelements()-1 ) 937 d = pointingDirection.xyPlane( idx ) ; 938 else 939 d = interp( pointingTime[idx], pointingTime[idx+1], t, 940 pointingDirection.xyPlane( idx ), pointingDirection.xyPlane( idx+1 ) ) ; 941 } 942 else { 943 if ( idx == 0 ) 944 d = pointingDirection.xyPlane( idx ) ; 945 else 946 d = interp( pointingTime[idx-1], pointingTime[idx], t, 947 pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ; 948 } 949 mf.resetEpoch( currentTime ) ; 950 Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ; 951 if ( dirType != MDirection::J2000 ) { 952 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 953 dir = toj2000( tmp ).getAngle( "rad" ).getValue() ; 954 } 955 else { 956 dir = d.column( 0 ) ; 957 } 958 if ( dirType != MDirection::AZELGEO ) { 959 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ; 960 //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 961 azel = toazel( tmp ).getAngle( "rad" ).getValue() ; 962 } 963 else { 964 azel = d.column( 0 ) ; 965 } 966 if ( d.ncolumn() > 1 ) 967 srate = d.column( 1 ) ; 968 } 969 void getSourceDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate ) 970 { 971 dir = sourceDir.getAngle( "rad" ).getValue() ; 972 mf.resetEpoch( currentTime ) ; 973 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ; 974 azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ; 975 if ( dirType != MDirection::J2000 ) { 976 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 977 dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ; 978 } 979 } 980 String detectSeparator( String &s ) 981 { 982 String tmp = s.substr( 0, s.find_first_of( "," ) ) ; 983 Char *separators[] = { ":", "#", ".", "_" } ; 984 uInt nsep = 4 ; 985 for ( uInt i = 0 ; i < nsep ; i++ ) { 986 if ( tmp.find( separators[i] ) != String::npos ) 987 return separators[i] ; 988 } 989 return "" ; 990 } 991 Int getSrcType( Int stateId ) 992 { 993 // get values 994 Bool sig ; 995 getScalar( "SIG", stateId, statetab, sig ) ; 996 Bool ref ; 997 getScalar( "REF", stateId, statetab, ref ) ; 998 Double cal ; 999 getScalar( "CAL", stateId, statetab, cal ) ; 1000 String obsmode ; 1001 getScalar( "OBS_MODE", stateId, statetab, obsmode ) ; 1002 String sep = detectSeparator( obsmode ) ; 1003 1004 Int srcType = SrcType::NOTYPE ; 1005 if ( sep == ":" ) 1006 srcTypeGBT( srcType, sep, obsmode, sig, ref, cal ) ; 1007 else if ( sep == "." || sep == "#" ) 1008 srcTypeALMA( srcType, sep, obsmode ) ; 1009 else if ( sep == "_" ) 1010 srcTypeOldALMA( srcType, sep, obsmode, sig, ref ) ; 1011 else 1012 srcTypeDefault( srcType, sig, ref ) ; 1013 1014 return srcType ; 1015 } 1016 void srcTypeDefault( Int &st, Bool &sig, Bool &ref ) 1017 { 1018 if ( sig ) st = SrcType::SIG ; 1019 else if ( ref ) st = SrcType::REF ; 1020 } 1021 void srcTypeGBT( Int &st, String &sep, String &mode, Bool &sig, Bool &ref, Double &cal ) 1022 { 1023 Int epos = mode.find_first_of( sep ) ; 1024 Int nextpos = mode.find_first_of( sep, epos+1 ) ; 1025 String m1 = mode.substr( 0, epos ) ; 1026 String m2 = mode.substr( epos+1, nextpos-epos-1 ) ; 1027 if ( m1 == "Nod" ) { 1028 st = SrcType::NOD ; 1029 } 1030 else if ( m1 == "OffOn" ) { 1031 if ( m2 == "PSWITCHON" ) st = SrcType::PSON ; 1032 if ( m2 == "PSWITCHOFF" ) st = SrcType::PSOFF ; 1033 } 1034 else { 1035 if ( m2 == "FSWITCH" ) { 1036 if ( sig ) st = SrcType::FSON ; 1037 else if ( ref ) st = SrcType::FSOFF ; 1038 } 1039 } 1040 if ( cal > 0.0 ) { 1041 if ( st == SrcType::NOD ) 1042 st = SrcType::NODCAL ; 1043 else if ( st == SrcType::PSON ) 1044 st = SrcType::PONCAL ; 1045 else if ( st == SrcType::PSOFF ) 1046 st = SrcType::POFFCAL ; 1047 else if ( st == SrcType::FSON ) 1048 st = SrcType::FONCAL ; 1049 else if ( st == SrcType::FSOFF ) 1050 st = SrcType::FOFFCAL ; 1051 else 1052 st = SrcType::CAL ; 1053 } 1054 } 1055 void srcTypeALMA( Int &st, String &sep, String &mode ) 1056 { 1057 Int epos = mode.find_first_of( "," ) ; 1058 String first = mode.substr( 0, epos ) ; 1059 epos = first.find_first_of( sep ) ; 1060 Int nextpos = first.find_first_of( sep, epos+1 ) ; 1061 String m1 = first.substr( 0, epos ) ; 1062 String m2 = first.substr( epos+1, nextpos-epos-1 ) ; 1063 if ( m1.find( "CALIBRATE_" ) == 0 ) { 1064 if ( m2.find( "ON_SOURCE" ) == 0 ) 1065 st = SrcType::PONCAL ; 1066 else if ( m2.find( "OFF_SOURCE" ) == 0 ) 1067 st = SrcType::POFFCAL ; 1068 } 1069 else if ( m1.find( "OBSERVE_TARGET" ) == 0 ) { 1070 if ( m2.find( "ON_SOURCE" ) == 0 ) 1071 st = SrcType::PSON ; 1072 else if ( m2.find( "OFF_SOURCE" ) == 0 ) 1073 st = SrcType::PSOFF ; 1074 } 1075 } 1076 void srcTypeOldALMA( Int &st, String &sep, String &mode, Bool &sig, Bool &ref ) 1077 { 1078 Int epos = mode.find_first_of( "," ) ; 1079 String first = mode.substr( 0, epos ) ; 1080 string substr[4] ; 1081 int numSubstr = split( first, substr, 4, sep ) ; 1082 String m1( substr[0] ) ; 1083 String m2( substr[2] ) ; 1084 if ( numSubstr == 4 ) { 1085 if ( m1.find( "CALIBRATE" ) == 0 ) { 1086 if ( m2.find( "ON" ) == 0 ) 1087 st = SrcType::PONCAL ; 1088 else if ( m2.find( "OFF" ) == 0 ) 1089 st = SrcType::POFFCAL ; 1090 } 1091 else if ( m1.find( "OBSERVE" ) == 0 ) { 1092 if ( m2.find( "ON" ) == 0 ) 1093 st = SrcType::PSON ; 1094 else if ( m2.find( "OFF" ) == 0 ) 1095 st = SrcType::PSOFF ; 1096 } 1097 } 1098 else { 1099 if ( sig ) st = SrcType::SIG ; 1100 else if ( ref ) st = SrcType::REF ; 1101 } 1102 } 1103 Block<uInt> getPolNos( Vector<Int> &corr ) 1104 { 1105 Block<uInt> polnos( npol ) ; 1106 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 1107 if ( corr[ipol] == Stokes::I || corr[ipol] == Stokes::RR || corr[ipol] == Stokes::XX ) 1108 polnos[ipol] = 0 ; 1109 else if ( corr[ipol] == Stokes::Q || corr[ipol] == Stokes::LL || corr[ipol] == Stokes::YY ) 1110 polnos[ipol] = 1 ; 1111 else if ( corr[ipol] == Stokes::U || corr[ipol] == Stokes::RL || corr[ipol] == Stokes::XY ) 1112 polnos[ipol] = 2 ; 1113 else if ( corr[ipol] == Stokes::V || corr[ipol] == Stokes::LR || corr[ipol] == Stokes::YX ) 1114 polnos[ipol] = 3 ; 1115 } 1116 return polnos ; 1117 } 1118 String getPolType( Int &corr ) 1119 { 1120 String poltype = "" ; 1121 if ( corr == Stokes::I || corr == Stokes::Q || corr == Stokes::U || corr == Stokes::V ) 1122 poltype = "stokes" ; 1123 else if ( corr == Stokes::XX || corr == Stokes::YY || corr == Stokes::XY || corr == Stokes::YX ) 1124 poltype = "linear" ; 1125 else if ( corr == Stokes::RR || corr == Stokes::LL || corr == Stokes::RL || corr == Stokes::LR ) 1126 poltype = "circular" ; 1127 else if ( corr == Stokes::Plinear || corr == Stokes::Pangle ) 1128 poltype = "linpol" ; 1129 return poltype ; 1130 } 1131 uInt getWeatherId() 1132 { 1133 // if only one row, return 0 1134 if ( weatherTime.nelements() == 1 ) 1135 return 0 ; 1136 1137 // @todo At the moment, do binary search every time 1138 // if this is bottleneck, frequency of binary search must be reduced 1139 Double t = currentTime.get( "s" ).getValue() ; 1140 uInt idx = min( binarySearch( weatherTime, t ), weatherTime.nelements()-1 ) ; 1141 if ( weatherTime[idx] < t ) { 1142 if ( idx != weatherTime.nelements()-1 ) { 1143 if ( weatherTime[idx+1] - t < 0.5 * weatherInterval[idx+1] ) 1144 idx++ ; 1145 } 1146 } 1147 else if ( weatherTime[idx] > t ) { 1148 if ( idx != 0 ) { 1149 if ( weatherTime[idx] - t > 0.5 * weatherInterval[idx] ) 1150 idx-- ; 1151 } 1152 } 1153 return idx ; 1154 } 1155 void processSysCal( Int &spwId ) 1156 { 1157 // get feedId from row 1158 Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ; 1159 1160 uInt nrow = sctab.nrow() ; 1161 ROScalarColumn<Int> col( sctab, "ANTENNA_ID" ) ; 1162 Vector<Int> aids = col.getColumn() ; 1163 col.attach( sctab, "FEED_ID" ) ; 1164 Vector<Int> fids = col.getColumn() ; 1165 col.attach( sctab, "SPECTRAL_WINDOW_ID" ) ; 1166 Vector<Int> sids = col.getColumn() ; 1167 ROScalarColumn<Double> timeCol( sctab, "TIME" ) ; 1168 ROScalarColumn<Double> intCol( sctab, "INTERVAL" ) ; 1169 for ( uInt irow = 0 ; irow < nrow ; irow++ ) { 1170 if ( aids[irow] == antennaId 1171 && fids[irow] == feedId 1172 && sids[irow] == spwId ) { 1173 syscalRow[numSysCalRow] = irow ; 1174 syscalTime[numSysCalRow] = timeCol( irow ) ; 1175 syscalInterval[numSysCalRow] = intCol( irow ) ; 1176 numSysCalRow++ ; 1177 } 1178 } 1179 } 1180 uInt getSysCalIndex() 1181 { 1182 // if only one row, return 0 1183 if ( numSysCalRow == 1 || !isSysCal ) 1184 return 0 ; 1185 1186 // @todo At the moment, do binary search every time 1187 // if this is bottleneck, frequency of binary search must be reduced 1188 Double t = currentTime.get( "s" ).getValue() ; 1189 Vector<Double> tslice = syscalTime( Slice(0, numSysCalRow) ) ; 1190 uInt idx = min( binarySearch( tslice, t ), numSysCalRow-1 ) ; 1191 if ( syscalTime[idx] < t ) { 1192 if ( idx != numSysCalRow-1 ) { 1193 if ( syscalTime[idx+1] - t < 0.5 * syscalInterval[idx+1] ) 1194 idx++ ; 1195 } 1196 } 1197 else if ( syscalTime[idx] > t ) { 1198 if ( idx != 0 ) { 1199 if ( syscalTime[idx] - t > 0.5 * syscalInterval[idx] ) 1200 idx-- ; 1201 } 1202 } 1203 return idx ; 1204 } 1205 Block<uInt> getTcalId( Double &t ) 1206 { 1207 // return 0 if no SysCal table 1208 if ( !isSysCal ) { 1209 return Block<uInt>( 4, 0 ) ; 1210 } 1211 1212 // get feedId from row 1213 Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ; 1214 1215 // key 1216 String key = keyTcal( feedId, spwId, t ) ; 1217 1218 // retrieve ids 1219 Vector<uInt> ids = syscalRecord.asArrayuInt( key ) ; 1220 //Vector<uInt> ids = syscalRecord[key] ; 1221 uInt np = ids[1] - ids[0] + 1 ; 1222 Block<uInt> tcalids( np ) ; 1223 if ( np > 0 ) { 1224 tcalids[0] = ids[0] ; 1225 if ( np > 1 ) { 1226 tcalids[1] = ids[1] ; 1227 for ( uInt ip = 2 ; ip < np ; ip++ ) 1228 tcalids[ip] = ids[0] + ip - 1 ; 1229 } 1230 } 1231 return tcalids ; 1232 } 1233 uInt maxNumPol() 1234 { 1235 ROScalarColumn<Int> numCorrCol( poltab, "NUM_CORR" ) ; 1236 return max( numCorrCol.getColumn() ) ; 1237 } 1238 1239 Scantable &scantable; 1240 Int antennaId; 1241 uInt rowidx; 1242 String dataColumnName; 1243 TableRow tablerow; 1244 STHeader header; 1245 Vector<Int> feedEntry; 1246 uInt nbeam; 1247 Int npol; 1248 Int nchan; 1249 Int sourceId; 1250 Int polId; 1251 Int spwId; 1252 uInt cycleNo; 1253 MDirection sourceDir; 1254 MPosition antpos; 1255 MEpoch currentTime; 1256 map<Int,uInt> ifmap; 1257 Block<uInt> polnos; 1258 Bool getpt; 1259 Vector<Double> pointingTime; 1260 Cube<Double> pointingDirection; 1261 MDirection::Types dirType; 1262 Bool isWeather; 1263 Vector<Double> weatherTime; 1264 Vector<Double> weatherInterval; 1265 Bool isSysCal; 1266 Record syscalRecord; 1267 //map< String,Vector<uInt> > syscalRecord; 1268 uInt numSysCalRow ; 1269 Vector<uInt> syscalRow; 1270 Vector<Double> syscalTime; 1271 Vector<Double> syscalInterval; 1272 //String tsysCol; 1273 //String tcalCol; 1274 1275 // MS subtables 1276 Table obstab; 1277 Table sctab; 1278 Table spwtab; 1279 Table statetab; 1280 Table ddtab; 1281 Table poltab; 1282 Table fieldtab; 1283 Table anttab; 1284 Table srctab; 1285 Matrix<Float> sp; 1286 Matrix<uChar> fl; 1287 MeasFrame mf; 1288 1289 // MS MAIN columns 1290 ROTableColumn intervalCol; 1291 ROTableColumn flagRowCol; 1292 ROArrayColumn<Float> floatDataCol; 1293 ROArrayColumn<Complex> dataCol; 1294 ROArrayColumn<Bool> flagCol; 1295 1296 // MS SYSCAL columns 1297 ROArrayColumn<Float> sysCalTsysCol; 1298 1299 // Scantable MAIN columns 1300 RecordFieldPtr<Double> timeRF,intervalRF,sourceVelocityRF; 1301 RecordFieldPtr< Vector<Double> > directionRF,scanRateRF, 1302 sourceProperMotionRF,sourceDirectionRF; 1303 RecordFieldPtr<Float> azimuthRF,elevationRF; 1304 RecordFieldPtr<uInt> weatherIdRF,cycleNoRF,flagRowRF,polNoRF,tcalIdRF, 1305 ifNoRF,freqIdRF,moleculeIdRF,beamNoRF,focusIdRF,scanNoRF; 1306 RecordFieldPtr< Vector<Float> > spectraRF,tsysRF; 1307 RecordFieldPtr< Vector<uChar> > flagtraRF; 1308 RecordFieldPtr<String> sourceNameRF,fieldNameRF; 1309 RecordFieldPtr<Int> sourceTypeRF; 1310 }; 1311 1312 class BaseTcalVisitor: public TableVisitor { 1313 uInt lastRecordNo ; 1314 Int lastAntennaId ; 1315 Int lastFeedId ; 1316 Int lastSpwId ; 1317 Double lastTime ; 1318 protected: 1319 const Table &table; 1320 uInt count; 1321 public: 1322 BaseTcalVisitor(const Table &table) 1323 : table(table) 1324 { 1325 count = 0; 1326 } 1327 1328 virtual void enterAntennaId(const uInt recordNo, Int columnValue) { } 1329 virtual void leaveAntennaId(const uInt recordNo, Int columnValue) { } 1330 virtual void enterFeedId(const uInt recordNo, Int columnValue) { } 1331 virtual void leaveFeedId(const uInt recordNo, Int columnValue) { } 1332 virtual void enterSpwId(const uInt recordNo, Int columnValue) { } 1333 virtual void leaveSpwId(const uInt recordNo, Int columnValue) { } 1334 virtual void enterTime(const uInt recordNo, Double columnValue) { } 1335 virtual void leaveTime(const uInt recordNo, Double columnValue) { } 1336 1337 virtual Bool visitRecord(const uInt recordNo, 1338 const Int antennaId, 1339 const Int feedId, 1340 const Int spwId, 1341 const Double time) { return True ; } 1342 1343 virtual Bool visit(Bool isFirst, const uInt recordNo, 1344 const uInt nCols, void const *const colValues[]) { 1345 Int antennaId, feedId, spwId; 1346 Double time; 1347 { // prologue 1348 uInt i = 0; 1349 { 1350 const Int *col = (const Int *)colValues[i++]; 1351 antennaId = col[recordNo]; 1352 } 1353 { 1354 const Int *col = (const Int *)colValues[i++]; 1355 feedId = col[recordNo]; 1356 } 1357 { 1358 const Int *col = (const Int *)colValues[i++]; 1359 spwId = col[recordNo]; 1360 } 1361 { 1362 const Double *col = (const Double *)colValues[i++]; 1363 time = col[recordNo]; 1364 } 1365 assert(nCols == i); 1366 } 1367 1368 if (isFirst) { 1369 enterAntennaId(recordNo, antennaId); 1370 enterFeedId(recordNo, feedId); 1371 enterSpwId(recordNo, spwId); 1372 enterTime(recordNo, time); 1373 } else { 1374 if ( lastAntennaId != antennaId ) { 1375 leaveTime(lastRecordNo, lastTime); 1376 leaveSpwId(lastRecordNo, lastSpwId); 1377 leaveFeedId(lastRecordNo, lastFeedId); 1378 leaveAntennaId(lastRecordNo, lastAntennaId); 1379 1380 enterAntennaId(recordNo, antennaId); 1381 enterFeedId(recordNo, feedId); 1382 enterSpwId(recordNo, spwId); 1383 enterTime(recordNo, time); 1384 } 1385 else if (lastFeedId != feedId) { 1386 leaveTime(lastRecordNo, lastTime); 1387 leaveSpwId(lastRecordNo, lastSpwId); 1388 leaveFeedId(lastRecordNo, lastFeedId); 1389 1390 enterFeedId(recordNo, feedId); 1391 enterSpwId(recordNo, spwId); 1392 enterTime(recordNo, time); 1393 } else if (lastSpwId != spwId) { 1394 leaveTime(lastRecordNo, lastTime); 1395 leaveSpwId(lastRecordNo, lastSpwId); 1396 1397 enterSpwId(recordNo, spwId); 1398 enterTime(recordNo, time); 1399 } else if (lastTime != time) { 1400 leaveTime(lastRecordNo, lastTime); 1401 enterTime(recordNo, time); 1402 } 1403 } 1404 count++; 1405 Bool result = visitRecord(recordNo, antennaId, feedId, spwId, time); 1406 1407 { // epilogue 1408 lastRecordNo = recordNo; 1409 1410 lastAntennaId = antennaId; 1411 lastFeedId = feedId; 1412 lastSpwId = spwId; 1413 lastTime = time; 1414 } 1415 return result ; 1416 } 1417 1418 virtual void finish() { 1419 if (count > 0) { 1420 leaveTime(lastRecordNo, lastTime); 1421 leaveSpwId(lastRecordNo, lastSpwId); 1422 leaveFeedId(lastRecordNo, lastFeedId); 1423 leaveAntennaId(lastRecordNo, lastAntennaId); 1424 } 1425 } 1426 }; 1427 1428 class TcalVisitor: public BaseTcalVisitor, public MSFillerUtils { 1429 public: 1430 TcalVisitor(const Table &table, Table &tcaltab, Record &r, Int aid ) 1431 //TcalVisitor(const Table &table, Table &tcaltab, map< String,Vector<uInt> > &r, Int aid ) 1432 : BaseTcalVisitor( table ), 1433 tcal(tcaltab), 1434 rec(r), 1435 antenna(aid) 1436 { 1437 process = False ; 1438 rowidx = 0 ; 1439 1440 // attach to SYSCAL columns 1441 timeCol.attach( table, "TIME" ) ; 1442 1443 // add rows 1444 uInt addrow = table.nrow() * 4 ; 1445 tcal.addRow( addrow ) ; 1446 1447 // attach to TCAL columns 1448 row = TableRow( tcal ) ; 1449 TableRecord &trec = row.record() ; 1450 idRF.attachToRecord( trec, "ID" ) ; 1451 timeRF.attachToRecord( trec, "TIME" ) ; 1452 tcalRF.attachToRecord( trec, "TCAL" ) ; 1453 } 1454 1455 virtual void enterAntennaId(const uInt recordNo, Int columnValue) { 1456 if ( columnValue == antenna ) 1457 process = True ; 1458 } 1459 virtual void leaveAntennaId(const uInt recordNo, Int columnValue) { 1460 process = False ; 1461 } 1462 virtual void enterFeedId(const uInt recordNo, Int columnValue) { } 1463 virtual void leaveFeedId(const uInt recordNo, Int columnValue) { } 1464 virtual void enterSpwId(const uInt recordNo, Int columnValue) { } 1465 virtual void leaveSpwId(const uInt recordNo, Int columnValue) { } 1466 virtual void enterTime(const uInt recordNo, Double columnValue) { 1467 qtime = timeCol( recordNo ) ; 1468 } 1469 virtual void leaveTime(const uInt recordNo, Double columnValue) { } 1470 virtual Bool visitRecord(const uInt recordNo, 1471 const Int antennaId, 1472 const Int feedId, 1473 const Int spwId, 1474 const Double time) 1475 { 1476 //cout << "(" << recordNo << "," << antennaId << "," << feedId << "," << spwId << ")" << endl ; 1477 if ( process ) { 1478 String sTime = MVTime( qtime ).string( MVTime::YMD ) ; 1479 *timeRF = sTime ; 1480 uInt oldidx = rowidx ; 1481 Matrix<Float> subtcal = tcalCol( recordNo ) ; 1482 Vector<uInt> idminmax( 2 ) ; 1483 for ( uInt ipol = 0 ; ipol < subtcal.nrow() ; ipol++ ) { 1484 *idRF = rowidx ; 1485 tcalRF.define( subtcal.row( ipol ) ) ; 1486 1487 // commit row 1488 row.put( rowidx ) ; 1489 rowidx++ ; 1490 } 1491 1492 idminmax[0] = oldidx ; 1493 idminmax[1] = rowidx - 1 ; 1494 1495 String key = keyTcal( feedId, spwId, sTime ) ; 1496 rec.define( key, idminmax ) ; 1497 //rec[key] = idminmax ; 1498 } 1499 return True ; 1500 } 1501 virtual void finish() 1502 { 1503 BaseTcalVisitor::finish() ; 1504 1505 if ( tcal.nrow() > rowidx ) { 1506 uInt numRemove = tcal.nrow() - rowidx ; 1507 //cout << "numRemove = " << numRemove << endl ; 1508 Vector<uInt> rows( numRemove ) ; 1509 indgen( rows, rowidx ) ; 1510 tcal.removeRow( rows ) ; 1511 } 1512 1513 } 1514 void setTcalColumn( String &col ) 1515 { 1516 //colName = col ; 1517 tcalCol.attach( table, col ) ; 1518 } 1519 private: 1520 Table &tcal; 1521 Record &rec; 1522 //map< String,Vector<uInt> > &rec; 1523 Int antenna; 1524 uInt rowidx; 1525 Bool process; 1526 Quantum<Double> qtime; 1527 TableRow row; 1528 String colName; 1529 1530 // MS SYSCAL columns 1531 ROScalarQuantColumn<Double> timeCol; 1532 ROArrayColumn<Float> tcalCol; 1533 1534 // TCAL columns 1535 RecordFieldPtr<uInt> idRF; 1536 RecordFieldPtr<String> timeRF; 1537 RecordFieldPtr< Vector<Float> > tcalRF; 1538 }; 68 1539 69 1540 MSFiller::MSFiller( casa::CountedPtr<Scantable> stable ) … … 122 1593 123 1594 MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ; 124 //mstable_ = (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_125 // && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ;126 1595 tablename_ = tmpMS->tableName() ; 127 1596 if ( antenna_ == -1 && antennaStr_.size() > 0 ) { … … 139 1608 mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_ 140 1609 && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ; 141 // stringstream ss ; 142 // ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 == " << antenna_ ; 143 // String taql( ss.str() ) ; 144 // mstable_ = MeasurementSet( tableCommand( taql, *tmpMS ) ) ; 1610 145 1611 delete tmpMS ; 146 1612 … … 156 1622 void MSFiller::fill() 157 1623 { 158 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;159 1624 //double startSec = mathutil::gettimeofday_sec() ; 160 1625 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ; 161 1626 162 //double time0 = mathutil::gettimeofday_sec() ; 163 //os_ << "start init fill: " << time0 << LogIO::POST ; 1627 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ; 164 1628 165 1629 // Initialize header 166 1630 STHeader sdh ; 167 1631 initHeader( sdh ) ; 1632 table_->setHeader( sdh ) ; 168 1633 169 1634 // check if optional table exists 170 //const TableRecord msrec = tablesel_.keywordSet() ; 171 const TableRecord msrec = mstable_.keywordSet() ; 1635 const TableRecord &msrec = mstable_.keywordSet() ; 172 1636 isDoppler_ = msrec.isDefined( "DOPPLER" ) ; 173 1637 if ( isDoppler_ ) … … 199 1663 isWeather_ = False ; 200 1664 201 // Access to MS subtables 202 MSField fieldtab = mstable_.field() ; 203 MSPolarization poltab = mstable_.polarization() ; 204 MSDataDescription ddtab = mstable_.dataDescription() ; 205 MSObservation obstab = mstable_.observation() ; 206 MSSource srctab = mstable_.source() ; 207 MSSpectralWindow spwtab = mstable_.spectralWindow() ; 208 MSSysCal caltab = mstable_.sysCal() ; 209 if ( caltab.nrow() == 0 ) 210 isSysCal_ = False ; 211 else { 1665 // column name for Tsys and Tcal 1666 if ( isSysCal_ ) { 1667 const MSSysCal &caltab = mstable_.sysCal() ; 212 1668 if ( !caltab.tableDesc().isColumn( colTcal_ ) ) { 213 1669 colTcal_ = "TCAL" ; … … 221 1677 } 222 1678 } 223 // colTcal_ = "TCAL" ; 224 // colTsys_ = "TSYS" ; 225 MSPointing pointtab = mstable_.pointing() ; 226 if ( mstable_.weather().nrow() == 0 ) 227 isWeather_ = False ; 228 MSState stattab = mstable_.state() ; 229 MSAntenna anttab = mstable_.antenna() ; 230 231 // TEST 232 // memory allocation by boost::object_pool 233 boost::object_pool<ROTableColumn> *tpoolr = new boost::object_pool<ROTableColumn> ; 234 // 1679 else { 1680 colTcal_ = "NONE" ; 1681 colTsys_ = "NONE" ; 1682 } 1683 1684 // Access to MS subtables 1685 //MSField &fieldtab = mstable_.field() ; 1686 //MSPolarization &poltab = mstable_.polarization() ; 1687 //MSDataDescription &ddtab = mstable_.dataDescription() ; 1688 //MSObservation &obstab = mstable_.observation() ; 1689 //MSSource &srctab = mstable_.source() ; 1690 //MSSpectralWindow &spwtab = mstable_.spectralWindow() ; 1691 //MSSysCal &caltab = mstable_.sysCal() ; 1692 MSPointing &pointtab = mstable_.pointing() ; 1693 //MSState &stattab = mstable_.state() ; 1694 //MSAntenna &anttab = mstable_.antenna() ; 235 1695 236 1696 // SUBTABLES: FREQUENCIES … … 247 1707 248 1708 // SUBTABLES: TCAL 249 fillTcal( tpoolr) ;1709 fillTcal() ; 250 1710 251 1711 // SUBTABLES: FIT … … 255 1715 //fillHistory() ; 256 1716 257 // MAIN 258 // Iterate over several ids 259 map<Int, uInt> ifmap ; // (IFNO, FREQ_ID) pair 260 ROArrayQuantColumn<Double> *sharedQDArrCol = new ROArrayQuantColumn<Double>( anttab, "POSITION" ) ; 261 Vector< Quantum<Double> > antpos = (*sharedQDArrCol)( antenna_ ) ; 262 delete sharedQDArrCol ; 263 MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ; 264 Vector<Double> pt ; 265 ROArrayColumn<Double> pdcol ; 266 Vector<Double> defaultScanrate( 2, 0.0 ) ; 267 if ( getPt_ ) { 268 pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ; 269 ROScalarColumn<Double> ptcol( pointtab, "TIME" ) ; 270 ptcol.getColumn( pt ) ; 271 TableRecord trec = ptcol.keywordSet() ; 272 String tUnit = trec.asArrayString( "QuantumUnits" ).data()[0] ; 273 if ( tUnit == "d" ) 274 pt *= 86400.0 ; 275 pdcol.attach( pointtab, "DIRECTION" ) ; 276 } 277 String stationName = asString( "STATION", antenna_, anttab, tpoolr ) ; 278 String antennaName = asString( "NAME", antenna_, anttab, tpoolr ) ; 279 sdh.antennaposition.resize( 3 ) ; 280 for ( int i = 0 ; i < 3 ; i++ ) 281 sdh.antennaposition[i] = antpos[i].getValue( "m" ) ; 282 String telescopeName = "" ; 283 284 //double time1 = mathutil::gettimeofday_sec() ; 285 //os_ << "end init fill: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 286 287 // row based 288 Table &stab = table_->table() ; 289 TableRow row( stab ) ; 290 TableRecord &trec = row.record() ; 291 RecordFieldPtr< Array<Float> > spRF( trec, "SPECTRA" ) ; 292 RecordFieldPtr< Array<uChar> > ucarrRF( trec, "FLAGTRA" ) ; 293 RecordFieldPtr<Double> timeRF( trec, "TIME" ) ; 294 RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ; 295 RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ; 296 RecordFieldPtr< Array<Double> > dirRF( trec, "DIRECTION" ) ; 297 RecordFieldPtr<Float> azRF( trec, "AZIMUTH" ) ; 298 RecordFieldPtr<Float> elRF( trec, "ELEVATION" ) ; 299 RecordFieldPtr< Array<Double> > scrRF( trec, "SCANRATE" ) ; 300 RecordFieldPtr<uInt> cycleRF( trec, "CYCLENO" ) ; 301 RecordFieldPtr<uInt> flrRF( trec, "FLAGROW" ) ; 302 RecordFieldPtr<uInt> tcalidRF( trec, "TCAL_ID" ) ; 303 RecordFieldPtr<uInt> widRF( trec, "WEATHER_ID" ) ; 304 RecordFieldPtr<uInt> polnoRF( trec, "POLNO" ) ; 305 RecordFieldPtr<Int> refbRF( trec, "REFBEAMNO" ) ; 306 RecordFieldPtr<Int> fitidRF( trec, "FIT_ID" ) ; 307 RecordFieldPtr<Float> tauRF( trec, "OPACITY" ) ; 308 RecordFieldPtr<uInt> beamRF( trec, "BEAMNO" ) ; 309 RecordFieldPtr<uInt> focusidRF( trec, "FOCUS_ID" ) ; 310 RecordFieldPtr<uInt> ifnoRF( trec, "IFNO" ) ; 311 RecordFieldPtr<uInt> molidRF( trec, "MOLECULE_ID" ) ; 312 RecordFieldPtr<uInt> freqidRF( trec, "FREQ_ID" ) ; 313 RecordFieldPtr<uInt> scanRF( trec, "SCANNO" ) ; 314 RecordFieldPtr<Int> srctypeRF( trec, "SRCTYPE" ) ; 315 RecordFieldPtr<String> srcnameRF( trec, "SRCNAME" ) ; 316 RecordFieldPtr<String> fieldnameRF( trec, "FIELDNAME" ) ; 317 RecordFieldPtr< Array<Double> > srcpmRF( trec, "SRCPROPERMOTION" ) ; 318 RecordFieldPtr< Array<Double> > srcdirRF( trec, "SRCDIRECTION" ) ; 319 RecordFieldPtr<Double> sysvelRF( trec, "SRCVELOCITY" ) ; 320 321 // REFBEAMNO 322 *refbRF = -1 ; 323 324 // FIT_ID 325 *fitidRF = -1 ; 326 327 // OPACITY 328 *tauRF = 0.0 ; 329 330 // 331 // ITERATION: OBSERVATION_ID 332 // 333 TableIterator iter0( mstable_, "OBSERVATION_ID" ) ; 334 while( !iter0.pastEnd() ) { 335 //time0 = mathutil::gettimeofday_sec() ; 336 //os_ << "start 0th iteration: " << time0 << LogIO::POST ; 337 Table t0 = iter0.table() ; 338 Int obsId = asInt( "OBSERVATION_ID", 0, t0, tpoolr ) ; 339 if ( sdh.observer == "" ) { 340 sdh.observer = asString( "OBSERVER", obsId, obstab, tpoolr ) ; 341 } 342 if ( sdh.project == "" ) { 343 sdh.project = asString( "PROJECT", obsId, obstab, tpoolr ) ; 344 } 345 ROArrayMeasColumn<MEpoch> *tmpMeasCol = new ROArrayMeasColumn<MEpoch>( obstab, "TIME_RANGE" ) ; 346 MEpoch me = (*tmpMeasCol)( obsId )( IPosition(1,0) ) ; 347 delete tmpMeasCol ; 348 if ( sdh.utc == 0.0 ) { 349 sdh.utc = me.get( "d" ).getValue() ; 350 } 351 if ( telescopeName == "" ) { 352 telescopeName = asString( "TELESCOPE_NAME", obsId, obstab, tpoolr ) ; 353 } 354 Int nbeam = 0 ; 355 //time1 = mathutil::gettimeofday_sec() ; 356 //os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 357 // 358 // ITERATION: FEED1 359 // 360 TableIterator iter1( t0, "FEED1" ) ; 361 while( !iter1.pastEnd() ) { 362 //time0 = mathutil::gettimeofday_sec() ; 363 //os_ << "start 1st iteration: " << time0 << LogIO::POST ; 364 Table t1 = iter1.table() ; 365 // assume FEED1 == FEED2 366 Int feedId = asInt( "FEED1", 0, t1, tpoolr ) ; 367 nbeam++ ; 368 369 // BEAMNO 370 *beamRF = feedId ; 371 372 // FOCUS_ID 373 *focusidRF = 0 ; 374 375 //time1 = mathutil::gettimeofday_sec() ; 376 //os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 377 // 378 // ITERATION: FIELD_ID 379 // 380 TableIterator iter2( t1, "FIELD_ID" ) ; 381 while( !iter2.pastEnd() ) { 382 //time0 = mathutil::gettimeofday_sec() ; 383 //os_ << "start 2nd iteration: " << time0 << LogIO::POST ; 384 Table t2 = iter2.table() ; 385 Int fieldId = asInt( "FIELD_ID", 0, t2, tpoolr ) ; 386 Int srcId = asInt( "SOURCE_ID", fieldId, fieldtab, tpoolr ) ; 387 String fieldName = asString( "NAME", fieldId, fieldtab, tpoolr ) ; 388 fieldName += "__" + String::toString(fieldId) ; 389 ROArrayMeasColumn<MDirection> *delayDirCol = new ROArrayMeasColumn<MDirection>( fieldtab, "DELAY_DIR" ) ; 390 Vector<MDirection> delayDir = (*delayDirCol)( fieldId ) ; 391 delete delayDirCol ; 392 393 // FIELDNAME 394 *fieldnameRF = fieldName ; 395 396 397 //time1 = mathutil::gettimeofday_sec() ; 398 //os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 399 // 400 // ITERATION: DATA_DESC_ID 401 // 402 TableIterator iter3( t2, "DATA_DESC_ID" ) ; 403 while( !iter3.pastEnd() ) { 404 //time0 = mathutil::gettimeofday_sec() ; 405 //os_ << "start 3rd iteration: " << time0 << LogIO::POST ; 406 Table t3 = iter3.table() ; 407 Int ddId = asInt( "DATA_DESC_ID", 0, t3, tpoolr ) ; 408 Int polId = asInt( "POLARIZATION_ID", ddId, ddtab, tpoolr ) ; 409 Int spwId = asInt( "SPECTRAL_WINDOW_ID", ddId, ddtab, tpoolr ) ; 410 411 // IFNO 412 *ifnoRF = (uInt)spwId ; 413 414 // polarization information 415 Int npol = asInt( "NUM_CORR", polId, poltab, tpoolr ) ; 416 ROArrayColumn<Int> *roArrICol = new ROArrayColumn<Int>( poltab, "CORR_TYPE" ) ; 417 Vector<Int> corrtype = (*roArrICol)( polId ) ; 418 delete roArrICol ; 419 String srcName( "" ) ; 420 Vector<Double> srcPM( 2, 0.0 ) ; 421 Vector<Double> srcDir( 2, 0.0 ) ; 422 MDirection md ; 423 Vector<Double> restFreqs ; 424 Vector<String> transitionName ; 425 Vector<Double> sysVels ; 426 // os_ << "npol = " << npol << LogIO::POST ; 427 // os_ << "corrtype = " << corrtype << LogIO::POST ; 428 429 // source and molecular transition 430 sourceInfo( srcId, spwId, srcName, md, srcPM, restFreqs, transitionName, sysVels, tpoolr ) ; 431 // os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ; 432 433 // SRCNAME 434 *srcnameRF = srcName ; 435 436 // os_ << "srcName = " << srcName << LogIO::POST ; 437 438 // SRCPROPERMOTION 439 *srcpmRF = srcPM ; 440 441 //os_ << "srcPM = " << srcPM << LogIO::POST ; 442 443 // SRCDIRECTION 444 *srcdirRF = md.getAngle().getValue( "rad" ) ; 445 446 //os_ << "srcDir = " << srcDir << LogIO::POST ; 447 448 // SRCVELOCITY 449 Double sysVel = 0.0 ; 450 if ( !sysVels.empty() ) 451 sysVel = sysVels[0] ; 452 *sysvelRF = sysVel ; 453 454 // os_ << "sysVel = " << sysVel << LogIO::POST ; 455 456 uInt molId = table_->molecules().addEntry( restFreqs, transitionName, transitionName ) ; 457 458 // MOLECULE_ID 459 *molidRF = molId ; 460 461 // spectral setup 462 uInt freqId ; 463 Int nchan = asInt( "NUM_CHAN", spwId, spwtab, tpoolr ) ; 464 Bool iswvr = False ; 465 if ( nchan == 4 ) iswvr = True ; 466 sdh.nchan = max( sdh.nchan, nchan ) ; 467 map<Int,uInt>::iterator iter = ifmap.find( spwId ) ; 468 if ( iter == ifmap.end() ) { 469 ROScalarMeasColumn<MEpoch> *tmpMeasCol = new ROScalarMeasColumn<MEpoch>( t3, "TIME" ) ; 470 me = (*tmpMeasCol)( 0 ) ; 471 delete tmpMeasCol ; 472 Double refpix ; 473 Double refval ; 474 Double increment ; 475 spectralSetup( spwId, 476 me, 477 mp, 478 md, 479 refpix, 480 refval, 481 increment, 482 nchan, 483 sdh.freqref, 484 sdh.reffreq, 485 sdh.bandwidth, 486 tpoolr ) ; 487 freqId = table_->frequencies().addEntry( refpix, refval, increment ) ; 488 ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ; 489 //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ; 490 } 491 else { 492 freqId = iter->second ; 493 } 494 495 // FREQ_ID 496 *freqidRF = freqId ; 497 498 // for TSYS and TCAL 499 Vector<MEpoch> scTime ; 500 Vector<Double> scInterval ; 501 ROArrayColumn<Float> scTsysCol ; 502 MSSysCal caltabsel ; 503 if ( isSysCal_ ) { 504 caltabsel = caltab( caltab.col("ANTENNA_ID") == antenna_ && caltab.col("FEED_ID") == feedId && caltab.col("SPECTRAL_WINDOW_ID") == spwId ).sort("TIME") ; 505 ROScalarMeasColumn<MEpoch> scTimeCol( caltabsel, "TIME" ) ; 506 scTime.resize( caltabsel.nrow() ) ; 507 for ( uInt irow = 0 ; irow < caltabsel.nrow() ; irow++ ) 508 scTime[irow] = scTimeCol( irow ) ; 509 ROScalarColumn<Double> scIntervalCol( caltabsel, "INTERVAL" ) ; 510 scIntervalCol.getColumn( scInterval ) ; 511 if ( colTsys_ != "NONE" ) 512 scTsysCol.attach( caltabsel, colTsys_ ) ; 513 } 514 515 sdh.npol = max( sdh.npol, npol ) ; 516 if ( !iswvr && sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ; 517 518 //time1 = mathutil::gettimeofday_sec() ; 519 //os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 520 // 521 // ITERATION: SCAN_NUMBER 522 // 523 TableIterator iter4( t3, "SCAN_NUMBER" ) ; 524 while( !iter4.pastEnd() ) { 525 //time0 = mathutil::gettimeofday_sec() ; 526 //os_ << "start 4th iteration: " << time0 << LogIO::POST ; 527 Table t4 = iter4.table() ; 528 Int scanNum = asInt( "SCAN_NUMBER", 0, t4, tpoolr ) ; 529 530 // SCANNO 531 *scanRF = scanNum - 1 ; 532 533 uInt cycle = 0 ; 534 535 //time1 = mathutil::gettimeofday_sec() ; 536 //os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 537 // 538 // ITERATION: STATE_ID 539 // 540 TableIterator iter5( t4, "STATE_ID" ) ; 541 while( !iter5.pastEnd() ) { 542 //time0 = mathutil::gettimeofday_sec() ; 543 //os_ << "start 5th iteration: " << time0 << LogIO::POST ; 544 Table t5 = iter5.table() ; 545 Int stateId = asInt( "STATE_ID", 0, t5, tpoolr ) ; 546 String obstype = asString( "OBS_MODE", 0, stattab, tpoolr ) ; 547 if ( sdh.obstype == "" ) sdh.obstype = obstype ; 548 549 Int nrow = t5.nrow() ; 550 //time1 = mathutil::gettimeofday_sec() ; 551 //os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 552 553 Cube<Float> spArr ; 554 Cube<Bool> flArr ; 555 reshapeSpectraAndFlagtra( spArr, 556 flArr, 557 t5, 558 npol, 559 nchan, 560 nrow, 561 corrtype ) ; 562 if ( sdh.fluxunit == "" ) { 563 String colName = "FLOAT_DATA" ; 564 if ( isData_ ) colName = "DATA" ; 565 ROTableColumn dataCol( t5, colName ) ; 566 const TableRecord &dataColKeys = dataCol.keywordSet() ; 567 if ( dataColKeys.isDefined( "UNIT" ) ) 568 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ; 569 else if ( dataColKeys.isDefined( "QuantumUnits" ) ) 570 sdh.fluxunit = dataColKeys.asString( "QuantumUnits" ) ; 571 } 572 573 ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ; 574 Block<MEpoch> mTimeB( nrow ) ; 575 for ( Int irow = 0 ; irow < nrow ; irow++ ) 576 mTimeB[irow] = (*mTimeCol)( irow ) ; 577 Block<Int> sysCalIdx( nrow, -1 ) ; 578 if ( isSysCal_ ) { 579 getSysCalTime( scTime, scInterval, mTimeB, sysCalIdx ) ; 580 } 581 delete mTimeCol ; 582 Matrix<Float> defaulttsys( npol, 1, 1.0 ) ; 583 Int srcType = getSrcType( stateId, tpoolr ) ; 584 uInt diridx = 0 ; 585 uInt wid = 0 ; 586 Int pidx = 0 ; 587 Bool crossOK = False ; 588 Block<uInt> polnos( npol, 99 ) ; 589 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 590 Block<uInt> p = getPolNo( corrtype[ipol] ) ; 591 if ( p.size() > 1 ) { 592 if ( crossOK ) continue ; 593 else { 594 polnos[pidx] = p[0] ; 595 pidx++ ; 596 polnos[pidx] = p[1] ; 597 pidx++ ; 598 crossOK = True ; 599 } 600 } 601 else { 602 polnos[pidx] = p[0] ; 603 pidx++ ; 604 } 605 } 606 607 // SRCTYPE 608 *srctypeRF = srcType ; 609 610 for ( Int irow = 0 ; irow < nrow ; irow++ ) { 611 // CYCLENO 612 *cycleRF = cycle ; 613 614 // FLAGROW 615 *flrRF = (uInt)asBool( "FLAG_ROW", irow, t5, tpoolr ) ; 616 617 // SPECTRA, FLAG 618 Matrix<Float> sp = spArr.xyPlane( irow ) ; 619 Matrix<Bool> flb = flArr.xyPlane( irow ) ; 620 Matrix<uChar> fl( flb.shape() ) ; 621 convertArray( fl, flb ) ; 622 623 // TIME 624 *timeRF = mTimeB[irow].get("d").getValue() ; 625 626 // INTERVAL 627 *intervalRF = asDouble( "INTERVAL", irow, t5, tpoolr ) ; 628 629 // TSYS 630 Matrix<Float> tsys ; 631 if ( sysCalIdx[irow] != -1 && colTsys_ != "NONE" ) 632 tsys = scTsysCol( sysCalIdx[irow] ) ; 633 else 634 tsys = defaulttsys ; 635 636 // TCAL_ID 637 Block<uInt> tcalids( npol, 0 ) ; 638 if ( sysCalIdx[irow] != -1 && colTcal_ != "NONE" ) { 639 tcalids = getTcalId( feedId, spwId, scTime[sysCalIdx[irow]] ) ; 640 } 641 642 // WEATHER_ID 643 if ( isWeather_ ) { 644 wid = getWeatherId( wid, mTimeB[irow].get("s").getValue() ) ; 645 *widRF = mwIndex_[wid] ; 646 } 647 else { 648 *widRF = wid ; 649 } 650 651 652 // DIRECTION, AZEL, SCANRATE 653 Vector<Double> dir ; 654 Vector<Double> azel ; 655 Vector<Double> scanrate = defaultScanrate ; 656 String refString ; 657 if ( getPt_ ) 658 diridx = getDirection( diridx, dir, azel, scanrate, pt, pdcol, mTimeB[irow], mp ) ; 659 else 660 getSourceDirection( dir, azel, scanrate, mTimeB[irow], mp, delayDir ) ; 661 *dirRF = dir ; 662 *azRF = azel[0] ; 663 *elRF = azel[1] ; 664 *scrRF = scanrate ; 665 666 // Polarization dependent things 667 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 668 // POLNO 669 *polnoRF = polnos[ipol] ; 670 671 spRF.define( sp.row( ipol ) ) ; 672 ucarrRF.define( fl.row( ipol ) ) ; 673 tsysRF.define( tsys.row( ipol ) ) ; 674 *tcalidRF = tcalids[ipol] ; 675 676 // Commit row 677 stab.addRow() ; 678 row.put( stab.nrow()-1 ) ; 679 } 680 681 cycle++ ; 682 } 683 684 //time1 = mathutil::gettimeofday_sec() ; 685 //os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 686 687 iter5.next() ; 688 } 689 iter4.next() ; 690 } 691 iter3.next() ; 692 } 693 iter2.next() ; 694 } 695 iter1.next() ; 696 } 697 if ( sdh.nbeam < nbeam ) sdh.nbeam = nbeam ; 698 699 iter0.next() ; 700 } 701 702 703 delete tpoolr ; 704 705 706 // Table Keywords 707 sdh.nif = ifmap.size() ; 708 if ( ( telescopeName == "" ) || ( antennaName == telescopeName ) ) { 709 sdh.antennaname = antennaName ; 710 } 711 else { 712 sdh.antennaname = telescopeName + "//" + antennaName ; 713 } 714 if ( stationName != "" && stationName != antennaName ) { 715 sdh.antennaname += "@" + stationName ; 716 } 717 ROArrayColumn<Double> pdirCol( pointtab, "DIRECTION" ) ; 718 String dirref = pdirCol.keywordSet().asRecord("MEASINFO").asString("Ref") ; 719 if ( dirref == "AZELGEO" || dirref == "AZEL" ) { 720 dirref = "J2000" ; 721 } 722 sscanf( dirref.chars()+1, "%f", &sdh.equinox ) ; 723 sdh.epoch = "UTC" ; 724 if (sdh.freqref == "TOPO") { 725 sdh.freqref = "TOPOCENT"; 726 } else if (sdh.freqref == "GEO") { 727 sdh.freqref = "GEOCENTR"; 728 } else if (sdh.freqref == "BARY") { 729 sdh.freqref = "BARYCENT"; 730 } else if (sdh.freqref == "GALACTO") { 731 sdh.freqref = "GALACTOC"; 732 } else if (sdh.freqref == "LGROUP") { 733 sdh.freqref = "LOCALGRP"; 734 } else if (sdh.freqref == "CMB") { 735 sdh.freqref = "CMBDIPOL"; 736 } else if (sdh.freqref == "REST") { 737 sdh.freqref = "SOURCE"; 738 } 739 740 if ( sdh.fluxunit == "" || sdh.fluxunit == "CNTS" ) 741 sdh.fluxunit = "K" ; 742 table_->setHeader( sdh ) ; 1717 /*** 1718 * Start iteration using TableVisitor 1719 ***/ 1720 Table stab = table_->table() ; 1721 { 1722 static const char *cols[] = { 1723 "OBSERVATION_ID", "FEED1", "FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER", 1724 "STATE_ID", "TIME", 1725 NULL 1726 }; 1727 static const TypeManagerImpl<Int> tmInt; 1728 static const TypeManagerImpl<Double> tmDouble; 1729 static const TypeManager *const tms[] = { 1730 &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmDouble, NULL 1731 }; 1732 //double t0 = mathutil::gettimeofday_sec() ; 1733 MSFillerVisitor myVisitor(mstable_, *table_ ); 1734 //double t1 = mathutil::gettimeofday_sec() ; 1735 //cout << "MSFillerVisitor(): elapsed time " << t1-t0 << " sec" << endl ; 1736 myVisitor.setAntenna( antenna_ ) ; 1737 //myVisitor.setHeader( sdh ) ; 1738 if ( getPt_ ) { 1739 Table ptsel = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort( "TIME" ) ; 1740 myVisitor.setPointingTable( ptsel ) ; 1741 } 1742 if ( isWeather_ ) 1743 myVisitor.setWeatherTime( mwTime_, mwInterval_ ) ; 1744 if ( isSysCal_ ) 1745 myVisitor.setSysCalRecord( tcalrec_ ) ; 1746 1747 //double t2 = mathutil::gettimeofday_sec() ; 1748 traverseTable(mstable_, cols, tms, &myVisitor); 1749 //double t3 = mathutil::gettimeofday_sec() ; 1750 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ; 1751 1752 sdh = myVisitor.getHeader() ; 1753 } 1754 /*** 1755 * End iteration using TableVisitor 1756 ***/ 1757 1758 // set header 1759 //sdh = myVisitor.getHeader() ; 1760 //table_->setHeader( sdh ) ; 743 1761 744 1762 // save path to POINTING table … … 753 1771 754 1772 // for GBT 755 if ( antennaName.contains( "GBT" ) ) {1773 if ( sdh.antennaname.contains( "GBT" ) ) { 756 1774 String goTabName = datapath.absoluteName() + "/GBT_GO" ; 757 1775 stab.rwKeywordSet().define( "GBT_GO", goTabName ) ; … … 759 1777 760 1778 // for MS created from ASDM 761 //mstable_.keywordSet().print(cout) ;762 1779 const TableRecord &msKeys = mstable_.keywordSet() ; 763 1780 uInt nfields = msKeys.nfields() ; … … 782 1799 //tablesel_.unlock() ; 783 1800 mstable_.unlock() ; 784 }785 786 Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool )787 {788 //double startSec = mathutil::gettimeofday_sec() ;789 //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;790 791 MSState statetab = mstable_.state() ;792 String obsMode = asString( "OBS_MODE", stateId, statetab, tpool ) ;793 Bool sig = asBool( "SIG", stateId, statetab, tpool ) ;794 Bool ref = asBool( "REF", stateId, statetab, tpool ) ;795 Double cal = asDouble( "CAL", stateId, statetab, tpool ) ;796 //os_ << "OBS_MODE = " << obsMode << LogIO::POST ;797 798 // determine separator799 String sep = "" ;800 String tmpStr = obsMode.substr( 0, obsMode.find_first_of( "," ) ) ;801 //os_ << "tmpStr = " << tmpStr << LogIO::POST ;802 //if ( obsMode.find( ":" ) != String::npos ) {803 if ( tmpStr.find( ":" ) != String::npos ) {804 sep = ":" ;805 }806 //else if ( obsMode.find( "." ) != String::npos ) {807 else if ( tmpStr.find( "." ) != String::npos ) {808 sep = "." ;809 }810 else if ( tmpStr.find( "#" ) != String::npos ) {811 sep = "#" ;812 }813 //else if ( obsMode.find( "_" ) != String::npos ) {814 else if ( tmpStr.find( "_" ) != String::npos ) {815 sep = "_" ;816 }817 //os_ << "separator = " << sep << LogIO::POST ;818 819 // determine SRCTYPE820 Int srcType = SrcType::NOTYPE ;821 if ( sep == ":" ) {822 // sep == ":"823 //824 // GBT case825 //826 // obsMode1=Nod827 // NOD828 // obsMode1=OffOn829 // obsMode2=PSWITCHON: PSON830 // obsMode2=PSWITCHOFF: PSOFF831 // obsMode1=??832 // obsMode2=FSWITCH:833 // SIG=1: FSON834 // REF=1: FSOFF835 // Calibration scan if CAL != 0836 Int epos = obsMode.find_first_of( sep ) ;837 Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;838 String obsMode1 = obsMode.substr( 0, epos ) ;839 String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;840 if ( obsMode1 == "Nod" ) {841 srcType = SrcType::NOD ;842 }843 else if ( obsMode1 == "OffOn" ) {844 if ( obsMode2 == "PSWITCHON" ) srcType = SrcType::PSON ;845 if ( obsMode2 == "PSWITCHOFF" ) srcType = SrcType::PSOFF ;846 }847 else {848 if ( obsMode2 == "FSWITCH" ) {849 if ( sig ) srcType = SrcType::FSON ;850 if ( ref ) srcType = SrcType::FSOFF ;851 }852 }853 if ( cal > 0.0 ) {854 if ( srcType == SrcType::NOD )855 srcType = SrcType::NODCAL ;856 else if ( srcType == SrcType::PSON )857 srcType = SrcType::PONCAL ;858 else if ( srcType == SrcType::PSOFF )859 srcType = SrcType::POFFCAL ;860 else if ( srcType == SrcType::FSON )861 srcType = SrcType::FONCAL ;862 else if ( srcType == SrcType::FSOFF )863 srcType = SrcType::FOFFCAL ;864 else865 srcType = SrcType::CAL ;866 }867 }868 else if ( sep == "." || sep == "#" ) {869 // sep == "." or "#"870 //871 // ALMA & EVLA case (MS via ASDM) before3.1872 //873 // obsMode1=CALIBRATE_*874 // obsMode2=ON_SOURCE: PONCAL875 // obsMode2=OFF_SOURCE: POFFCAL876 // obsMode1=OBSERVE_TARGET877 // obsMode2=ON_SOURCE: PON878 // obsMode2=OFF_SOURCE: POFF879 string substr[2] ;880 int numSubstr = split( obsMode, substr, 2, "," ) ;881 //os_ << "numSubstr = " << numSubstr << LogIO::POST ;882 //for ( int i = 0 ; i < numSubstr ; i++ )883 //os_ << "substr[" << i << "] = " << substr[i] << LogIO::POST ;884 String obsType( substr[0] ) ;885 //os_ << "obsType = " << obsType << LogIO::POST ;886 Int epos = obsType.find_first_of( sep ) ;887 Int nextpos = obsType.find_first_of( sep, epos+1 ) ;888 String obsMode1 = obsType.substr( 0, epos ) ;889 String obsMode2 = obsType.substr( epos+1, nextpos-epos-1 ) ;890 //os_ << "obsMode1 = " << obsMode1 << LogIO::POST ;891 //os_ << "obsMode2 = " << obsMode2 << LogIO::POST ;892 if ( obsMode1.find( "CALIBRATE_" ) == 0 ) {893 if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PONCAL ;894 if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::POFFCAL ;895 }896 else if ( obsMode1 == "OBSERVE_TARGET" ) {897 if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PSON ;898 if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::PSOFF ;899 }900 }901 else if ( sep == "_" ) {902 // sep == "_"903 //904 // ALMA & EVLA case (MS via ASDM) after 3.2905 //906 // obsMode1=CALIBRATE_*907 // obsMode2=ON_SOURCE: PONCAL908 // obsMode2=OFF_SOURCE: POFFCAL909 // obsMode1=OBSERVE_TARGET910 // obsMode2=ON_SOURCE: PON911 // obsMode2=OFF_SOURCE: POFF912 string substr[2] ;913 int numSubstr = split( obsMode, substr, 2, "," ) ;914 //os_ << "numSubstr = " << numSubstr << LogIO::POST ;915 //for ( int i = 0 ; i < numSubstr ; i++ )916 //os_ << "substr[" << i << "] = " << substr[i] << LogIO::POST ;917 String obsType( substr[0] ) ;918 //os_ << "obsType = " << obsType << LogIO::POST ;919 string substr2[4] ;920 int numSubstr2 = split( obsType, substr2, 4, sep ) ;921 //Int epos = obsType.find_first_of( sep ) ;922 //Int nextpos = obsType.find_first_of( sep, epos+1 ) ;923 //String obsMode1 = obsType.substr( 0, epos ) ;924 //String obsMode2 = obsType.substr( epos+1, nextpos-epos-1 ) ;925 String obsMode1( substr2[0] ) ;926 String obsMode2( substr2[2] ) ;927 //os_ << "obsMode1 = " << obsMode1 << LogIO::POST ;928 //os_ << "obsMode2 = " << obsMode2 << LogIO::POST ;929 if ( obsMode1.find( "CALIBRATE" ) == 0 ) {930 if ( obsMode2 == "ON" ) srcType = SrcType::PONCAL ;931 if ( obsMode2 == "OFF" ) srcType = SrcType::POFFCAL ;932 }933 else if ( obsMode1 == "OBSERVE" ) {934 if ( obsMode2 == "ON" ) srcType = SrcType::PSON ;935 if ( obsMode2 == "OFF" ) srcType = SrcType::PSOFF ;936 }937 }938 else {939 if ( sig ) srcType = SrcType::SIG ;940 if ( ref ) srcType = SrcType::REF ;941 }942 943 //os_ << "srcType = " << srcType << LogIO::POST ;944 //double endSec = mathutil::gettimeofday_sec() ;945 //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;946 return srcType ;947 }948 949 //Vector<uInt> MSFiller::getPolNo( Int corrType )950 Block<uInt> MSFiller::getPolNo( Int corrType )951 {952 //double startSec = mathutil::gettimeofday_sec() ;953 //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;954 Block<uInt> polno( 1 ) ;955 956 if ( corrType == Stokes::I || corrType == Stokes::RR || corrType == Stokes::XX ) {957 polno = 0 ;958 }959 else if ( corrType == Stokes::Q || corrType == Stokes::LL || corrType == Stokes::YY ) {960 polno = 1 ;961 }962 else if ( corrType == Stokes::U ) {963 polno = 2 ;964 }965 else if ( corrType == Stokes::V ) {966 polno = 3 ;967 }968 else if ( corrType == Stokes::RL || corrType == Stokes::XY || corrType == Stokes::LR || corrType == Stokes::RL ) {969 polno.resize( 2 ) ;970 polno[0] = 2 ;971 polno[1] = 3 ;972 }973 else if ( corrType == Stokes::Plinear ) {974 polno[0] = 1 ;975 }976 else if ( corrType == Stokes::Pangle ) {977 polno[0] = 2 ;978 }979 else {980 polno = 99 ;981 }982 //os_ << "polno = " << polno << LogIO::POST ;983 //double endSec = mathutil::gettimeofday_sec() ;984 //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;985 986 return polno ;987 }988 989 String MSFiller::getPolType( Int corrType )990 {991 //double startSec = mathutil::gettimeofday_sec() ;992 //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;993 String poltype = "" ;994 995 if ( corrType == Stokes::I || corrType == Stokes::Q || corrType == Stokes::U || corrType == Stokes::V )996 poltype = "stokes" ;997 else if ( corrType == Stokes::XX || corrType == Stokes::YY || corrType == Stokes::XY || corrType == Stokes::YX )998 poltype = "linear" ;999 else if ( corrType == Stokes::RR || corrType == Stokes::LL || corrType == Stokes::RL || corrType == Stokes::LR )1000 poltype = "circular" ;1001 else if ( corrType == Stokes::Plinear || corrType == Stokes::Pangle )1002 poltype = "linpol" ;1003 1004 //double endSec = mathutil::gettimeofday_sec() ;1005 //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1006 return poltype ;1007 1801 } 1008 1802 … … 1164 1958 } 1165 1959 1166 void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr)1960 void MSFiller::fillTcal() 1167 1961 { 1168 1962 //double startSec = mathutil::gettimeofday_sec() ; … … 1189 1983 } 1190 1984 1191 Table sctab = mstable_.sysCal() ;1985 Table &sctab = mstable_.sysCal() ; 1192 1986 if ( sctab.nrow() == 0 ) { 1193 1987 os_ << "No SYSCAL rows" << LogIO::POST ; 1194 1988 return ; 1195 1989 } 1196 Table sctabsel( sctab( sctab.col("ANTENNA_ID") == antenna_ ) ) ; 1197 if ( sctabsel.nrow() == 0 ) { 1990 ROScalarColumn<Int> antCol( sctab, "ANTENNA_ID" ) ; 1991 Vector<Int> ant = antCol.getColumn() ; 1992 if ( allNE( ant, antenna_ ) ) { 1198 1993 os_ << "No SYSCAL rows" << LogIO::POST ; 1199 1994 return ; 1200 1995 } 1201 ROArrayColumn<Float> *tmpTcalCol = new ROArrayColumn<Float>( sctabsel, colTcal_ ) ; 1202 // return if any rows without Tcal value exists 1996 ROTableColumn tcalCol( sctab, colTcal_ ) ; 1203 1997 Bool notDefined = False ; 1204 for ( uInt irow = 0 ; irow < sctab sel.nrow() ; irow++ ) {1205 if ( !tmpTcalCol->isDefined( irow ) ) {1998 for ( uInt irow = 0 ; irow < sctab.nrow() ; irow++ ) { 1999 if ( ant[irow] == antenna_ && !tcalCol.isDefined( irow ) ) { 1206 2000 notDefined = True ; 1207 2001 break ; … … 1210 2004 if ( notDefined ) { 1211 2005 os_ << "No TCAL value" << LogIO::POST ; 1212 delete tmpTcalCol ;1213 2006 table_->tcal().table().addRow(1,True) ; 1214 2007 Vector<Float> defaultTcal( 1, 1.0 ) ; … … 1217 2010 return ; 1218 2011 } 1219 uInt npol = tmpTcalCol->shape( 0 )(0) ; 1220 delete tmpTcalCol ; 1221 //os_ << "fillTcal(): npol = " << npol << LogIO::POST ; 2012 2013 static const char *cols[] = { 2014 "ANTENNA_ID", "FEED_ID", "SPECTRAL_WINDOW_ID", "TIME", 2015 NULL 2016 }; 2017 static const TypeManagerImpl<Int> tmInt; 2018 static const TypeManagerImpl<Double> tmDouble; 2019 static const TypeManager *const tms[] = { 2020 &tmInt, &tmInt, &tmInt, &tmDouble, NULL 2021 }; 1222 2022 Table tab = table_->tcal().table() ; 1223 ArrayColumn<Float> tcalCol( tab, "TCAL" ) ; 1224 uInt oldnr = 0 ; 1225 uInt newnr = 0 ; 1226 TableRow row( tab ) ; 1227 TableRecord &trec = row.record() ; 1228 RecordFieldPtr<uInt> idRF( trec, "ID" ) ; 1229 RecordFieldPtr<String> timeRF( trec, "TIME" ) ; 1230 RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ; 1231 TableIterator iter0( sctabsel, "FEED_ID" ) ; 1232 while( !iter0.pastEnd() ) { 1233 Table t0 = iter0.table() ; 1234 Int feedId = asInt( "FEED_ID", 0, t0, tpoolr ) ; 1235 TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ; 1236 while( !iter1.pastEnd() ) { 1237 Table t1 = iter1.table() ; 1238 Int spwId = asInt( "SPECTRAL_WINDOW_ID", 0, t1, tpoolr ) ; 1239 tmpTcalCol = new ROArrayColumn<Float>( t1, colTcal_ ) ; 1240 ROScalarQuantColumn<Double> scTimeCol( t1, "TIME" ) ; 1241 Vector<uInt> idminmax( 2, oldnr ) ; 1242 for ( uInt irow = 0 ; irow < t1.nrow() ; irow++ ) { 1243 String sTime = MVTime( scTimeCol(irow) ).string( MVTime::YMD ) ; 1244 *timeRF = sTime ; 1245 uInt idx = oldnr ; 1246 Matrix<Float> subtcal = (*tmpTcalCol)( irow ) ; 1247 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) { 1248 *idRF = idx++ ; 1249 //*tcalRF = subtcal.row( ipol ) ; 1250 tcalRF.define( subtcal.row( ipol ) ) ; 1251 1252 // commit row 1253 tab.addRow() ; 1254 row.put( tab.nrow()-1 ) ; 1255 1256 newnr++ ; 1257 } 1258 idminmax[0] = oldnr ; 1259 idminmax[1] = newnr - 1 ; 1260 oldnr = newnr ; 1261 1262 String key = keyTcal( feedId, spwId, sTime ) ; 1263 tcalrec_.define( key, idminmax ) ; 1264 } 1265 delete tmpTcalCol ; 1266 iter1++ ; 1267 } 1268 iter0++ ; 1269 } 2023 TcalVisitor visitor( sctab, tab, tcalrec_, antenna_ ) ; 2024 visitor.setTcalColumn( colTcal_ ) ; 2025 2026 traverseTable(sctab, cols, tms, &visitor); 1270 2027 1271 2028 //tcalrec_.print( std::cout ) ; 1272 2029 //double endSec = mathutil::gettimeofday_sec() ; 1273 2030 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1274 }1275 1276 uInt MSFiller::getWeatherId( uInt idx, Double wtime )1277 {1278 //double startSec = mathutil::gettimeofday_sec() ;1279 //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;1280 uInt nrow = mwTime_.size() ;1281 if ( nrow < 2 )1282 return 0 ;1283 uInt wid = nrow ;1284 if ( idx == 0 ) {1285 wid = 0 ;1286 Double tStart = mwTime_[wid]-0.5*mwInterval_[wid] ;1287 if ( wtime < tStart )1288 return wid ;1289 }1290 for ( uInt i = idx ; i < nrow-1 ; i++ ) {1291 Double tStart = mwTime_[i]-0.5*mwInterval_[i] ;1292 // use of INTERVAL column is problematic1293 // since there are "blank" time of weather monitoring1294 //Double tEnd = tStart + mwInterval_[i] ;1295 Double tEnd = mwTime_[i+1]-0.5*mwInterval_[i+1] ;1296 //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;1297 if ( wtime >= tStart && wtime <= tEnd ) {1298 wid = i ;1299 break ;1300 }1301 }1302 if ( wid == nrow ) {1303 uInt i = nrow - 1 ;1304 Double tStart = mwTime_[i-1]+0.5*mwInterval_[i-1] ;1305 Double tEnd = mwTime_[i]+0.5*mwInterval_[i] ;1306 //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;1307 if ( wtime >= tStart && wtime <= tEnd )1308 wid = i-1 ;1309 else1310 wid = i ;1311 }1312 1313 //if ( wid == nrow )1314 //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ;1315 1316 //double endSec = mathutil::gettimeofday_sec() ;1317 //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1318 return wid ;1319 }1320 1321 void MSFiller::getSysCalTime( Vector<MEpoch> &scTime, Vector<Double> &scInterval, Block<MEpoch> &tcol, Block<Int> &tidx )1322 {1323 //double startSec = mathutil::gettimeofday_sec() ;1324 //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;1325 1326 if ( !isSysCal_ )1327 return ;1328 1329 uInt nrow = tidx.nelements() ;1330 if ( scTime.nelements() == 0 )1331 return ;1332 else if ( scTime.nelements() == 1 ) {1333 tidx[0] = 0 ;1334 return ;1335 }1336 uInt scnrow = scTime.nelements() ;1337 uInt idx = 0 ;1338 const Double half = 0.5e0 ;1339 // execute binary search1340 idx = binarySearch( scTime, tcol[0].get( "s" ).getValue() ) ;1341 if ( idx != 0 )1342 idx -= 1 ;1343 for ( uInt i = 0 ; i < nrow ; i++ ) {1344 Double t = tcol[i].get( "s" ).getValue() ;1345 Double tsc = scTime[0].get( "s" ).getValue() ;1346 if ( t < tsc ) {1347 tidx[i] = 0 ;1348 continue ;1349 }1350 for ( uInt j = idx ; j < scnrow-1 ; j++ ) {1351 Double tsc1 = scTime[j].get( "s" ).getValue() ;1352 Double dt1 = scInterval[j] ;1353 Double tsc2 = scTime[j+1].get( "s" ).getValue() ;1354 Double dt2 = scInterval[j+1] ;1355 if ( t > tsc1-half*dt1 && t <= tsc2-half*dt2 ) {1356 tidx[i] = j ;1357 idx = j ;1358 break ;1359 }1360 }1361 if ( tidx[i] == -1 ) {1362 // Double tsc = scTime[scnrow-1].get( "s" ).getValue() ;1363 // Double dt = scInterval[scnrow-1] ;1364 // if ( t <= tsc+0.5*dt ) {1365 // tidx[i] = scnrow-1 ;1366 // }1367 tidx[i] = scnrow-1 ;1368 }1369 }1370 //double endSec = mathutil::gettimeofday_sec() ;1371 //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;1372 return ;1373 }1374 1375 Block<uInt> MSFiller::getTcalId( Int fid, Int spwid, MEpoch &t )1376 {1377 //double startSec = mathutil::gettimeofday_sec() ;1378 //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;1379 //if ( table_->tcal().table().nrow() == 0 ) {1380 if ( !isSysCal_ ) {1381 os_ << "No TCAL rows" << LogIO::POST ;1382 Block<uInt> tcalids( 4, 0 ) ;1383 return tcalids ;1384 }1385 //String sctime = MVTime( Quantum<Double>(t,"s") ).string(MVTime::YMD) ;1386 String sctime = MVTime( t.getValue() ).string(MVTime::YMD) ;1387 String key = keyTcal( fid, spwid, sctime ) ;1388 if ( !tcalrec_.isDefined( key ) ) {1389 os_ << "No TCAL rows" << LogIO::POST ;1390 Block<uInt> tcalids( 4, 0 ) ;1391 return tcalids ;1392 }1393 Vector<uInt> ids = tcalrec_.asArrayuInt( key ) ;1394 uInt npol = ids[1] - ids[0] + 1 ;1395 Block<uInt> tcalids( npol ) ;1396 tcalids[0] = ids[0] ;1397 tcalids[1] = ids[1] ;1398 for ( uInt ipol = 2 ; ipol < npol ; ipol++ )1399 tcalids[ipol] = ids[0] + ipol - 1 ;1400 1401 //double endSec = mathutil::gettimeofday_sec() ;1402 //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1403 return tcalids ;1404 }1405 1406 uInt MSFiller::getDirection( uInt idx,1407 Vector<Double> &dir,1408 Vector<Double> &srate,1409 String &ref,1410 Vector<Double> &tcol,1411 ROArrayColumn<Double> &dcol,1412 Double t )1413 {1414 //double startSec = mathutil::gettimeofday_sec() ;1415 //os_ << "start MSFiller::getDirection1() startSec=" << startSec << LogIO::POST ;1416 //double time0 = mathutil::gettimeofday_sec() ;1417 //os_ << "start getDirection 1st stage startSec=" << time0 << LogIO::POST ;1418 // assume that cols is sorted by TIME1419 Bool doInterp = False ;1420 uInt nrow = dcol.nrow() ;1421 if ( nrow == 0 )1422 return 0 ;1423 if ( idx == 0 ) {1424 uInt nrowb = 1000 ;1425 if ( nrow > nrowb ) {1426 uInt nblock = nrow / nrowb + 1 ;1427 for ( uInt iblock = 0 ; iblock < nblock ; iblock++ ) {1428 uInt high = min( nrowb, nrow-iblock*nrowb ) ;1429 1430 if ( tcol( high-1 ) < t ) {1431 idx = iblock * nrowb ;1432 continue ;1433 }1434 1435 Slice slice( iblock*nrowb, nrowb ) ;1436 Vector<Double> tarr = tcol( slice ) ;1437 1438 uInt bidx = binarySearch( tarr, t ) ;1439 1440 idx = iblock * nrowb + bidx ;1441 break ;1442 }1443 }1444 else {1445 idx = binarySearch( tcol, t ) ;1446 }1447 }1448 //double time1 = mathutil::gettimeofday_sec() ;1449 //os_ << "end getDirection 1st stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;1450 // ensure that tcol(idx) < t1451 //os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ;1452 //time0 = mathutil::gettimeofday_sec() ;1453 //os_ << "start getDirection 2nd stage startSec=" << time0 << LogIO::POST ;1454 //while( tcol( idx ) * factor > t && idx > 0 )1455 while( tcol[idx] > t && idx > 0 )1456 idx-- ;1457 //os_ << "idx = " << idx << LogIO::POST ;1458 1459 // index search1460 for ( uInt i = idx ; i < nrow ; i++ ) {1461 Double tref = tcol[i] ;1462 if ( tref == t ) {1463 idx = i ;1464 break ;1465 }1466 else if ( tref > t ) {1467 if ( i == 0 ) {1468 idx = i ;1469 }1470 else {1471 idx = i-1 ;1472 doInterp = True ;1473 }1474 break ;1475 }1476 else {1477 idx = nrow - 1 ;1478 }1479 }1480 //time1 = mathutil::gettimeofday_sec() ;1481 //os_ << "end getDirection 2nd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;1482 //os_ << "searched idx = " << idx << LogIO::POST ;1483 1484 //time0 = mathutil::gettimeofday_sec() ;1485 //os_ << "start getDirection 3rd stage startSec=" << time0 << LogIO::POST ;1486 //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;1487 //IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;1488 IPosition ip( dcol(idx).shape().nelements(), 0 ) ;1489 //os_ << "ip = " << ip << LogIO::POST ;1490 //ref = dmcol(idx)(ip).getRefString() ;1491 TableRecord trec = dcol.keywordSet() ;1492 Record rec = trec.asRecord( "MEASINFO" ) ;1493 ref = rec.asString( "Ref" ) ;1494 //os_ << "ref = " << ref << LogIO::POST ;1495 if ( doInterp ) {1496 //os_ << "do interpolation" << LogIO::POST ;1497 //os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ;1498 Double tref0 = tcol[idx] ;1499 Double tref1 = tcol[idx+1] ;1500 Matrix<Double> mdir0 = dcol( idx ) ;1501 Matrix<Double> mdir1 = dcol( idx+1 ) ;1502 Vector<Double> dir0 = mdir0.column( 0 ) ;1503 //os_ << "dir0 = " << dir0 << LogIO::POST ;1504 Vector<Double> dir1 = mdir1.column( 0 ) ;1505 //os_ << "dir1 = " << dir1 << LogIO::POST ;1506 Double dt0 = t - tref0 ;1507 Double dt1 = tref1 - t ;1508 dir.reference( (dt0*dir1+dt1*dir0)/(dt0+dt1) ) ;1509 if ( mdir0.ncolumn() > 1 ) {1510 if ( dt0 >= dt1 )1511 srate.reference( mdir0.column( 1 ) ) ;1512 else1513 srate.reference( mdir1.column( 1 ) ) ;1514 }1515 //os_ << "dir = " << dir << LogIO::POST ;1516 }1517 else {1518 //os_ << "no interpolation" << LogIO::POST ;1519 Matrix<Double> mdir0 = dcol( idx ) ;1520 dir.reference( mdir0.column( 0 ) ) ;1521 if ( mdir0.ncolumn() > 1 )1522 srate.reference( mdir0.column( 1 ) ) ;1523 }1524 1525 //time1 = mathutil::gettimeofday_sec() ;1526 //os_ << "end getDirection 3rd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;1527 //double endSec = mathutil::gettimeofday_sec() ;1528 //os_ << "end MSFiller::getDirection1() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1529 return idx ;1530 }1531 1532 String MSFiller::keyTcal( Int feedid, Int spwid, String stime )1533 {1534 String sfeed = "FEED" + String::toString( feedid ) ;1535 String sspw = "SPW" + String::toString( spwid ) ;1536 return sfeed+":"+sspw+":"+stime ;1537 }1538 1539 uInt MSFiller::binarySearch( Vector<MEpoch> &timeList, Double target )1540 {1541 Int low = 0 ;1542 Int high = timeList.nelements() ;1543 uInt idx = 0 ;1544 1545 while ( low <= high ) {1546 idx = (Int)( 0.5 * ( low + high ) ) ;1547 Double t = timeList[idx].get( "s" ).getValue() ;1548 if ( t < target )1549 low = idx + 1 ;1550 else if ( t > target )1551 high = idx - 1 ;1552 else {1553 return idx ;1554 }1555 }1556 1557 idx = max( 0, min( low, high ) ) ;1558 1559 return idx ;1560 }1561 1562 uInt MSFiller::binarySearch( Vector<Double> &timeList, Double target )1563 {1564 Int low = 0 ;1565 Int high = timeList.nelements() ;1566 uInt idx = 0 ;1567 1568 while ( low <= high ) {1569 idx = (Int)( 0.5 * ( low + high ) ) ;1570 Double t = timeList[idx] ;1571 if ( t < target )1572 low = idx + 1 ;1573 else if ( t > target )1574 high = idx - 1 ;1575 else {1576 return idx ;1577 }1578 }1579 1580 idx = max( 0, min( low, high ) ) ;1581 1582 return idx ;1583 2031 } 1584 2032 … … 1603 2051 } 1604 2052 1605 void MSFiller::reshapeSpectraAndFlagtra( Cube<Float> &sp,1606 Cube<Bool> &fl,1607 Table &tab,1608 Int &npol,1609 Int &nchan,1610 Int &nrow,1611 Vector<Int> &corrtype )1612 {1613 //double startSec = mathutil::gettimeofday_sec() ;1614 //os_ << "start MSFiller::reshapeSpectraAndFlagtra() startSec=" << startSec << LogIO::POST ;1615 if ( isFloatData_ ) {1616 ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;1617 ROArrayColumn<Float> mFloatDataCol( tab, "FLOAT_DATA" ) ;1618 mFloatDataCol.getColumn( sp ) ;1619 mFlagCol.getColumn( fl ) ;1620 }1621 else if ( isData_ ) {1622 sp.resize( npol, nchan, nrow ) ;1623 fl.resize( npol, nchan, nrow ) ;1624 ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;1625 ROArrayColumn<Complex> mDataCol( tab, "DATA" ) ;1626 if ( npol < 3 ) {1627 Cube<Float> tmp = ComplexToReal( mDataCol.getColumn() ) ;1628 IPosition start( 3, 0, 0, 0 ) ;1629 IPosition end( 3, 2*npol-1, nchan-1, nrow-1 ) ;1630 IPosition inc( 3, 2, 1, 1 ) ;1631 sp = tmp( start, end, inc ) ;1632 fl = mFlagCol.getColumn() ;1633 }1634 else {1635 for ( Int irow = 0 ; irow < nrow ; irow++ ) {1636 Bool crossOK = False ;1637 Matrix<Complex> mSp = mDataCol( irow ) ;1638 Matrix<Bool> mFl = mFlagCol( irow ) ;1639 Matrix<Float> spxy = sp.xyPlane( irow ) ;1640 Matrix<Bool> flxy = fl.xyPlane( irow ) ;1641 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {1642 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX1643 || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {1644 if ( !crossOK ) {1645 Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ;1646 IPosition start( 1, 0 ) ;1647 IPosition end( 1, 2*nchan-1 ) ;1648 IPosition inc( 1, 2 ) ;1649 spxy.row( ipol ) = tmp( start, end, inc ) ;1650 flxy.row( ipol ) = mFl.row( ipol ) ;1651 start = IPosition( 1, 1 ) ;1652 spxy.row( ipol+1 ) = tmp( start, end, inc ) ;1653 flxy.row( ipol+1 ) = mFl.row( ipol ) ;1654 if ( corrtype[ipol] == Stokes::YX || corrtype[ipol] == Stokes::LR ) {1655 spxy.row( ipol+1 ) = spxy.row( ipol+1 ) * (Float)-1.0 ;1656 }1657 crossOK = True ;1658 }1659 }1660 else {1661 Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ;1662 IPosition start( 1, 0 ) ;1663 IPosition end( 1, 2*nchan-1 ) ;1664 IPosition inc( 1, 2 ) ;1665 spxy.row( ipol ) = tmp( start, end, inc ) ;1666 flxy.row( ipol ) = mFl.row( ipol ) ;1667 }1668 }1669 }1670 }1671 }1672 //double endSec = mathutil::gettimeofday_sec() ;1673 //os_ << "end MSFiller::reshapeSpectraAndFlagtra() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1674 }1675 1676 uInt MSFiller::getDirection( uInt idx,1677 Vector<Double> &dir,1678 Vector<Double> &azel,1679 Vector<Double> &srate,1680 Vector<Double> &ptcol,1681 ROArrayColumn<Double> &pdcol,1682 MEpoch &t,1683 MPosition &antpos )1684 {1685 //double startSec = mathutil::gettimeofday_sec() ;1686 //os_ << "start MSFiller::getDirection2() startSec=" << startSec << LogIO::POST ;1687 String refString ;1688 MDirection::Types dirType ;1689 uInt diridx = getDirection( idx, dir, srate, refString, ptcol, pdcol, t.get("s").getValue() ) ;1690 MDirection::getType( dirType, refString ) ;1691 MeasFrame mf( t, antpos ) ;1692 if ( refString == "J2000" ) {1693 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;1694 azel = toazel( dir ).getAngle("rad").getValue() ;1695 }1696 else if ( refString(0,4) == "AZEL" ) {1697 azel = dir.copy() ;1698 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;1699 dir = toj2000( dir ).getAngle("rad").getValue() ;1700 }1701 else {1702 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;1703 azel = toazel( dir ).getAngle("rad").getValue() ;1704 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;1705 dir = toj2000( dir ).getAngle("rad").getValue() ;1706 }1707 if ( srate.size() == 0 ) {1708 srate.resize( 2 ) ;1709 srate = 0.0 ;1710 }1711 //double endSec = mathutil::gettimeofday_sec() ;1712 //os_ << "end MSFiller::getDirection2() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1713 return diridx ;1714 }1715 1716 void MSFiller::getSourceDirection( Vector<Double> &dir,1717 Vector<Double> &azel,1718 Vector<Double> &srate,1719 MEpoch &t,1720 MPosition &antpos,1721 Vector<MDirection> &srcdir )1722 {1723 //double startSec = mathutil::gettimeofday_sec() ;1724 //os_ << "start MSFiller::getSourceDirection() startSec=" << startSec << LogIO::POST ;1725 Vector<Double> defaultDir = srcdir[0].getAngle( "rad" ).getValue() ;1726 if ( srcdir.nelements() > 1 )1727 srate = srcdir[1].getAngle( "rad" ).getValue() ;1728 String ref = srcdir[0].getRefString() ;1729 MDirection::Types dirType ;1730 MDirection::getType( dirType, ref ) ;1731 MeasFrame mf( t, antpos ) ;1732 if ( ref != "J2000" ) {1733 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;1734 dir = toj2000( defaultDir ).getAngle("rad").getValue() ;1735 }1736 else1737 dir = defaultDir ;1738 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;1739 azel = toazel( defaultDir ).getAngle("rad").getValue() ;1740 //double endSec = mathutil::gettimeofday_sec() ;1741 //os_ << "end MSFiller::getSourceDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1742 }1743 1744 2053 void MSFiller::initHeader( STHeader &header ) 1745 2054 { … … 1752 2061 header.obstype = "" ; 1753 2062 header.antennaname = "" ; 1754 header.antennaposition.resize( 0) ;2063 header.antennaposition.resize( 3 ) ; 1755 2064 header.equinox = 0.0 ; 1756 2065 header.freqref = "" ; … … 1763 2072 } 1764 2073 1765 String MSFiller::asString( String name, 1766 uInt idx, 1767 Table tab, 1768 boost::object_pool<ROTableColumn> *pool ) 1769 { 1770 ROTableColumn *col = pool->construct( tab, name ) ; 1771 String v = col->asString( idx ) ; 1772 pool->destroy( col ) ; 1773 return v ; 1774 } 1775 1776 Bool MSFiller::asBool( String name, 1777 uInt idx, 1778 Table &tab, 1779 boost::object_pool<ROTableColumn> *pool ) 1780 { 1781 ROTableColumn *col = pool->construct( tab, name ) ; 1782 Bool v = col->asBool( idx ) ; 1783 pool->destroy( col ) ; 1784 return v ; 1785 } 1786 1787 uInt MSFiller::asuInt( String name, 1788 uInt idx, 1789 Table &tab, 1790 boost::object_pool<ROTableColumn> *pool ) 1791 { 1792 ROTableColumn *col = pool->construct( tab, name ) ; 1793 uInt v = col->asuInt( idx ) ; 1794 pool->destroy( col ) ; 1795 return v ; 1796 } 1797 1798 Int MSFiller::asInt( String name, 1799 uInt idx, 1800 Table &tab, 1801 boost::object_pool<ROTableColumn> *pool ) 1802 { 1803 ROTableColumn *col = pool->construct( tab, name ) ; 1804 Int v = col->asInt( idx ) ; 1805 pool->destroy( col ) ; 1806 return v ; 1807 } 1808 1809 Float MSFiller::asFloat( String name, 1810 uInt idx, 1811 Table &tab, 1812 boost::object_pool<ROTableColumn> *pool ) 1813 { 1814 ROTableColumn *col = pool->construct( tab, name ) ; 1815 Float v = col->asfloat( idx ) ; 1816 pool->destroy( col ) ; 1817 return v ; 1818 } 1819 1820 Double MSFiller::asDouble( String name, 1821 uInt idx, 1822 Table &tab, 1823 boost::object_pool<ROTableColumn> *pool ) 1824 { 1825 ROTableColumn *col = pool->construct( tab, name ) ; 1826 Double v = col->asdouble( idx ) ; 1827 pool->destroy( col ) ; 1828 return v ; 1829 } 1830 1831 void MSFiller::sourceInfo( Int sourceId, 1832 Int spwId, 1833 String &name, 1834 MDirection &direction, 1835 Vector<casa::Double> &properMotion, 1836 Vector<casa::Double> &restFreqs, 1837 Vector<casa::String> &transitions, 1838 Vector<casa::Double> &sysVels, 1839 boost::object_pool<ROTableColumn> *tpoolr ) 1840 { 1841 //double startSec = mathutil::gettimeofday_sec() ; 1842 //os_ << "start MSFiller::sourceInfo() startSec=" << startSec << LogIO::POST ; 1843 1844 MSSource srctab = mstable_.source() ; 1845 MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ; 1846 if ( srctabSel.nrow() == 0 ) { 1847 srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ; 1848 } 1849 Int numLines = 0 ; 1850 if ( srctabSel.nrow() > 0 ) { 1851 // source name 1852 name = asString( "NAME", 0, srctabSel, tpoolr ) ; 1853 1854 // source proper motion 1855 ROArrayColumn<Double> roArrDCol( srctabSel, "PROPER_MOTION" ) ; 1856 properMotion = roArrDCol( 0 ) ; 1857 1858 // source direction as MDirection object 1859 ROScalarMeasColumn<MDirection> tmpMeasCol( srctabSel, "DIRECTION" ) ; 1860 direction = tmpMeasCol( 0 ) ; 1861 1862 // number of lines 1863 numLines = asInt( "NUM_LINES", 0, srctabSel, tpoolr ) ; 1864 } 1865 else { 1866 name = "" ; 1867 properMotion = Vector<Double>( 2, 0.0 ) ; 1868 direction = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ; 1869 } 1870 1871 restFreqs.resize( numLines ) ; 1872 transitions.resize( numLines ) ; 1873 sysVels.resize( numLines ) ; 1874 if ( numLines > 0 ) { 1875 if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) { 1876 ROArrayQuantColumn<Double> quantArrCol( srctabSel, "REST_FREQUENCY" ) ; 1877 Array< Quantum<Double> > qRestFreqs = quantArrCol( 0 ) ; 1878 for ( int i = 0 ; i < numLines ; i++ ) { 1879 restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ; 1880 } 1881 } 1882 //os_ << "restFreqs = " << restFreqs << LogIO::POST ; 1883 if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) { 1884 ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ; 1885 if ( transitionCol.isDefined( 0 ) ) 1886 transitions = transitionCol( 0 ) ; 1887 //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ; 1888 } 1889 if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) { 1890 ROArrayColumn<Double> roArrDCol( srctabSel, "SYSVEL" ) ; 1891 sysVels = roArrDCol( 0 ) ; 1892 } 1893 } 1894 1895 //double endSec = mathutil::gettimeofday_sec() ; 1896 //os_ << "end MSFiller::sourceInfo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1897 } 1898 1899 void MSFiller::spectralSetup( Int spwId, 1900 MEpoch &me, 1901 MPosition &mp, 1902 MDirection &md, 1903 Double &refpix, 1904 Double &refval, 1905 Double &increment, 1906 Int &nchan, 1907 String &freqref, 1908 Double &reffreq, 1909 Double &bandwidth, 1910 boost::object_pool<ROTableColumn> *tpoolr ) 1911 { 1912 //double startSec = mathutil::gettimeofday_sec() ; 1913 //os_ << "start MSFiller::spectralSetup() startSec=" << startSec << LogIO::POST ; 1914 1915 MSSpectralWindow spwtab = mstable_.spectralWindow() ; 1916 MeasFrame mf( me, mp, md ) ; 1917 MFrequency::Types freqRef = MFrequency::castType( (uInt)asInt( "MEAS_FREQ_REF", spwId, spwtab, tpoolr ) ) ; 1918 Bool even = False ; 1919 if ( (nchan/2)*2 == nchan ) even = True ; 1920 ROScalarQuantColumn<Double> tmpQuantCol( spwtab, "TOTAL_BANDWIDTH" ) ; 1921 Double totbw = tmpQuantCol( spwId ).getValue( "Hz" ) ; 1922 if ( nchan != 4 ) 1923 bandwidth = max( bandwidth, totbw ) ; 1924 if ( freqref == "" && nchan != 4) 1925 //sdh.freqref = MFrequency::showType( freqRef ) ; 1926 freqref = "LSRK" ; 1927 if ( reffreq == -1.0 && nchan != 4 ) { 1928 tmpQuantCol.attach( spwtab, "REF_FREQUENCY" ) ; 1929 Quantum<Double> qreffreq = tmpQuantCol( spwId ) ; 1930 if ( freqRef == MFrequency::LSRK ) { 1931 reffreq = qreffreq.getValue("Hz") ; 1932 } 1933 else { 1934 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 1935 reffreq = tolsr( qreffreq ).get("Hz").getValue() ; 1936 } 1937 } 1938 Int refchan = nchan / 2 ; 1939 IPosition refip( 1, refchan ) ; 1940 refpix = 0.5*(nchan-1) ; 1941 refval = 0.0 ; 1942 ROArrayQuantColumn<Double> sharedQDArrCol( spwtab, "CHAN_WIDTH" ) ; 1943 ROTableColumn netSidebandCol( spwtab, "NET_SIDEBAND" ) ; 1944 Int netSideband = netSidebandCol.asInt( spwId ) ; 1945 increment = sharedQDArrCol( spwId )( refip ).getValue( "Hz" ) ; 1946 // os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ; 1947 sharedQDArrCol.attach( spwtab, "CHAN_FREQ" ) ; 1948 Vector< Quantum<Double> > chanFreqs = sharedQDArrCol( spwId ) ; 1949 if ( ( nchan > 1 && 1950 chanFreqs[0].getValue("Hz") > chanFreqs[1].getValue("Hz") ) 1951 || ( nchan == 1 && netSideband == 1 ) ) 1952 increment *= -1.0 ; 1953 if ( freqRef == MFrequency::LSRK ) { 1954 if ( even ) { 1955 IPosition refip0( 1, refchan-1 ) ; 1956 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 1957 Double refval1 = chanFreqs(refip).getValue("Hz") ; 1958 refval = 0.5 * ( refval0 + refval1 ) ; 1959 } 1960 else { 1961 refval = chanFreqs(refip).getValue("Hz") ; 1962 } 1963 } 1964 else { 1965 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 1966 if ( even ) { 1967 IPosition refip0( 1, refchan-1 ) ; 1968 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 1969 Double refval1 = chanFreqs(refip).getValue("Hz") ; 1970 refval = 0.5 * ( refval0 + refval1 ) ; 1971 refval = tolsr( refval ).get("Hz").getValue() ; 1972 } 1973 else { 1974 refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ; 1975 } 1976 } 1977 1978 //double endSec = mathutil::gettimeofday_sec() ; 1979 //os_ << "end MSFiller::spectralSetup() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1980 } 1981 1982 } ; 1983 2074 }; -
trunk/src/MSFiller.h
r2289 r2291 16 16 // STL 17 17 #include <string> 18 19 // Boost20 #include <boost/pool/object_pool.hpp>21 18 22 19 // AIPS++ … … 27 24 #include <casa/Arrays/Cube.h> 28 25 #include <casa/Logging/LogIO.h> 29 26 #include <casa/Containers/RecordField.h> 30 27 #include <casa/Containers/Record.h> 31 28 #include <casa/Containers/Block.h> 29 #include <casa/Quanta/MVTime.h> 32 30 33 31 #include <ms/MeasurementSets/MeasurementSet.h> 34 32 #include <ms/MeasurementSets/MSPointing.h> 35 33 34 #include <tables/Tables/ScalarColumn.h> 35 #include <tables/Tables/ArrayColumn.h> 36 #include <tables/Tables/TableRow.h> 37 38 #include <measures/TableMeasures/ScalarMeasColumn.h> 39 #include <measures/TableMeasures/ArrayMeasColumn.h> 40 #include <measures/TableMeasures/ScalarQuantColumn.h> 41 #include <measures/TableMeasures/ArrayQuantColumn.h> 42 43 #include "TableTraverse.h" 36 44 #include "Scantable.h" 45 #include "MathUtils.h" 46 47 using namespace casa; 37 48 38 49 namespace asap 39 50 { 51 52 class MSFillerUtils { 53 protected: 54 template<class T> void getScalar( const String &name, 55 const uInt &idx, 56 const Table &tab, 57 T &val ) 58 { 59 ROScalarColumn<T> col( tab, name ) ; 60 val = col( idx ) ; 61 } 62 template<class T> void getArray( const String &name, 63 const uInt &idx, 64 const Table &tab, 65 Array<T> &val ) 66 { 67 ROArrayColumn<T> col( tab, name ) ; 68 val = col( idx ) ; 69 } 70 template<class T> void getScalarMeas( const String &name, 71 const uInt &idx, 72 const Table &tab, 73 T &val ) 74 { 75 ROScalarMeasColumn<T> measCol( tab, name ) ; 76 val = measCol( idx ) ; 77 } 78 template<class T> void getArrayMeas( const String &name, 79 const uInt &idx, 80 const Table &tab, 81 Array<T> &val ) 82 { 83 ROArrayMeasColumn<T> measCol( tab, name ) ; 84 val = measCol( idx ) ; 85 } 86 template<class T> void getScalarQuant( const String &name, 87 const uInt &idx, 88 const Table &tab, 89 Quantum<T> &val ) 90 { 91 ROScalarQuantColumn<T> quantCol( tab, name ) ; 92 val = quantCol( idx ) ; 93 } 94 template<class T> void getArrayQuant( const String &name, 95 const uInt &idx, 96 const Table &tab, 97 Array< Quantum<T> > &val ) 98 { 99 ROArrayQuantColumn<T> quantCol( tab, name ) ; 100 val = quantCol( idx ) ; 101 } 102 // template<class T> void putField( const String &name, 103 // TableRecord &rec, 104 // T &val ) 105 // { 106 // RecordFieldPtr<T> rf( rec, name ) ; 107 // *rf = val ; 108 // } 109 // template<class T> void defineField( const String &name, 110 // TableRecord &rec, 111 // T &val ) 112 // { 113 // RecordFieldPtr<T> rf( rec, name ) ; 114 // rf.define( val ) ; 115 // } 116 template<class T> T interp( Double x0, Double x1, Double x, T y0, T y1 ) 117 { 118 Double dx0 = x - x0 ; 119 Double dx1 = x1 - x ; 120 return ( y0 * dx1 + y1 * dx0 ) / ( x1 - x0 ) ; 121 } 122 String keyTcal( const Int &feedid, const Int &spwid, const Double &time ) 123 { 124 String stime = MVTime( Quantity(time,Unit("s")) ).string( MVTime::YMD ) ; 125 String sfeed = "FEED" + String::toString( feedid ) ; 126 String sspw = "SPW" + String::toString( spwid ) ; 127 return sfeed+":"+sspw+":"+stime ; 128 } 129 String keyTcal( const Int &feedid, const Int &spwid, const String &stime ) 130 { 131 String sfeed = "FEED" + String::toString( feedid ) ; 132 String sspw = "SPW" + String::toString( spwid ) ; 133 return sfeed+":"+sspw+":"+stime ; 134 } 135 }; 40 136 41 137 class MSFiller 42 138 { 43 139 public: 44 explicit MSFiller( casa::CountedPtr<Scantable> stable) ;140 explicit MSFiller(CountedPtr<Scantable> stable) ; 45 141 virtual ~MSFiller() ; 46 142 47 virtual bool open(const std::string& filename, const casa::Record& rec) ;143 virtual bool open(const std::string& filename, const Record& rec) ; 48 144 virtual void fill() ; 49 145 virtual void close() ; … … 65 161 //void fillHistory() ; 66 162 //void fillFit() ; 67 void fillTcal( boost::object_pool<casa::ROTableColumn> *poolr ) ; 68 69 // get SRCTYPE from STATE_ID 70 casa::Int getSrcType( casa::Int stateId, boost::object_pool<casa::ROTableColumn> *pool ) ; 71 72 // get POLNO from CORR_TYPE 73 casa::Block<casa::uInt> getPolNo( casa::Int corrType ) ; 74 75 // get poltype from CORR_TYPE 76 casa::String getPolType( casa::Int corrType ) ; 77 78 // get WEATHER_ID 79 casa::uInt getWeatherId( casa::uInt idx, casa::Double wtime ) ; 80 81 // get time stamp in SYSCAL table 82 // assume that tab is selected by ANTENNA_ID, FEED_ID, SPECTRAL_WINDOW_ID 83 // and sorted by TIME 84 void getSysCalTime( casa::Vector<casa::MEpoch> &scTimeIn, casa::Vector<casa::Double> &scInterval, casa::Block<casa::MEpoch> &tcol, casa::Block<casa::Int> &tidx ) ; 85 86 // get TCAL_ID 87 casa::Block<casa::uInt> getTcalId( casa::Int feedId, casa::Int spwId, casa::MEpoch &t ) ; 88 89 // get direction for DIRECTION, AZIMUTH, and ELEVATION columns 90 casa::uInt getDirection( casa::uInt idx, 91 casa::Vector<casa::Double> &dir, 92 casa::Vector<casa::Double> &srate, 93 casa::String &ref, 94 casa::Vector<casa::Double> &ptcol, 95 casa::ROArrayColumn<casa::Double> &pdcol, 96 casa::Double t ) ; 97 casa::uInt getDirection( casa::uInt idx, 98 casa::Vector<casa::Double> &dir, 99 casa::Vector<casa::Double> &azel, 100 casa::Vector<casa::Double> &srate, 101 casa::Vector<casa::Double> &ptcol, 102 casa::ROArrayColumn<casa::Double> &pdcol, 103 casa::MEpoch &t, casa::MPosition &antpos ) ; 104 void getSourceDirection( casa::Vector<casa::Double> &dir, 105 casa::Vector<casa::Double> &azel, 106 casa::Vector<casa::Double> &srate, 107 casa::MEpoch &t, 108 casa::MPosition &antpos, 109 casa::Vector<casa::MDirection> &srcdir ) ; 163 void fillTcal() ; 110 164 111 165 // create key for TCAL table 112 casa::String keyTcal( casa::Int feedid, casa::Int spwid, casa::String stime ) ; 113 114 // binary search 115 casa::uInt binarySearch( casa::Vector<casa::MEpoch> &timeList, casa::Double target ) ; 116 casa::uInt binarySearch( casa::Vector<casa::Double> &timeList, casa::Double target ) ; 117 118 // tool for HPC 119 // double gettimeofday_sec() ; 166 String keyTcal( Int feedid, Int spwid, String stime ) ; 120 167 121 168 // get frequency frame 122 169 std::string getFrame() ; 123 170 124 // reshape SPECTRA and FLAGTRA125 void reshapeSpectraAndFlagtra( casa::Cube<casa::Float> &sp,126 casa::Cube<casa::Bool> &fl,127 casa::Table &tab,128 casa::Int &npol,129 casa::Int &nchan,130 casa::Int &nrow,131 casa::Vector<casa::Int> &corrtype ) ;132 133 171 // initialize header 134 172 void initHeader( STHeader &header ) ; 135 173 136 // get value from table column using object pool 137 casa::String asString( casa::String name, 138 casa::uInt idx, 139 casa::Table tab, 140 boost::object_pool<casa::ROTableColumn> *pool ) ; 141 casa::Bool asBool( casa::String name, 142 casa::uInt idx, 143 casa::Table &tab, 144 boost::object_pool<casa::ROTableColumn> *pool ) ; 145 casa::Int asInt( casa::String name, 146 casa::uInt idx, 147 casa::Table &tab, 148 boost::object_pool<casa::ROTableColumn> *pool ) ; 149 casa::uInt asuInt( casa::String name, 150 casa::uInt idx, 151 casa::Table &tab, 152 boost::object_pool<casa::ROTableColumn> *pool ) ; 153 casa::Float asFloat( casa::String name, 154 casa::uInt idx, 155 casa::Table &tab, 156 boost::object_pool<casa::ROTableColumn> *pool ) ; 157 casa::Double asDouble( casa::String name, 158 casa::uInt idx, 159 casa::Table &tab, 160 boost::object_pool<casa::ROTableColumn> *pool ) ; 161 162 void sourceInfo( casa::Int sourceId, 163 casa::Int spwId, 164 casa::String &name, 165 casa::MDirection &direction, 166 casa::Vector<casa::Double> &properMotion, 167 casa::Vector<casa::Double> &restFreqs, 168 casa::Vector<casa::String> &transitions, 169 casa::Vector<casa::Double> &sysVels, 170 boost::object_pool<casa::ROTableColumn> *pool ) ; 171 172 void spectralSetup( casa::Int spwId, 173 casa::MEpoch &me, 174 casa::MPosition &mp, 175 casa::MDirection &md, 176 casa::Double &refpix, 177 casa::Double &refval, 178 casa::Double &increment, 179 casa::Int &nchan, 180 casa::String &freqref, 181 casa::Double &reffreq, 182 casa::Double &bandwidth, 183 boost::object_pool<casa::ROTableColumn> *pool ) ; 184 185 casa::CountedPtr<Scantable> table_ ; 186 casa::MeasurementSet mstable_ ; 187 casa::String tablename_ ; 188 casa::Int antenna_ ; 189 casa::String antennaStr_ ; 190 casa::Bool getPt_ ; 191 192 casa::Bool isFloatData_ ; 193 casa::Bool isData_ ; 194 195 casa::Bool isDoppler_ ; 196 casa::Bool isFlagCmd_ ; 197 casa::Bool isFreqOffset_ ; 198 casa::Bool isHistory_ ; 199 casa::Bool isProcessor_ ; 200 casa::Bool isSysCal_ ; 201 casa::Bool isWeather_ ; 202 203 casa::String colTsys_ ; 204 casa::String colTcal_ ; 205 206 casa::LogIO os_ ; 207 208 casa::Vector<casa::Double> mwTime_ ; 209 casa::Vector<casa::Double> mwInterval_ ; 210 casa::Vector<casa::uInt> mwIndex_ ; 174 CountedPtr<Scantable> table_ ; 175 MeasurementSet mstable_ ; 176 String tablename_ ; 177 Int antenna_ ; 178 String antennaStr_ ; 179 Bool getPt_ ; 180 181 Bool isFloatData_ ; 182 Bool isData_ ; 183 184 Bool isDoppler_ ; 185 Bool isFlagCmd_ ; 186 Bool isFreqOffset_ ; 187 Bool isHistory_ ; 188 Bool isProcessor_ ; 189 Bool isSysCal_ ; 190 Bool isWeather_ ; 191 192 String colTsys_ ; 193 String colTcal_ ; 194 195 LogIO os_ ; 196 197 Vector<Double> mwTime_ ; 198 Vector<Double> mwInterval_ ; 199 Vector<uInt> mwIndex_ ; 211 200 212 201 // Record for TCAL_ID … … 214 203 // "SPW1": Vector<uInt> 215 204 // ... 216 casa::Record tcalrec_ ; 217 218 //casa::ROTableColumn *scCol_ ; 205 Record tcalrec_ ; 206 //map< String,Vector<uInt> > tcalrec_ ; 219 207 }; 220 208 -
trunk/src/MSWriter.cpp
r2259 r2291 11 11 // 12 12 // 13 #include <assert.h> 13 14 14 15 #include <casa/OS/File.h> … … 17 18 #include <casa/OS/SymLink.h> 18 19 #include <casa/BasicSL/String.h> 19 #include <casa/Containers/RecordField.h>20 20 #include <casa/Arrays/Cube.h> 21 21 … … 39 39 #include "STTcal.h" 40 40 #include "MathUtils.h" 41 42 // #include <ctime> 43 // #include <sys/time.h> 44 41 #include "TableTraverse.h" 45 42 46 43 using namespace casa ; … … 48 45 49 46 namespace asap { 50 // double MSWriter::gettimeofday_sec() 51 // { 52 // struct timeval tv ; 53 // gettimeofday( &tv, NULL ) ; 54 // return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ; 55 // } 47 48 class PolarizedComponentHolder { 49 public: 50 PolarizedComponentHolder() 51 : nchan(0), 52 maxnpol(4) 53 { 54 reset() ; 55 } 56 PolarizedComponentHolder( uInt n ) 57 : nchan(n), 58 maxnpol(4) 59 { 60 reset() ; 61 } 62 63 void reset() 64 { 65 npol = 0 ; 66 data.clear() ; 67 flag.clear() ; 68 flagrow = False ; 69 polnos.resize() ; 70 } 71 72 void accumulate( uInt id, Vector<Float> sp, Vector<Bool> fl, Bool flr ) 73 { 74 map< uInt,Vector<Float> >::iterator itr = data.find( id ) ; 75 if ( id < maxnpol && itr == data.end() ) { 76 addPol( id ) ; 77 accumulateData( id, sp ) ; 78 accumulateFlag( id, fl ) ; 79 accumulateFlagRow( flr ) ; 80 npol++ ; 81 } 82 } 83 84 void setNchan( uInt n ) { nchan = n ; } 85 uInt nPol() { return npol ; } 86 uInt nChan() { return nchan ; } 87 Vector<uInt> polNos() { return polnos ; } 88 Vector<Float> getWeight() { return Vector<Float>( npol, 1.0 ) ; } 89 Vector<Float> getSigma() { return Vector<Float>( npol, 1.0 ) ; } 90 Bool getFlagRow() { return flagrow ; } 91 Cube<Bool> getFlagCategory() { return Cube<Bool>( npol, nchan, 1, False ) ; } 92 Matrix<Float> getData() 93 { 94 Matrix<Float> v( npol, nchan ) ; 95 for ( map< uInt,Vector<Float> >::iterator i = data.begin() ; i != data.end() ; i++ ) { 96 v.row( i->first ) = i->second ; 97 } 98 return v ; 99 } 100 Matrix<Bool> getFlag() 101 { 102 Matrix<Bool> v( npol, nchan ) ; 103 for ( map< uInt,Vector<Bool> >::iterator i = flag.begin() ; i != flag.end() ; i++ ) { 104 v.row( i->first ) = i->second ; 105 } 106 return v ; 107 } 108 Matrix<Complex> getComplexData() 109 { 110 Matrix<Complex> v( npol, nchan ) ; 111 Matrix<Float> dummy( 2, nchan, 0.0 ) ; 112 map< uInt,Vector<Float> >::iterator itr0 = data.find( 0 ) ; 113 map< uInt,Vector<Float> >::iterator itr1 = data.find( 1 ) ; 114 if ( itr0 != data.end() ) { 115 dummy.row( 0 ) = itr0->second ; 116 v.row( 0 ) = RealToComplex( dummy ) ; 117 } 118 if ( itr1 != data.end() ) { 119 dummy.row( 0 ) = itr1->second ; 120 v.row( npol-1 ) = RealToComplex( dummy ) ; 121 } 122 itr0 = data.find( 2 ) ; 123 itr1 = data.find( 3 ) ; 124 if ( itr0 != data.end() && itr1 != data.end() ) { 125 dummy.row( 0 ) = itr0->second ; 126 dummy.row( 1 ) = itr1->second ; 127 v.row( 1 ) = RealToComplex( dummy ) ; 128 v.row( 2 ) = conj( v.row( 1 ) ) ; 129 } 130 return v ; 131 } 132 133 Matrix<Bool> getComplexFlag() 134 { 135 Matrix<Bool> tmp = getFlag() ; 136 Matrix<Bool> v( npol, nchan ) ; 137 v.row( 0 ) = tmp.row( 0 ) ; 138 if ( npol == 2 ) { 139 v.row( npol-1 ) = tmp.row( 1 ) ; 140 } 141 else if ( npol > 2 ) { 142 v.row( npol-1 ) = tmp.row( 1 ) ; 143 v.row( 1 ) = tmp.row( 2 ) || tmp.row( 3 ) ; 144 v.row( 2 ) = v.row( 1 ) ; 145 } 146 return v ; 147 } 148 149 private: 150 void accumulateData( uInt &id, Vector<Float> &v ) 151 { 152 data.insert( pair< uInt,Vector<Float> >( id, v ) ) ; 153 } 154 void accumulateFlag( uInt &id, Vector<Bool> &v ) 155 { 156 flag.insert( pair< uInt,Vector<Bool> >( id, v ) ) ; 157 } 158 void accumulateFlagRow( Bool &v ) 159 { 160 flagrow |= v ; 161 } 162 void addPol( uInt id ) 163 { 164 uInt i = polnos.nelements() ; 165 polnos.resize( i+1, True ) ; 166 polnos[i] = id ; 167 } 168 169 uInt nchan; 170 const uInt maxnpol; 171 uInt npol; 172 Vector<uInt> polnos; 173 174 map< uInt,Vector<Float> > data; 175 map< uInt,Vector<Bool> > flag; 176 Bool flagrow; 177 }; 178 179 class BaseMSWriterVisitor: public TableVisitor { 180 const String *lastFieldName; 181 uInt lastRecordNo; 182 uInt lastBeamNo, lastScanNo, lastIfNo; 183 Int lastSrcType; 184 uInt lastCycleNo; 185 Double lastTime; 186 Int lastPolNo; 187 protected: 188 const Table &table; 189 uInt count; 190 public: 191 BaseMSWriterVisitor(const Table &table) 192 : table(table) 193 { 194 static const String dummy; 195 lastFieldName = &dummy; 196 count = 0; 197 } 198 199 virtual void enterFieldName(const uInt recordNo, const String &columnValue) { 200 } 201 virtual void leaveFieldName(const uInt recordNo, const String &columnValue) { 202 } 203 virtual void enterBeamNo(const uInt recordNo, uInt columnValue) { } 204 virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) { } 205 virtual void enterScanNo(const uInt recordNo, uInt columnValue) { } 206 virtual void leaveScanNo(const uInt recordNo, uInt columnValue) { } 207 virtual void enterIfNo(const uInt recordNo, uInt columnValue) { } 208 virtual void leaveIfNo(const uInt recordNo, uInt columnValue) { } 209 virtual void enterSrcType(const uInt recordNo, Int columnValue) { } 210 virtual void leaveSrcType(const uInt recordNo, Int columnValue) { } 211 virtual void enterCycleNo(const uInt recordNo, uInt columnValue) { } 212 virtual void leaveCycleNo(const uInt recordNo, uInt columnValue) { } 213 virtual void enterTime(const uInt recordNo, Double columnValue) { } 214 virtual void leaveTime(const uInt recordNo, Double columnValue) { } 215 virtual void enterPolNo(const uInt recordNo, uInt columnValue) { } 216 virtual void leavePolNo(const uInt recordNo, uInt columnValue) { } 217 218 virtual Bool visitRecord(const uInt recordNo, 219 const String &fieldName, 220 const uInt beamNo, 221 const uInt scanNo, 222 const uInt ifNo, 223 const Int srcType, 224 const uInt cycleNo, 225 const Double time, 226 const uInt polNo) { return True ;} 227 228 virtual Bool visit(Bool isFirst, const uInt recordNo, 229 const uInt nCols, void const *const colValues[]) { 230 const String *fieldName = NULL; 231 uInt beamNo, scanNo, ifNo; 232 Int srcType; 233 uInt cycleNo; 234 Double time; 235 Int polNo; 236 { // prologue 237 uInt i = 0; 238 { 239 const String *col = (const String*)colValues[i++]; 240 fieldName = &col[recordNo]; 241 } 242 { 243 const uInt *col = (const uInt *)colValues[i++]; 244 beamNo = col[recordNo]; 245 } 246 { 247 const uInt *col = (const uInt *)colValues[i++]; 248 scanNo = col[recordNo]; 249 } 250 { 251 const uInt *col = (const uInt *)colValues[i++]; 252 ifNo = col[recordNo]; 253 } 254 { 255 const Int *col = (const Int *)colValues[i++]; 256 srcType = col[recordNo]; 257 } 258 { 259 const uInt *col = (const uInt *)colValues[i++]; 260 cycleNo = col[recordNo]; 261 } 262 { 263 const Double *col = (const Double *)colValues[i++]; 264 time = col[recordNo]; 265 } 266 { 267 const Int *col = (const Int *)colValues[i++]; 268 polNo = col[recordNo]; 269 } 270 assert(nCols == i); 271 } 272 273 if (isFirst) { 274 enterFieldName(recordNo, *fieldName); 275 enterBeamNo(recordNo, beamNo); 276 enterScanNo(recordNo, scanNo); 277 enterIfNo(recordNo, ifNo); 278 enterSrcType(recordNo, srcType); 279 enterCycleNo(recordNo, cycleNo); 280 enterTime(recordNo, time); 281 enterPolNo(recordNo, polNo); 282 } else { 283 if (lastFieldName->compare(*fieldName) != 0) { 284 leavePolNo(lastRecordNo, lastPolNo); 285 leaveTime(lastRecordNo, lastTime); 286 leaveCycleNo(lastRecordNo, lastCycleNo); 287 leaveSrcType(lastRecordNo, lastSrcType); 288 leaveIfNo(lastRecordNo, lastIfNo); 289 leaveScanNo(lastRecordNo, lastScanNo); 290 leaveBeamNo(lastRecordNo, lastBeamNo); 291 leaveFieldName(lastRecordNo, *lastFieldName); 292 293 enterFieldName(recordNo, *fieldName); 294 enterBeamNo(recordNo, beamNo); 295 enterScanNo(recordNo, scanNo); 296 enterIfNo(recordNo, ifNo); 297 enterSrcType(recordNo, srcType); 298 enterCycleNo(recordNo, cycleNo); 299 enterTime(recordNo, time); 300 enterPolNo(recordNo, polNo); 301 } else if (lastBeamNo != beamNo) { 302 leavePolNo(lastRecordNo, lastPolNo); 303 leaveTime(lastRecordNo, lastTime); 304 leaveCycleNo(lastRecordNo, lastCycleNo); 305 leaveSrcType(lastRecordNo, lastSrcType); 306 leaveIfNo(lastRecordNo, lastIfNo); 307 leaveScanNo(lastRecordNo, lastScanNo); 308 leaveBeamNo(lastRecordNo, lastBeamNo); 309 310 enterBeamNo(recordNo, beamNo); 311 enterScanNo(recordNo, scanNo); 312 enterIfNo(recordNo, ifNo); 313 enterSrcType(recordNo, srcType); 314 enterCycleNo(recordNo, cycleNo); 315 enterTime(recordNo, time); 316 enterPolNo(recordNo, polNo); 317 } else if (lastScanNo != scanNo) { 318 leavePolNo(lastRecordNo, lastPolNo); 319 leaveTime(lastRecordNo, lastTime); 320 leaveCycleNo(lastRecordNo, lastCycleNo); 321 leaveSrcType(lastRecordNo, lastSrcType); 322 leaveIfNo(lastRecordNo, lastIfNo); 323 leaveScanNo(lastRecordNo, lastScanNo); 324 325 enterScanNo(recordNo, scanNo); 326 enterIfNo(recordNo, ifNo); 327 enterSrcType(recordNo, srcType); 328 enterCycleNo(recordNo, cycleNo); 329 enterTime(recordNo, time); 330 enterPolNo(recordNo, polNo); 331 } else if (lastIfNo != ifNo) { 332 leavePolNo(lastRecordNo, lastPolNo); 333 leaveTime(lastRecordNo, lastTime); 334 leaveCycleNo(lastRecordNo, lastCycleNo); 335 leaveSrcType(lastRecordNo, lastSrcType); 336 leaveIfNo(lastRecordNo, lastIfNo); 337 338 enterIfNo(recordNo, ifNo); 339 enterSrcType(recordNo, srcType); 340 enterCycleNo(recordNo, cycleNo); 341 enterTime(recordNo, time); 342 enterPolNo(recordNo, polNo); 343 } else if (lastSrcType != srcType) { 344 leavePolNo(lastRecordNo, lastPolNo); 345 leaveTime(lastRecordNo, lastTime); 346 leaveCycleNo(lastRecordNo, lastCycleNo); 347 leaveSrcType(lastRecordNo, lastSrcType); 348 349 enterSrcType(recordNo, srcType); 350 enterCycleNo(recordNo, cycleNo); 351 enterTime(recordNo, time); 352 enterPolNo(recordNo, polNo); 353 } else if (lastCycleNo != cycleNo) { 354 leavePolNo(lastRecordNo, lastPolNo); 355 leaveTime(lastRecordNo, lastTime); 356 leaveCycleNo(lastRecordNo, lastCycleNo); 357 358 enterCycleNo(recordNo, cycleNo); 359 enterTime(recordNo, time); 360 enterPolNo(recordNo, polNo); 361 } else if (lastTime != time) { 362 leavePolNo(lastRecordNo, lastPolNo); 363 leaveTime(lastRecordNo, lastTime); 364 365 enterTime(recordNo, time); 366 enterPolNo(recordNo, polNo); 367 } else if (lastPolNo != polNo) { 368 leavePolNo(lastRecordNo, lastPolNo); 369 enterPolNo(recordNo, polNo); 370 } 371 } 372 count++; 373 Bool result = visitRecord(recordNo, *fieldName, beamNo, scanNo, ifNo, srcType, 374 cycleNo, time, polNo); 375 376 { // epilogue 377 lastRecordNo = recordNo; 378 379 lastFieldName = fieldName; 380 lastBeamNo = beamNo; 381 lastScanNo = scanNo; 382 lastIfNo = ifNo; 383 lastSrcType = srcType; 384 lastCycleNo = cycleNo; 385 lastTime = time; 386 lastPolNo = polNo; 387 } 388 return result ; 389 } 390 391 virtual void finish() { 392 if (count > 0) { 393 leavePolNo(lastRecordNo, lastPolNo); 394 leaveTime(lastRecordNo, lastTime); 395 leaveCycleNo(lastRecordNo, lastCycleNo); 396 leaveSrcType(lastRecordNo, lastSrcType); 397 leaveIfNo(lastRecordNo, lastIfNo); 398 leaveScanNo(lastRecordNo, lastScanNo); 399 leaveBeamNo(lastRecordNo, lastBeamNo); 400 leaveFieldName(lastRecordNo, *lastFieldName); 401 } 402 } 403 }; 404 405 class MSWriterVisitor: public BaseMSWriterVisitor, public MSWriterUtils { 406 public: 407 MSWriterVisitor(const Table &table, Table &mstable) 408 : BaseMSWriterVisitor(table), 409 ms(mstable) 410 { 411 rowidx = 0 ; 412 fieldName = "" ; 413 defaultFieldId = 0 ; 414 spwId = -1 ; 415 subscan = 1 ; 416 ptName = "" ; 417 srcId = 0 ; 418 419 row = TableRow( ms ) ; 420 421 holder.reset() ; 422 423 makePolMap() ; 424 initFrequencies() ; 425 initTcal() ; 426 427 // 428 // add rows to MS 429 // 430 uInt addrow = table.nrow() ; 431 ms.addRow( addrow ) ; 432 433 // attach to Scantable columns 434 spectraCol.attach( table, "SPECTRA" ) ; 435 flagtraCol.attach( table, "FLAGTRA" ) ; 436 flagRowCol.attach( table, "FLAGROW" ) ; 437 tcalIdCol.attach( table, "TCAL_ID" ) ; 438 intervalCol.attach( table, "INTERVAL" ) ; 439 directionCol.attach( table, "DIRECTION" ) ; 440 scanRateCol.attach( table, "SCANRATE" ) ; 441 timeCol.attach( table, "TIME" ) ; 442 freqIdCol.attach( table, "FREQ_ID" ) ; 443 sourceNameCol.attach( table, "SRCNAME" ) ; 444 sourceDirectionCol.attach( table, "SRCDIRECTION" ) ; 445 fieldNameCol.attach( table, "FIELDNAME" ) ; 446 447 // MS subtables 448 attachSubtables() ; 449 450 // attach to MS columns 451 attachMain() ; 452 attachPointing() ; 453 } 454 455 virtual void enterFieldName(const uInt recordNo, const String &columnValue) { 456 //printf("%u: FieldName: %s\n", recordNo, columnValue.c_str()); 457 fieldName = fieldNameCol.asString( recordNo ) ; 458 String::size_type pos = fieldName.find( "__" ) ; 459 if ( pos != String::npos ) { 460 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ; 461 fieldName = fieldName.substr( 0, pos ) ; 462 } 463 else { 464 fieldId = defaultFieldId ; 465 defaultFieldId++ ; 466 } 467 Double tSec = timeCol.asdouble( recordNo ) * 86400.0 ; 468 Vector<Double> srcDir = sourceDirectionCol( recordNo ) ; 469 Vector<Double> srate = scanRateCol( recordNo ) ; 470 String srcName = sourceNameCol.asString( recordNo ) ; 471 472 addField( fieldId, fieldName, srcName, srcDir, srate, tSec ) ; 473 474 // put value 475 *fieldIdRF = fieldId ; 476 } 477 virtual void leaveFieldName(const uInt recordNo, const String &columnValue) { 478 } 479 virtual void enterBeamNo(const uInt recordNo, uInt columnValue) { 480 //printf("%u: BeamNo: %u\n", recordNo, columnValue); 481 482 feedId = (Int)columnValue ; 483 484 // put value 485 *feed1RF = feedId ; 486 *feed2RF = feedId ; 487 } 488 virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) { 489 } 490 virtual void enterScanNo(const uInt recordNo, uInt columnValue) { 491 //printf("%u: ScanNo: %u\n", recordNo, columnValue); 492 493 // put value 494 // SCAN_NUMBER is 0-based in Scantable while 1-based in MS 495 *scanNumberRF = (Int)columnValue + 1 ; 496 } 497 virtual void leaveScanNo(const uInt recordNo, uInt columnValue) { 498 subscan = 1 ; 499 } 500 virtual void enterIfNo(const uInt recordNo, uInt columnValue) { 501 //printf("%u: IfNo: %u\n", recordNo, columnValue); 502 503 spwId = (Int)columnValue ; 504 uInt freqId = freqIdCol.asuInt( recordNo ) ; 505 506 Vector<Float> sp = spectraCol( recordNo ) ; 507 uInt nchan = sp.nelements() ; 508 holder.setNchan( nchan ) ; 509 510 addSpectralWindow( spwId, freqId ) ; 511 512 addFeed( feedId, spwId ) ; 513 } 514 virtual void leaveIfNo(const uInt recordNo, uInt columnValue) { 515 } 516 virtual void enterSrcType(const uInt recordNo, Int columnValue) { 517 //printf("%u: SrcType: %d\n", recordNo, columnValue); 518 519 Int stateId = addState( columnValue ) ; 520 521 // put value 522 *stateIdRF = stateId ; 523 } 524 virtual void leaveSrcType(const uInt recordNo, Int columnValue) { 525 } 526 virtual void enterCycleNo(const uInt recordNo, uInt columnValue) { 527 //printf("%u: CycleNo: %u\n", recordNo, columnValue); 528 } 529 virtual void leaveCycleNo(const uInt recordNo, uInt columnValue) { 530 } 531 virtual void enterTime(const uInt recordNo, Double columnValue) { 532 //printf("%u: Time: %f\n", recordNo, columnValue); 533 534 Double timeSec = columnValue * 86400.0 ; 535 Double interval = intervalCol.asdouble( recordNo ) ; 536 537 if ( ptName.empty() ) { 538 Vector<Double> dir = directionCol( recordNo ) ; 539 Vector<Double> rate = scanRateCol( recordNo ) ; 540 if ( anyNE( rate, 0.0 ) ) { 541 Matrix<Double> msdir( 2, 2 ) ; 542 msdir.column( 0 ) = dir ; 543 msdir.column( 1 ) = rate ; 544 addPointing( timeSec, interval, msdir ) ; 545 } 546 else { 547 Matrix<Double> msdir( 2, 1 ) ; 548 msdir.column( 0 ) = dir ; 549 addPointing( timeSec, interval, msdir ) ; 550 } 551 } 552 553 // put value 554 *timeRF = timeSec ; 555 *timeCentroidRF = timeSec ; 556 *intervalRF = interval ; 557 *exposureRF = interval ; 558 } 559 virtual void leaveTime(const uInt recordNo, Double columnValue) { 560 if ( holder.nPol() > 0 ) { 561 Vector<Float> w = holder.getWeight() ; 562 Cube<Bool> c = holder.getFlagCategory() ; 563 Bool flr = holder.getFlagRow() ; 564 Matrix<Bool> fl = holder.getFlag() ; 565 Vector<uInt> polnos = holder.polNos() ; 566 Int polId = addPolarization( polnos ) ; 567 Int ddId = addDataDescription( polId, spwId ) ; 568 569 // put field 570 *dataDescIdRF = ddId ; 571 *flagRowRF = flr ; 572 weightRF.define( w ) ; 573 sigmaRF.define( w ) ; 574 flagCategoryRF.define( c ) ; 575 flagRF.define( fl ) ; 576 if ( useFloat ) { 577 Matrix<Float> sp = holder.getData() ; 578 floatDataRF.define( sp ) ; 579 } 580 else { 581 Matrix<Complex> sp = holder.getComplexData() ; 582 dataRF.define( sp ) ; 583 } 584 585 // commit row 586 row.put( rowidx ) ; 587 rowidx++ ; 588 589 // reset holder 590 holder.reset() ; 591 } 592 if ( tcalKey != -1 ) { 593 tcalNotYet[tcalKey] = False ; 594 tcalKey = -1 ; 595 } 596 } 597 virtual void enterPolNo(const uInt recordNo, uInt columnValue) { 598 //printf("%u: PolNo: %d\n", recordNo, columnValue); 599 uInt tcalId = tcalIdCol.asuInt( recordNo ) ; 600 if ( tcalKey == -1 ) { 601 tcalKey = tcalId ; 602 } 603 if ( tcalNotYet[tcalKey] ) { 604 map< Int,Vector<uInt> >::iterator itr = tcalIdRec.find( tcalKey ) ; 605 if ( itr != tcalIdRec.end() ) { 606 Vector<uInt> ids = itr->second ; 607 uInt nrow = ids.nelements() ; 608 ids.resize( nrow+1, True ) ; 609 ids[nrow] = tcalId ; 610 tcalIdRec.erase( tcalKey ) ; 611 tcalIdRec[tcalKey] = ids ; 612 } 613 else { 614 Vector<uInt> rows( 1, tcalId ) ; 615 tcalIdRec[tcalKey] = rows ; 616 } 617 } 618 map< Int,Vector<uInt> >::iterator itr = tcalRowRec.find( tcalKey ) ; 619 if ( itr != tcalRowRec.end() ) { 620 Vector<uInt> rows = itr->second ; 621 uInt nrow = rows.nelements() ; 622 rows.resize( nrow+1, True ) ; 623 rows[nrow] = recordNo ; 624 tcalRowRec.erase( tcalKey ) ; 625 tcalRowRec[tcalKey] = rows ; 626 } 627 else { 628 Vector<uInt> rows( 1, recordNo ) ; 629 tcalRowRec[tcalKey] = rows ; 630 } 631 } 632 virtual void leavePolNo(const uInt recordNo, uInt columnValue) { 633 } 634 635 virtual Bool visitRecord(const uInt recordNo, 636 const String &fieldName, 637 const uInt beamNo, 638 const uInt scanNo, 639 const uInt ifNo, 640 const Int srcType, 641 const uInt cycleNo, 642 const Double time, 643 const uInt polNo) { 644 //printf("%u: %s, %u, %u, %u, %d, %u, %f, %d\n", recordNo, 645 // fieldName.c_str(), beamNo, scanNo, ifNo, srcType, cycleNo, time, polNo); 646 647 Vector<Float> sp = spectraCol( recordNo ) ; 648 Vector<uChar> tmp = flagtraCol( recordNo ) ; 649 Vector<Bool> fl( tmp.shape() ) ; 650 convertArray( fl, tmp ) ; 651 Bool flr = (Bool)flagRowCol.asuInt( recordNo ) ; 652 holder.accumulate( polNo, sp, fl, flr ) ; 653 654 return True ; 655 } 656 657 virtual void finish() { 658 BaseMSWriterVisitor::finish(); 659 //printf("Total: %u\n", count); 660 661 // remove rows 662 if ( ms.nrow() > rowidx ) { 663 uInt numRemove = ms.nrow() - rowidx ; 664 //cout << "numRemove = " << numRemove << endl ; 665 Vector<uInt> rows( numRemove ) ; 666 indgen( rows, rowidx ) ; 667 ms.removeRow( rows ) ; 668 } 669 670 // fill empty SPECTRAL_WINDOW rows 671 infillSpectralWindow() ; 672 } 673 674 void dataColumnName( String name ) 675 { 676 TableRecord &r = row.record() ; 677 if ( name == "DATA" ) { 678 useFloat = False ; 679 dataRF.attachToRecord( r, name ) ; 680 } 681 else if ( name == "FLOAT_DATA" ) { 682 useFloat = True ; 683 floatDataRF.attachToRecord( r, name ) ; 684 } 685 } 686 void pointingTableName( String name ) { 687 ptName = name ; 688 } 689 void setSourceRecord( Record &r ) { 690 srcRec = r ; 691 } 692 map< Int,Vector<uInt> > &getTcalIdRecord() { return tcalIdRec ; } 693 map< Int,Vector<uInt> > &getTcalRowRecord() { return tcalRowRec ; } 694 private: 695 void addField( Int &fid, String &fname, String &srcName, 696 Vector<Double> &sdir, Vector<Double> &srate, 697 Double &tSec ) 698 { 699 uInt nrow = fieldtab.nrow() ; 700 while( (Int)nrow <= fid ) { 701 fieldtab.addRow( 1, True ) ; 702 nrow++ ; 703 } 704 705 Matrix<Double> dir ; 706 Int numPoly = 0 ; 707 if ( anyNE( srate, 0.0 ) ) { 708 dir.resize( 2, 2 ) ; 709 dir.column( 0 ) = sdir ; 710 dir.column( 1 ) = srate ; 711 numPoly = 1 ; 712 } 713 else { 714 dir.resize( 2, 1 ) ; 715 dir.column( 0 ) = sdir ; 716 } 717 srcId = srcRec.asInt( srcName ) ; 718 719 TableRow tr( fieldtab ) ; 720 TableRecord &r = tr.record() ; 721 putField( "NAME", r, fname ) ; 722 putField( "NUM_POLY", r, numPoly ) ; 723 putField( "TIME", r, tSec ) ; 724 putField( "SOURCE_ID", r, srcId ) ; 725 defineField( "DELAY_DIR", r, dir ) ; 726 defineField( "REFERENCE_DIR", r, dir ) ; 727 defineField( "PHASE_DIR", r, dir ) ; 728 tr.put( fid ) ; 729 730 // for POINTING table 731 *poNameRF = fname ; 732 } 733 Int addState( Int &id ) 734 { 735 String obsMode ; 736 Bool isSignal ; 737 Double tnoise ; 738 Double tload ; 739 queryType( id, obsMode, isSignal, tnoise, tload ) ; 740 741 String key = obsMode+"_"+String::toString( subscan ) ; 742 Int idx = -1 ; 743 uInt nEntry = stateEntry.nelements() ; 744 for ( uInt i = 0 ; i < nEntry ; i++ ) { 745 if ( stateEntry[i] == key ) { 746 idx = i ; 747 break ; 748 } 749 } 750 if ( idx == -1 ) { 751 uInt nrow = statetab.nrow() ; 752 statetab.addRow( 1, True ) ; 753 TableRow tr( statetab ) ; 754 TableRecord &r = tr.record() ; 755 putField( "OBS_MODE", r, obsMode ) ; 756 putField( "SIG", r, isSignal ) ; 757 isSignal = !isSignal ; 758 putField( "REF", r, isSignal ) ; 759 putField( "CAL", r, tnoise ) ; 760 putField( "LOAD", r, tload ) ; 761 tr.put( nrow ) ; 762 idx = nrow ; 763 764 stateEntry.resize( nEntry+1, True ) ; 765 stateEntry[nEntry] = key ; 766 } 767 subscan++ ; 768 769 return idx ; 770 } 771 void addPointing( Double &tSec, Double &interval, Matrix<Double> &dir ) 772 { 773 uInt nrow = potab.nrow() ; 774 potab.addRow( 1, True ) ; 775 776 *poNumPolyRF = dir.ncolumn() - 1 ; 777 *poTimeRF = tSec ; 778 *poTimeOriginRF = tSec ; 779 *poIntervalRF = interval ; 780 poDirectionRF.define( dir ) ; 781 poTargetRF.define( dir ) ; 782 porow.put( nrow ) ; 783 } 784 Int addPolarization( Vector<uInt> &nos ) 785 { 786 Int idx = -1 ; 787 uInt nEntry = polEntry.size() ; 788 for ( uInt i = 0 ; i < nEntry ; i++ ) { 789 if ( polEntry[i].conform( nos ) && allEQ( polEntry[i], nos ) ) { 790 idx = i ; 791 break ; 792 } 793 } 794 795 Int numCorr ; 796 Vector<Int> corrType ; 797 Matrix<Int> corrProduct ; 798 polProperty( nos, numCorr, corrType, corrProduct ) ; 799 800 if ( idx == -1 ) { 801 uInt nrow = poltab.nrow() ; 802 poltab.addRow( 1, True ) ; 803 TableRow tr( poltab ) ; 804 TableRecord &r = tr.record() ; 805 putField( "NUM_CORR", r, numCorr ) ; 806 defineField( "CORR_TYPE", r, corrType ) ; 807 defineField( "CORR_PRODUCT", r, corrProduct ) ; 808 tr.put( nrow ) ; 809 idx = nrow ; 810 811 polEntry.resize( nEntry+1 ) ; 812 polEntry[nEntry] = nos ; 813 } 814 815 return idx ; 816 } 817 Int addDataDescription( Int pid, Int sid ) 818 { 819 Int idx = -1 ; 820 uInt nEntry = ddEntry.nrow() ; 821 Vector<Int> key( 2 ) ; 822 key[0] = pid ; 823 key[1] = sid ; 824 for ( uInt i = 0 ; i < nEntry ; i++ ) { 825 if ( allEQ( ddEntry.row(i), key ) ) { 826 idx = i ; 827 break ; 828 } 829 } 830 831 if ( idx == -1 ) { 832 uInt nrow = ddtab.nrow() ; 833 ddtab.addRow( 1, True ) ; 834 TableRow tr( ddtab ) ; 835 TableRecord &r = tr.record() ; 836 putField( "POLARIZATION_ID", r, pid ) ; 837 putField( "SPECTRAL_WINDOW_ID", r, sid ) ; 838 tr.put( nrow ) ; 839 idx = nrow ; 840 841 ddEntry.resize( nEntry+1, 2, True ) ; 842 ddEntry.row(nEntry) = key ; 843 } 844 845 return idx ; 846 } 847 void infillSpectralWindow() 848 { 849 ROScalarColumn<Int> nchanCol( spwtab, "NUM_CHAN" ) ; 850 Vector<Int> nchan = nchanCol.getColumn() ; 851 TableRow tr( spwtab ) ; 852 TableRecord &r = tr.record() ; 853 Int mfr = 1 ; 854 Vector<Double> dummy( 1, 0.0 ) ; 855 putField( "MEAS_FREQ_REF", r, mfr ) ; 856 defineField( "CHAN_FREQ", r, dummy ) ; 857 defineField( "CHAN_WIDTH", r, dummy ) ; 858 defineField( "EFFECTIVE_BW", r, dummy ) ; 859 defineField( "RESOLUTION", r, dummy ) ; 860 861 for ( uInt i = 0 ; i < spwtab.nrow() ; i++ ) { 862 if ( nchan[i] == 0 ) 863 tr.put( i ) ; 864 } 865 } 866 void addSpectralWindow( Int sid, uInt fid ) 867 { 868 if ( !processedFreqId[fid] ) { 869 uInt nrow = spwtab.nrow() ; 870 while( (Int)nrow <= sid ) { 871 spwtab.addRow( 1, True ) ; 872 nrow++ ; 873 } 874 processedFreqId[fid] = True ; 875 } 876 877 Double rp = refpix[fid] ; 878 Double rv = refval[fid] ; 879 Double ic = increment[fid] ; 880 881 Int mfrInt = (Int)freqframe ; 882 Int nchan = holder.nChan() ; 883 Double bw = nchan * abs( ic ) ; 884 Double reffreq = rv - rp * ic ; 885 Int netsb = 0 ; // USB->0, LSB->1 886 if ( ic < 0 ) 887 netsb = 1 ; 888 Vector<Double> res( nchan, abs(ic) ) ; 889 Vector<Double> chanf( nchan ) ; 890 indgen( chanf, reffreq, ic ) ; 891 892 TableRow tr( spwtab ) ; 893 TableRecord &r = tr.record() ; 894 putField( "MEAS_FREQ_REF", r, mfrInt ) ; 895 putField( "NUM_CHAN", r, nchan ) ; 896 putField( "TOTAL_BANDWIDTH", r, bw ) ; 897 putField( "REF_FREQUENCY", r, reffreq ) ; 898 putField( "NET_SIDEBAND", r, netsb ) ; 899 defineField( "RESOLUTION", r, res ) ; 900 defineField( "CHAN_WIDTH", r, res ) ; 901 defineField( "EFFECTIVE_BW", r, res ) ; 902 defineField( "CHAN_FREQ", r, chanf ) ; 903 tr.put( sid ) ; 904 } 905 void addFeed( Int fid, Int sid ) 906 { 907 Int idx = -1 ; 908 uInt nEntry = feedEntry.nrow() ; 909 Vector<Int> key( 2 ) ; 910 key[0] = fid ; 911 key[1] = sid ; 912 for ( uInt i = 0 ; i < nEntry ; i++ ) { 913 if ( allEQ( feedEntry.row(i), key ) ) { 914 idx = i ; 915 break ; 916 } 917 } 918 919 if ( idx == -1 ) { 920 uInt nrow = feedtab.nrow() ; 921 feedtab.addRow( 1, True ) ; 922 Int numReceptors = 2 ; 923 Vector<String> polType( numReceptors ) ; 924 Matrix<Double> beamOffset( 2, numReceptors, 0.0 ) ; 925 Vector<Double> receptorAngle( numReceptors, 0.0 ) ; 926 if ( poltype == "linear" ) { 927 polType[0] = "X" ; 928 polType[1] = "Y" ; 929 } 930 else if ( poltype == "circular" ) { 931 polType[0] = "R" ; 932 polType[1] = "L" ; 933 } 934 else { 935 polType[0] = "X" ; 936 polType[1] = "Y" ; 937 } 938 Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ; 939 940 TableRow tr( feedtab ) ; 941 TableRecord &r = tr.record() ; 942 putField( "FEED_ID", r, fid ) ; 943 putField( "BEAM_ID", r, fid ) ; 944 Int tmp = 0 ; 945 putField( "ANTENNA_ID", r, tmp ) ; 946 putField( "SPECTRAL_WINDOW_ID", r, sid ) ; 947 putField( "NUM_RECEPTORS", r, numReceptors ) ; 948 defineField( "POLARIZATION_TYPE", r, polType ) ; 949 defineField( "BEAM_OFFSET", r, beamOffset ) ; 950 defineField( "RECEPTOR_ANGLE", r, receptorAngle ) ; 951 defineField( "POL_RESPONSE", r, polResponse ) ; 952 tr.put( nrow ) ; 953 954 feedEntry.resize( nEntry+1, 2, True ) ; 955 feedEntry.row( nEntry ) = key ; 956 } 957 } 958 void makePolMap() 959 { 960 const TableRecord &keys = table.keywordSet() ; 961 poltype = keys.asString( "POLTYPE" ) ; 962 963 if ( poltype == "stokes" ) { 964 polmap.resize( 4 ) ; 965 polmap[0] = Stokes::I ; 966 polmap[1] = Stokes::Q ; 967 polmap[2] = Stokes::U ; 968 polmap[3] = Stokes::V ; 969 } 970 else if ( poltype == "linear" ) { 971 polmap.resize( 4 ) ; 972 polmap[0] = Stokes::XX ; 973 polmap[1] = Stokes::YY ; 974 polmap[2] = Stokes::XY ; 975 polmap[3] = Stokes::YX ; 976 } 977 else if ( poltype == "circular" ) { 978 polmap.resize( 4 ) ; 979 polmap[0] = Stokes::RR ; 980 polmap[1] = Stokes::LL ; 981 polmap[2] = Stokes::RL ; 982 polmap[3] = Stokes::LR ; 983 } 984 else if ( poltype == "linpol" ) { 985 polmap.resize( 2 ) ; 986 polmap[0] = Stokes::Plinear ; 987 polmap[1] = Stokes::Pangle ; 988 } 989 else { 990 polmap.resize( 0 ) ; 991 } 992 } 993 void initFrequencies() 994 { 995 const TableRecord &keys = table.keywordSet() ; 996 Table tab = keys.asTable( "FREQUENCIES" ) ; 997 ROScalarColumn<uInt> idcol( tab, "ID" ) ; 998 ROScalarColumn<Double> rpcol( tab, "REFPIX" ) ; 999 ROScalarColumn<Double> rvcol( tab, "REFVAL" ) ; 1000 ROScalarColumn<Double> iccol( tab, "INCREMENT" ) ; 1001 Vector<uInt> id = idcol.getColumn() ; 1002 Vector<Double> rp = rpcol.getColumn() ; 1003 Vector<Double> rv = rvcol.getColumn() ; 1004 Vector<Double> ic = iccol.getColumn() ; 1005 for ( uInt i = 0 ; i < id.nelements() ; i++ ) { 1006 processedFreqId.insert( pair<uInt,Bool>( id[i], False ) ) ; 1007 refpix.insert( pair<uInt,Double>( id[i], rp[i] ) ) ; 1008 refval.insert( pair<uInt,Double>( id[i], rv[i] ) ) ; 1009 increment.insert( pair<uInt,Double>( id[i], ic[i] ) ) ; 1010 } 1011 String frameStr = tab.keywordSet().asString( "BASEFRAME" ) ; 1012 MFrequency::getType( freqframe, frameStr ) ; 1013 } 1014 void attachSubtables() 1015 { 1016 const TableRecord &keys = table.keywordSet() ; 1017 TableRecord &mskeys = ms.rwKeywordSet() ; 1018 1019 // FIELD table 1020 fieldtab = mskeys.asTable( "FIELD" ) ; 1021 1022 // SPECTRAL_WINDOW table 1023 spwtab = mskeys.asTable( "SPECTRAL_WINDOW" ) ; 1024 1025 // POINTING table 1026 potab = mskeys.asTable( "POINTING" ) ; 1027 1028 // POLARIZATION table 1029 poltab = mskeys.asTable( "POLARIZATION" ) ; 1030 1031 // DATA_DESCRIPTION table 1032 ddtab = mskeys.asTable( "DATA_DESCRIPTION" ) ; 1033 1034 // STATE table 1035 statetab = mskeys.asTable( "STATE" ) ; 1036 1037 // FEED table 1038 feedtab = mskeys.asTable( "FEED" ) ; 1039 } 1040 void attachMain() 1041 { 1042 TableRecord &r = row.record() ; 1043 dataDescIdRF.attachToRecord( r, "DATA_DESC_ID" ) ; 1044 flagRowRF.attachToRecord( r, "FLAG_ROW" ) ; 1045 weightRF.attachToRecord( r, "WEIGHT" ) ; 1046 sigmaRF.attachToRecord( r, "SIGMA" ) ; 1047 flagCategoryRF.attachToRecord( r, "FLAG_CATEGORY" ) ; 1048 flagRF.attachToRecord( r, "FLAG" ) ; 1049 timeRF.attachToRecord( r, "TIME" ) ; 1050 timeCentroidRF.attachToRecord( r, "TIME_CENTROID" ) ; 1051 intervalRF.attachToRecord( r, "INTERVAL" ) ; 1052 exposureRF.attachToRecord( r, "EXPOSURE" ) ; 1053 fieldIdRF.attachToRecord( r, "FIELD_ID" ) ; 1054 feed1RF.attachToRecord( r, "FEED1" ) ; 1055 feed2RF.attachToRecord( r, "FEED2" ) ; 1056 scanNumberRF.attachToRecord( r, "SCAN_NUMBER" ) ; 1057 stateIdRF.attachToRecord( r, "STATE_ID" ) ; 1058 1059 // constant values 1060 Int id = 0 ; 1061 RecordFieldPtr<Int> intRF( r, "OBSERVATION_ID" ) ; 1062 *intRF = 0 ; 1063 intRF.attachToRecord( r, "ANTENNA1" ) ; 1064 *intRF = 0 ; 1065 intRF.attachToRecord( r, "ANTENNA2" ) ; 1066 *intRF = 0 ; 1067 intRF.attachToRecord( r, "ARRAY_ID" ) ; 1068 *intRF = 0 ; 1069 intRF.attachToRecord( r, "PROCESSOR_ID" ) ; 1070 *intRF = 0 ; 1071 RecordFieldPtr< Vector<Double> > arrayRF( r, "UVW" ) ; 1072 arrayRF.define( Vector<Double>( 3, 0.0 ) ) ; 1073 } 1074 void attachPointing() 1075 { 1076 porow = TableRow( potab ) ; 1077 TableRecord &r = porow.record() ; 1078 poNumPolyRF.attachToRecord( r, "NUM_POLY" ) ; 1079 poTimeRF.attachToRecord( r, "TIME" ) ; 1080 poTimeOriginRF.attachToRecord( r, "TIME_ORIGIN" ) ; 1081 poIntervalRF.attachToRecord( r, "INTERVAL" ) ; 1082 poNameRF.attachToRecord( r, "NAME" ) ; 1083 poDirectionRF.attachToRecord( r, "DIRECTION" ) ; 1084 poTargetRF.attachToRecord( r, "TARGET" ) ; 1085 1086 // constant values 1087 RecordFieldPtr<Int> antIdRF( r, "ANTENNA_ID" ) ; 1088 *antIdRF = 0 ; 1089 RecordFieldPtr<Bool> trackingRF( r, "TRACKING" ) ; 1090 *trackingRF = True ; 1091 } 1092 void queryType( Int type, String &stype, Bool &b, Double &t, Double &l ) 1093 { 1094 t = 0.0 ; 1095 l = 0.0 ; 1096 1097 String sep1="#" ; 1098 String sep2="," ; 1099 String target="OBSERVE_TARGET" ; 1100 String atmcal="CALIBRATE_TEMPERATURE" ; 1101 String onstr="ON_SOURCE" ; 1102 String offstr="OFF_SOURCE" ; 1103 String pswitch="POSITION_SWITCH" ; 1104 String nod="NOD" ; 1105 String fswitch="FREQUENCY_SWITCH" ; 1106 String sigstr="SIG" ; 1107 String refstr="REF" ; 1108 String unspecified="UNSPECIFIED" ; 1109 String ftlow="LOWER" ; 1110 String fthigh="HIGHER" ; 1111 switch ( type ) { 1112 case SrcType::PSON: 1113 stype = target+sep1+onstr+sep2+pswitch ; 1114 b = True ; 1115 break ; 1116 case SrcType::PSOFF: 1117 stype = target+sep1+offstr+sep2+pswitch ; 1118 b = False ; 1119 break ; 1120 case SrcType::NOD: 1121 stype = target+sep1+onstr+sep2+nod ; 1122 b = True ; 1123 break ; 1124 case SrcType::FSON: 1125 stype = target+sep1+onstr+sep2+fswitch+sep1+sigstr ; 1126 b = True ; 1127 break ; 1128 case SrcType::FSOFF: 1129 stype = target+sep1+onstr+sep2+fswitch+sep1+refstr ; 1130 b = False ; 1131 break ; 1132 case SrcType::SKY: 1133 stype = atmcal+sep1+offstr+sep2+unspecified ; 1134 b = False ; 1135 break ; 1136 case SrcType::HOT: 1137 stype = atmcal+sep1+offstr+sep2+unspecified ; 1138 b = False ; 1139 break ; 1140 case SrcType::WARM: 1141 stype = atmcal+sep1+offstr+sep2+unspecified ; 1142 b = False ; 1143 break ; 1144 case SrcType::COLD: 1145 stype = atmcal+sep1+offstr+sep2+unspecified ; 1146 b = False ; 1147 break ; 1148 case SrcType::PONCAL: 1149 stype = atmcal+sep1+onstr+sep2+pswitch ; 1150 b = True ; 1151 break ; 1152 case SrcType::POFFCAL: 1153 stype = atmcal+sep1+offstr+sep2+pswitch ; 1154 b = False ; 1155 break ; 1156 case SrcType::NODCAL: 1157 stype = atmcal+sep1+onstr+sep2+nod ; 1158 b = True ; 1159 break ; 1160 case SrcType::FONCAL: 1161 stype = atmcal+sep1+onstr+sep2+fswitch+sep1+sigstr ; 1162 b = True ; 1163 break ; 1164 case SrcType::FOFFCAL: 1165 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+refstr ; 1166 b = False ; 1167 break ; 1168 case SrcType::FSLO: 1169 stype = target+sep1+onstr+sep2+fswitch+sep1+ftlow ; 1170 b = True ; 1171 break ; 1172 case SrcType::FLOOFF: 1173 stype = target+sep1+offstr+sep2+fswitch+sep1+ftlow ; 1174 b = False ; 1175 break ; 1176 case SrcType::FLOSKY: 1177 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ; 1178 b = False ; 1179 break ; 1180 case SrcType::FLOHOT: 1181 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ; 1182 b = False ; 1183 break ; 1184 case SrcType::FLOWARM: 1185 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ; 1186 b = False ; 1187 break ; 1188 case SrcType::FLOCOLD: 1189 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ; 1190 b = False ; 1191 break ; 1192 case SrcType::FSHI: 1193 stype = target+sep1+onstr+sep2+fswitch+sep1+fthigh ; 1194 b = True ; 1195 break ; 1196 case SrcType::FHIOFF: 1197 stype = target+sep1+offstr+sep2+fswitch+sep1+fthigh ; 1198 b = False ; 1199 break ; 1200 case SrcType::FHISKY: 1201 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ; 1202 b = False ; 1203 break ; 1204 case SrcType::FHIHOT: 1205 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ; 1206 b = False ; 1207 break ; 1208 case SrcType::FHIWARM: 1209 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ; 1210 b = False ; 1211 break ; 1212 case SrcType::FHICOLD: 1213 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ; 1214 b = False ; 1215 break ; 1216 case SrcType::SIG: 1217 stype = target+sep1+onstr+sep2+unspecified ; 1218 b = True ; 1219 break ; 1220 case SrcType::REF: 1221 stype = target+sep1+offstr+sep2+unspecified ; 1222 b = False ; 1223 break ; 1224 default: 1225 stype = unspecified ; 1226 b = True ; 1227 break ; 1228 } 1229 } 1230 void polProperty( Vector<uInt> &nos, Int &n, Vector<Int> &c, Matrix<Int> &cp ) 1231 { 1232 n = nos.nelements() ; 1233 c.resize( n ) ; 1234 for ( Int i = 0 ; i < n ; i++ ) 1235 c[i] = (Int)polmap[nos[i]] ; 1236 cp.resize( 2, n ) ; 1237 if ( n == 1 ) 1238 cp = 0 ; 1239 else if ( n == 2 ) { 1240 cp.column( 0 ) = 0 ; 1241 cp.column( 1 ) = 1 ; 1242 } 1243 else { 1244 cp.column( 0 ) = 0 ; 1245 cp.column( 1 ) = 1 ; 1246 cp( 0, 1 ) = 0 ; 1247 cp( 1, 1 ) = 1 ; 1248 cp( 0, 2 ) = 1 ; 1249 cp( 1, 2 ) = 0 ; 1250 } 1251 } 1252 void initTcal() 1253 { 1254 const TableRecord &rec = table.keywordSet() ; 1255 Table tcalTable = rec.asTable( "TCAL" ) ; 1256 ROScalarColumn<uInt> idCol( tcalTable, "ID" ) ; 1257 Vector<uInt> id = idCol.getColumn() ; 1258 uInt maxId = max( id ) ; 1259 tcalNotYet.resize( maxId+1 ) ; 1260 tcalNotYet = True ; 1261 tcalKey = -1 ; 1262 } 1263 1264 Table &ms; 1265 TableRow row; 1266 uInt rowidx; 1267 String fieldName; 1268 Int fieldId; 1269 Int srcId; 1270 Int defaultFieldId; 1271 Int spwId; 1272 Int feedId; 1273 Int subscan; 1274 PolarizedComponentHolder holder; 1275 String ptName; 1276 Bool useFloat; 1277 String poltype; 1278 Vector<Stokes::StokesTypes> polmap; 1279 1280 // MS subtables 1281 Table spwtab; 1282 Table statetab; 1283 Table ddtab; 1284 Table poltab; 1285 Table fieldtab; 1286 Table feedtab; 1287 Table potab; 1288 1289 // Scantable MAIN columns 1290 ROArrayColumn<Float> spectraCol; 1291 ROArrayColumn<Double> directionCol,scanRateCol,sourceDirectionCol; 1292 ROArrayColumn<uChar> flagtraCol; 1293 ROTableColumn tcalIdCol,intervalCol,flagRowCol,timeCol,freqIdCol, 1294 sourceNameCol,fieldNameCol; 1295 1296 // MS MAIN columns 1297 RecordFieldPtr<Int> dataDescIdRF,fieldIdRF,feed1RF,feed2RF, 1298 scanNumberRF,stateIdRF; 1299 RecordFieldPtr<Bool> flagRowRF; 1300 RecordFieldPtr<Double> timeRF,timeCentroidRF,intervalRF,exposureRF; 1301 RecordFieldPtr< Vector<Float> > weightRF,sigmaRF; 1302 RecordFieldPtr< Cube<Bool> > flagCategoryRF; 1303 RecordFieldPtr< Matrix<Bool> > flagRF; 1304 RecordFieldPtr< Matrix<Float> > floatDataRF; 1305 RecordFieldPtr< Matrix<Complex> > dataRF; 1306 1307 // MS POINTING columns 1308 TableRow porow; 1309 RecordFieldPtr<Int> poNumPolyRF ; 1310 RecordFieldPtr<Double> poTimeRF, 1311 poTimeOriginRF, 1312 poIntervalRF ; 1313 RecordFieldPtr<String> poNameRF ; 1314 RecordFieldPtr< Matrix<Double> > poDirectionRF, 1315 poTargetRF ; 1316 1317 Vector<String> stateEntry; 1318 Matrix<Int> ddEntry; 1319 Matrix<Int> feedEntry; 1320 vector< Vector<uInt> > polEntry; 1321 map<uInt,Bool> processedFreqId; 1322 map<uInt,Double> refpix; 1323 map<uInt,Double> refval; 1324 map<uInt,Double> increment; 1325 MFrequency::Types freqframe; 1326 Record srcRec; 1327 map< Int,Vector<uInt> > tcalIdRec; 1328 map< Int,Vector<uInt> > tcalRowRec; 1329 Int tcalKey; 1330 Vector<Bool> tcalNotYet; 1331 }; 56 1332 57 1333 MSWriter::MSWriter(CountedPtr<Scantable> stable) … … 79 1355 delete mstable_ ; 80 1356 } 81 1357 82 1358 bool MSWriter::write(const string& filename, const Record& rec) 83 1359 { 84 1360 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ; 85 //double startSec = mathutil::gettimeofday_sec() ;86 //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;1361 //double startSec = mathutil::gettimeofday_sec() ; 1362 //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ; 87 1363 88 1364 filename_ = filename ; … … 134 1410 fillWeather() ; 135 1411 136 // MAIN 137 // Iterate over several ids 138 Vector<uInt> processedFreqId( 0 ) ; 139 Int defaultFieldId = 0 ; 140 141 // row based 142 TableRow row( *mstable_ ) ; 143 TableRecord &trec = row.record() ; 144 NoticeTarget *dataRF = 0 ; 145 if ( useFloatData_ ) 146 dataRF = new RecordFieldPtr< Array<Float> >( trec, "FLOAT_DATA" ) ; 147 else if ( useData_ ) 148 dataRF = new RecordFieldPtr< Array<Complex> >( trec, "DATA" ) ; 149 RecordFieldPtr< Array<Bool> > flagRF( trec, "FLAG" ) ; 150 RecordFieldPtr<Bool> flagrowRF( trec, "FLAG_ROW" ) ; 151 RecordFieldPtr<Double> timeRF( trec, "TIME" ) ; 152 RecordFieldPtr<Double> timecRF( trec, "TIME_CENTROID" ) ; 153 RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ; 154 RecordFieldPtr<Double> exposureRF( trec, "EXPOSURE" ) ; 155 RecordFieldPtr< Array<Float> > weightRF( trec, "WEIGHT" ) ; 156 RecordFieldPtr< Array<Float> > sigmaRF( trec, "SIGMA" ) ; 157 RecordFieldPtr<Int> ddidRF( trec, "DATA_DESC_ID" ) ; 158 RecordFieldPtr<Int> stateidRF( trec, "STATE_ID" ) ; 159 RecordFieldPtr< Array<Bool> > flagcatRF( trec, "FLAG_CATEGORY" ) ; 160 161 // OBSERVATION_ID is always 0 162 RecordFieldPtr<Int> intRF( trec, "OBSERVATION_ID" ) ; 163 *intRF = 0 ; 164 165 // ANTENNA1 and ANTENNA2 are always 0 166 intRF.attachToRecord( trec, "ANTENNA1" ) ; 167 *intRF = 0 ; 168 intRF.attachToRecord( trec, "ANTENNA2" ) ; 169 *intRF = 0 ; 170 171 // ARRAY_ID is tentatively set to 0 172 intRF.attachToRecord( trec, "ARRAY_ID" ) ; 173 *intRF = 0 ; 174 175 // PROCESSOR_ID is tentatively set to 0 176 intRF.attachToRecord( trec, "PROCESSOR_ID" ) ; 177 *intRF = 0 ; 178 179 // UVW is always [0,0,0] 180 RecordFieldPtr< Array<Double> > uvwRF( trec, "UVW" ) ; 181 *uvwRF = Vector<Double>( 3, 0.0 ) ; 182 183 // 184 // ITERATION: FIELDNAME 185 // 186 TableIterator iter0( table_->table(), "FIELDNAME" ) ; 187 while( !iter0.pastEnd() ) { 188 //Table t0( iter0.table() ) ; 189 Table t0 = iter0.table() ; 190 ROTableColumn sharedCol( t0, "FIELDNAME" ) ; 191 String fieldName = sharedCol.asString( 0 ) ; 192 sharedCol.attach( t0, "SRCNAME" ) ; 193 String srcName = sharedCol.asString( 0 ) ; 194 sharedCol.attach( t0, "TIME" ) ; 195 Double minTime = (Double)sharedCol.asdouble( 0 ) * 86400.0 ; // day->sec 196 ROArrayColumn<Double> scanRateCol( t0, "SCANRATE" ) ; 197 Vector<Double> scanRate = scanRateCol( 0 ) ; 198 String::size_type pos = fieldName.find( "__" ) ; 199 Int fieldId = -1 ; 200 if ( pos != String::npos ) { 201 // os_ << "fieldName.substr( pos+2 )=" << fieldName.substr( pos+2 ) << LogIO::POST ; 202 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ; 203 fieldName = fieldName.substr( 0, pos ) ; 204 } 205 else { 206 // os_ << "use default field id" << LogIO::POST ; 207 fieldId = defaultFieldId ; 208 defaultFieldId++ ; 209 } 210 // os_ << "fieldId" << fieldId << ": " << fieldName << LogIO::POST ; 211 212 // FIELD_ID 213 intRF.attachToRecord( trec, "FIELD_ID" ) ; 214 *intRF = fieldId ; 215 216 // 217 // ITERATION: BEAMNO 218 // 219 TableIterator iter1( t0, "BEAMNO" ) ; 220 while( !iter1.pastEnd() ) { 221 Table t1 = iter1.table() ; 222 sharedCol.attach( t1, "BEAMNO" ) ; 223 uInt beamNo = sharedCol.asuInt( 0 ) ; 224 // os_ << "beamNo = " << beamNo << LogIO::POST ; 225 226 // FEED1 and FEED2 227 intRF.attachToRecord( trec, "FEED1" ) ; 228 *intRF = beamNo ; 229 intRF.attachToRecord( trec, "FEED2" ) ; 230 *intRF = beamNo ; 231 232 // 233 // ITERATION: SCANNO 234 // 235 TableIterator iter2( t1, "SCANNO" ) ; 236 while( !iter2.pastEnd() ) { 237 Table t2 = iter2.table() ; 238 sharedCol.attach( t2, "SCANNO" ) ; 239 uInt scanNo = sharedCol.asuInt( 0 ) ; 240 // os_ << "scanNo = " << scanNo << LogIO::POST ; 241 242 // SCAN_NUMBER 243 // MS: 1-based 244 // Scantable: 0-based 245 intRF.attachToRecord( trec, "SCAN_NUMBER" ) ; 246 *intRF = scanNo + 1 ; 247 248 // 249 // ITERATION: IFNO 250 // 251 TableIterator iter3( t2, "IFNO" ) ; 252 while( !iter3.pastEnd() ) { 253 Table t3 = iter3.table() ; 254 sharedCol.attach( t3, "IFNO" ) ; 255 uInt ifNo = sharedCol.asuInt( 0 ) ; 256 // os_ << "ifNo = " << ifNo << LogIO::POST ; 257 sharedCol.attach( t3, "FREQ_ID" ) ; 258 uInt freqId = sharedCol.asuInt( 0 ) ; 259 // os_ << "freqId = " << freqId << LogIO::POST ; 260 Int subscan = 1 ; // 1-base 261 // 262 // ITERATION: SRCTYPE 263 // 264 TableIterator iter4( t3, "SRCTYPE" ) ; 265 while( !iter4.pastEnd() ) { 266 Table t4 = iter4.table() ; 267 sharedCol.attach( t4, "SRCTYPE" ) ; 268 Int srcType = sharedCol.asInt( 0 ) ; 269 Int stateId = addState( srcType, subscan ) ; 270 *stateidRF = stateId ; 271 // 272 // ITERATION: CYCLENO and TIME 273 // 274 Block<String> cols( 2 ) ; 275 cols[0] = "CYCLENO" ; 276 cols[1] = "TIME" ; 277 TableIterator iter5( t4, cols ) ; 278 while( !iter5.pastEnd() ) { 279 Table t5 = iter5.table().sort("POLNO") ; 280 //sharedCol.attach( t5, "CYCLENO" ) ; 281 //uInt cycleNo = sharedCol.asuInt( 0 ) ; 282 Int nrow = t5.nrow() ; 283 // os_ << "nrow = " << nrow << LogIO::POST ; 284 285 Vector<Int> polnos( nrow ) ; 286 indgen( polnos, 0 ) ; 287 Int polid = addPolarization( polnos ) ; 288 // os_ << "polid = " << polid << LogIO::POST ; 289 290 // DATA/FLOAT_DATA 291 ROArrayColumn<Float> specCol( t5, "SPECTRA" ) ; 292 ROArrayColumn<uChar> flagCol( t5, "FLAGTRA" ) ; 293 uInt nchan = specCol( 0 ).size() ; 294 IPosition cellshape( 2, nrow, nchan ) ; 295 if ( useFloatData_ ) { 296 // FLOAT_DATA 297 Matrix<Float> dataArr( cellshape ) ; 298 Matrix<Bool> flagArr( cellshape ) ; 299 Vector<Bool> tmpB ; 300 for ( Int ipol = 0 ; ipol < nrow ; ipol++ ) { 301 dataArr.row( ipol ) = specCol( ipol ) ; 302 tmpB.reference( flagArr.row( ipol ) ) ; 303 convertArray( tmpB, flagCol( ipol ) ) ; 304 } 305 ((RecordFieldPtr< Array<Float> > *)dataRF)->define( dataArr ) ; 306 307 // FLAG 308 flagRF.define( flagArr ) ; 309 } 310 else if ( useData_ ) { 311 // DATA 312 // assume nrow = 4 313 Matrix<Complex> dataArr( cellshape ) ; 314 Vector<Float> zeroIm( nchan, 0 ) ; 315 Matrix<Float> dummy( IPosition( 2, 2, nchan ) ) ; 316 dummy.row( 0 ) = specCol( 0 ) ; 317 dummy.row( 1 ) = zeroIm ; 318 dataArr.row( 0 ) = RealToComplex( dummy ) ; 319 dummy.row( 0 ) = specCol( 1 ) ; 320 dataArr.row( 3 ) = RealToComplex( dummy ) ; 321 dummy.row( 0 ) = specCol( 2 ) ; 322 dummy.row( 1 ) = specCol( 3 ) ; 323 dataArr.row( 1 ) = RealToComplex( dummy ) ; 324 dataArr.row( 2 ) = conj( dataArr.row( 1 ) ) ; 325 ((RecordFieldPtr< Array<Complex> > *)dataRF)->define( dataArr ) ; 326 327 328 // FLAG 329 Matrix<Bool> flagArr( cellshape ) ; 330 Vector<Bool> tmpB ; 331 tmpB.reference( flagArr.row( 0 ) ) ; 332 convertArray( tmpB, flagCol( 0 ) ) ; 333 tmpB.reference( flagArr.row( 3 ) ) ; 334 convertArray( tmpB, flagCol( 3 ) ) ; 335 tmpB.reference( flagArr.row( 1 ) ) ; 336 convertArray( tmpB, ( flagCol( 2 ) | flagCol( 3 ) ) ) ; 337 flagArr.row( 2 ) = flagArr.row( 1 ) ; 338 flagRF.define( flagArr ) ; 339 } 340 341 // FLAG_ROW 342 sharedCol.attach( t5, "FLAGROW" ) ; 343 Vector<uInt> flagRowArr( nrow ) ; 344 for ( Int irow = 0 ; irow < nrow ; irow++ ) 345 flagRowArr[irow] = sharedCol.asuInt( irow ) ; 346 *flagrowRF = anyNE( flagRowArr, (uInt)0 ) ; 347 348 // TIME and TIME_CENTROID 349 sharedCol.attach( t5, "TIME" ) ; 350 Double mTimeV = (Double)sharedCol.asdouble( 0 ) * 86400.0 ; // day -> sec 351 *timeRF = mTimeV ; 352 *timecRF = mTimeV ; 353 354 // INTERVAL and EXPOSURE 355 sharedCol.attach( t5, "INTERVAL" ) ; 356 Double interval = (Double)sharedCol.asdouble( 0 ) ; 357 *intervalRF = interval ; 358 *exposureRF = interval ; 359 360 // WEIGHT and SIGMA 361 // always 1 at the moment 362 Vector<Float> wArr( nrow, 1.0 ) ; 363 weightRF.define( wArr ) ; 364 sigmaRF.define( wArr ) ; 365 366 // add DATA_DESCRIPTION row 367 Int ddid = addDataDescription( polid, ifNo ) ; 368 // os_ << "ddid = " << ddid << LogIO::POST ; 369 *ddidRF = ddid ; 370 371 // for SYSCAL table 372 sharedCol.attach( t5, "TCAL_ID" ) ; 373 Vector<uInt> tcalIdArr( nrow ) ; 374 for ( Int irow = 0 ; irow < nrow ; irow++ ) 375 tcalIdArr[irow] = sharedCol.asuInt( irow ) ; 376 // os_ << "tcalIdArr = " << tcalIdArr << LogIO::POST ; 377 String key = String::toString( tcalIdArr[0] ) ; 378 if ( !tcalIdRec_.isDefined( key ) ) { 379 tcalIdRec_.define( key, tcalIdArr ) ; 380 tcalRowRec_.define( key, t5.rowNumbers() ) ; 381 } 382 else { 383 Vector<uInt> pastrows = tcalRowRec_.asArrayuInt( key ) ; 384 tcalRowRec_.define( key, concatenateArray( pastrows, t5.rowNumbers() ) ) ; 385 } 386 387 // for POINTING table 388 if ( ptTabName_ == "" ) { 389 ROArrayColumn<Double> dirCol( t5, "DIRECTION" ) ; 390 Vector<Double> dir = dirCol( 0 ) ; 391 dirCol.attach( t5, "SCANRATE" ) ; 392 Vector<Double> rate = dirCol( 0 ) ; 393 Matrix<Double> msDir( 2, 1 ) ; 394 msDir.column( 0 ) = dir ; 395 if ( anyNE( rate, 0.0 ) ) { 396 msDir.resize( 2, 2, True ) ; 397 msDir.column( 1 ) = rate ; 398 } 399 addPointing( fieldName, mTimeV, interval, msDir ) ; 400 } 401 402 // FLAG_CATEGORY is tentatively set 403 flagcatRF.define( Cube<Bool>( nrow, nchan, 1, False ) ) ; 404 405 // add row 406 mstable_->addRow( 1, True ) ; 407 row.put( mstable_->nrow()-1 ) ; 408 409 iter5.next() ; 410 } 411 412 iter4.next() ; 413 } 414 415 // add SPECTRAL_WINDOW row 416 if ( allNE( processedFreqId, freqId ) ) { 417 uInt vsize = processedFreqId.size() ; 418 processedFreqId.resize( vsize+1, True ) ; 419 processedFreqId[vsize] = freqId ; 420 addSpectralWindow( ifNo, freqId ) ; 421 } 422 423 iter3.next() ; 424 } 425 426 iter2.next() ; 427 } 428 429 // add FEED row 430 addFeed( beamNo ) ; 431 432 iter1.next() ; 433 } 434 435 // add FIELD row 436 addField( fieldId, fieldName, srcName, minTime, scanRate ) ; 437 438 iter0.next() ; 439 } 440 441 // delete tpoolr ; 442 delete dataRF ; 443 444 // SYSCAL 445 if ( isTcal_ ) 446 fillSysCal() ; 447 448 // fill empty SPECTRAL_WINDOW rows 449 infillSpectralWindow() ; 1412 /*** 1413 * Start iteration using TableVisitor 1414 ***/ 1415 { 1416 static const char *cols[] = { 1417 "FIELDNAME", "BEAMNO", "SCANNO", "IFNO", "SRCTYPE", "CYCLENO", "TIME", 1418 "POLNO", 1419 NULL 1420 }; 1421 static const TypeManagerImpl<uInt> tmUInt; 1422 static const TypeManagerImpl<Int> tmInt; 1423 static const TypeManagerImpl<Double> tmDouble; 1424 static const TypeManagerImpl<String> tmString; 1425 static const TypeManager *const tms[] = { 1426 &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmInt, NULL 1427 }; 1428 //double t0 = mathutil::gettimeofday_sec() ; 1429 MSWriterVisitor myVisitor(table_->table(),*mstable_); 1430 //double t1 = mathutil::gettimeofday_sec() ; 1431 //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ; 1432 String dataColName = "FLOAT_DATA" ; 1433 if ( useData_ ) 1434 dataColName = "DATA" ; 1435 myVisitor.dataColumnName( dataColName ) ; 1436 myVisitor.pointingTableName( ptTabName_ ) ; 1437 myVisitor.setSourceRecord( srcRec_ ) ; 1438 //double t2 = mathutil::gettimeofday_sec() ; 1439 traverseTable(table_->table(), cols, tms, &myVisitor); 1440 //double t3 = mathutil::gettimeofday_sec() ; 1441 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ; 1442 map< Int,Vector<uInt> > &idRec = myVisitor.getTcalIdRecord() ; 1443 map< Int,Vector<uInt> > &rowRec = myVisitor.getTcalRowRecord() ; 1444 1445 // SYSCAL 1446 if ( isTcal_ ) 1447 fillSysCal( idRec, rowRec ) ; 1448 } 1449 /*** 1450 * End iteration using TableVisitor 1451 ***/ 450 1452 451 1453 // ASDM tables … … 475 1477 } 476 1478 477 //double endSec = mathutil::gettimeofday_sec() ;478 //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1479 //double endSec = mathutil::gettimeofday_sec() ; 1480 //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 479 1481 480 1482 return True ; 481 1483 } 482 1484 483 1485 void MSWriter::init() 484 1486 { … … 697 1699 void MSWriter::fillObservation() 698 1700 { 699 //double startSec = mathutil::gettimeofday_sec() ;700 //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;1701 //double startSec = mathutil::gettimeofday_sec() ; 1702 //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ; 701 1703 702 1704 // only 1 row … … 726 1728 msObsCols.timeRangeMeas().put( 0, trange ) ; 727 1729 728 //double endSec = mathutil::gettimeofday_sec() ;729 //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1730 //double endSec = mathutil::gettimeofday_sec() ; 1731 //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 730 1732 } 1733 1734 void MSWriter::antennaProperty( String &name, String &m, String &t, Double &d ) 1735 { 1736 name.upcase() ; 1737 1738 m = "ALT-AZ" ; 1739 t = "GROUND-BASED" ; 1740 if ( name.matches( Regex( "DV[0-9]+$" ) ) 1741 || name.matches( Regex( "DA[0-9]+$" ) ) 1742 || name.matches( Regex( "PM[0-9]+$" ) ) ) 1743 d = 12.0 ; 1744 else if ( name.matches( Regex( "CM[0-9]+$" ) ) ) 1745 d = 7.0 ; 1746 else if ( name.contains( "GBT" ) ) 1747 d = 104.9 ; 1748 else if ( name.contains( "MOPRA" ) ) 1749 d = 22.0 ; 1750 else if ( name.contains( "PKS" ) || name.contains( "PARKS" ) ) 1751 d = 64.0 ; 1752 else if ( name.contains( "TIDBINBILLA" ) ) 1753 d = 70.0 ; 1754 else if ( name.contains( "CEDUNA" ) ) 1755 d = 30.0 ; 1756 else if ( name.contains( "HOBART" ) ) 1757 d = 26.0 ; 1758 else if ( name.contains( "APEX" ) ) 1759 d = 12.0 ; 1760 else if ( name.contains( "ASTE" ) ) 1761 d = 10.0 ; 1762 else if ( name.contains( "NRO" ) ) 1763 d = 45.0 ; 1764 else 1765 d = 1.0 ; 1766 } 731 1767 732 1768 void MSWriter::fillAntenna() 733 1769 { 734 //double startSec = mathutil::gettimeofday_sec() ;735 //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;1770 //double startSec = mathutil::gettimeofday_sec() ; 1771 //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ; 736 1772 737 1773 // only 1 row 738 mstable_->antenna().addRow( 1, True ) ; 739 MSAntennaColumns msAntCols( mstable_->antenna() ) ; 740 741 String hAntennaName = header_.antennaname ; 742 String::size_type pos = hAntennaName.find( "//" ) ; 1774 Table anttab = mstable_->antenna() ; 1775 anttab.addRow( 1, True ) ; 1776 1777 Table &table = table_->table() ; 1778 const TableRecord &keys = table.keywordSet() ; 1779 String hAntName = keys.asString( "AntennaName" ) ; 1780 String::size_type pos = hAntName.find( "//" ) ; 743 1781 String antennaName ; 744 1782 String stationName ; 745 1783 if ( pos != String::npos ) { 746 hAntennaName = hAntennaName.substr( pos+2 ) ; 747 } 748 pos = hAntennaName.find( "@" ) ; 1784 stationName = hAntName.substr( 0, pos ) ; 1785 hAntName = hAntName.substr( pos+2 ) ; 1786 } 1787 pos = hAntName.find( "@" ) ; 749 1788 if ( pos != String::npos ) { 750 antennaName = hAnt ennaName.substr( 0, pos ) ;751 stationName = hAnt ennaName.substr( pos+1 ) ;1789 antennaName = hAntName.substr( 0, pos ) ; 1790 stationName = hAntName.substr( pos+1 ) ; 752 1791 } 753 1792 else { 754 antennaName = hAntennaName ; 755 stationName = hAntennaName ; 756 } 757 // os_ << "antennaName = " << antennaName << LogIO::POST ; 758 // os_ << "stationName = " << stationName << LogIO::POST ; 1793 antennaName = hAntName ; 1794 } 1795 Vector<Double> antpos = keys.asArrayDouble( "AntennaPosition" ) ; 759 1796 760 msAntCols.name().put( 0, antennaName ) ; 761 msAntCols.station().put( 0, stationName ) ; 762 763 // os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ; 1797 String mount, atype ; 1798 Double diameter ; 1799 antennaProperty( antennaName, mount, atype, diameter ) ; 764 1800 765 msAntCols.position().put( 0, header_.antennaposition ) ; 766 767 // MOUNT is set to "ALT-AZ" 768 msAntCols.mount().put( 0, "ALT-AZ" ) ; 769 770 Double diameter = getDishDiameter( antennaName ) ; 771 msAntCols.dishDiameterQuant().put( 0, Quantity( diameter, "m" ) ) ; 772 773 // double endSec = mathutil::gettimeofday_sec() ; 774 // os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1801 TableRow tr( anttab ) ; 1802 TableRecord &r = tr.record() ; 1803 RecordFieldPtr<String> nameRF( r, "NAME" ) ; 1804 RecordFieldPtr<String> stationRF( r, "STATION" ) ; 1805 RecordFieldPtr<String> mountRF( r, "NAME" ) ; 1806 RecordFieldPtr<String> typeRF( r, "TYPE" ) ; 1807 RecordFieldPtr<Double> dishDiameterRF( r, "DISH_DIAMETER" ) ; 1808 RecordFieldPtr< Vector<Double> > positionRF( r, "POSITION" ) ; 1809 *nameRF = antennaName ; 1810 *mountRF = mount ; 1811 *typeRF = atype ; 1812 *dishDiameterRF = diameter ; 1813 *positionRF = antpos ; 1814 1815 tr.put( 0 ) ; 1816 1817 //double endSec = mathutil::gettimeofday_sec() ; 1818 //os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 775 1819 } 1820 1821 // void MSWriter::fillAntenna() 1822 // { 1823 // // double startSec = mathutil::gettimeofday_sec() ; 1824 // // os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ; 1825 1826 // // only 1 row 1827 // mstable_->antenna().addRow( 1, True ) ; 1828 // MSAntennaColumns msAntCols( mstable_->antenna() ) ; 1829 1830 // String hAntennaName = header_.antennaname ; 1831 // String::size_type pos = hAntennaName.find( "//" ) ; 1832 // String antennaName ; 1833 // String stationName ; 1834 // if ( pos != String::npos ) { 1835 // hAntennaName = hAntennaName.substr( pos+2 ) ; 1836 // } 1837 // pos = hAntennaName.find( "@" ) ; 1838 // if ( pos != String::npos ) { 1839 // antennaName = hAntennaName.substr( 0, pos ) ; 1840 // stationName = hAntennaName.substr( pos+1 ) ; 1841 // } 1842 // else { 1843 // antennaName = hAntennaName ; 1844 // stationName = hAntennaName ; 1845 // } 1846 // // os_ << "antennaName = " << antennaName << LogIO::POST ; 1847 // // os_ << "stationName = " << stationName << LogIO::POST ; 1848 1849 // msAntCols.name().put( 0, antennaName ) ; 1850 // msAntCols.station().put( 0, stationName ) ; 1851 1852 // // os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ; 1853 1854 // msAntCols.position().put( 0, header_.antennaposition ) ; 1855 1856 // // MOUNT is set to "ALT-AZ" 1857 // msAntCols.mount().put( 0, "ALT-AZ" ) ; 1858 1859 // Double diameter = getDishDiameter( antennaName ) ; 1860 // msAntCols.dishDiameterQuant().put( 0, Quantity( diameter, "m" ) ) ; 1861 1862 // // double endSec = mathutil::gettimeofday_sec() ; 1863 // // os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1864 // } 776 1865 777 1866 void MSWriter::fillProcessor() … … 836 1925 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ; 837 1926 Double srcVel = srcVelCol( 0 ) ; 1927 srcRec_.define( srcName, srcId ) ; 838 1928 839 1929 // NAME … … 990 2080 } 991 2081 992 void MSWriter::fillSysCal() 2082 void MSWriter::fillSysCal( map< Int,Vector<uInt> > &idRec, 2083 map< Int,Vector<uInt> > &rowRec ) 993 2084 { 994 //double startSec = mathutil::gettimeofday_sec() ;995 //os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ;996 997 // tcalIdRec_.print( cout ) ;2085 //double startSec = mathutil::gettimeofday_sec() ; 2086 //os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ; 2087 2088 //idRec.print( cout ) ; 998 2089 999 2090 // access to MS SYSCAL subtable … … 1017 2108 return ; 1018 2109 1019 nrow = tcalIdRec_.nfields() ; 2110 //nrow = idRec.nfields() ; 2111 nrow = idRec.size() ; 1020 2112 1021 2113 Double midTime ; … … 1052 2144 ROTableColumn beamnoCol( tab, "BEAMNO" ) ; 1053 2145 ROTableColumn ifnoCol( tab, "IFNO" ) ; 2146 map< Int,Vector<uInt> >::iterator itr0 = idRec.begin() ; 2147 map< Int,Vector<uInt> >::iterator itr1 = rowRec.begin() ; 1054 2148 for ( uInt irow = 0 ; irow < nrow ; irow++ ) { 1055 2149 // double t1 = mathutil::gettimeofday_sec() ; 1056 Vector<uInt> ids = tcalIdRec_.asArrayuInt( irow ) ; 2150 Vector<uInt> ids = itr0->second ; 2151 itr0++ ; 1057 2152 // os_ << "ids = " << ids << LogIO::POST ; 1058 2153 uInt npol = ids.size() ; 1059 Vector<uInt> rows = tcalRowRec_.asArrayuInt( irow ) ; 2154 Vector<uInt> rows = itr1->second ; 2155 itr1++ ; 1060 2156 // os_ << "rows = " << rows << LogIO::POST ; 1061 2157 Vector<Double> atime( rows.nelements() ) ; … … 1153 2249 if ( tcalSpec_ ) { 1154 2250 // put TCAL_SPECTRUM 1155 //*tcalspRF = tcal ;1156 2251 tcalspRF.define( tcal ) ; 1157 2252 // set TCAL (mean of TCAL_SPECTRUM) … … 1161 2256 } 1162 2257 // put TCAL 1163 *tcalRF = tcalMean;2258 tcalRF.define( tcalMean ) ; 1164 2259 } 1165 2260 else { 1166 2261 // put TCAL 1167 *tcalRF = tcal;2262 tcalRF.define( tcal ) ; 1168 2263 } 1169 2264 1170 2265 if ( tsysSpec_ ) { 1171 2266 // put TSYS_SPECTRUM 1172 //*tsysspRF = tsys ;1173 2267 tsysspRF.define( tsys ) ; 1174 2268 // set TSYS (mean of TSYS_SPECTRUM) … … 1178 2272 } 1179 2273 // put TSYS 1180 *tsysRF = tsysMean;2274 tsysRF.define( tsysMean ) ; 1181 2275 } 1182 2276 else { 1183 2277 // put TSYS 1184 *tsysRF = tsys;2278 tsysRF.define( tsys ) ; 1185 2279 } 1186 2280 … … 1193 2287 } 1194 2288 1195 //double endSec = mathutil::gettimeofday_sec() ;1196 //os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;2289 //double endSec = mathutil::gettimeofday_sec() ; 2290 //os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1197 2291 } 1198 2292 -
trunk/src/MSWriter.h
r2258 r2291 23 23 #include <casa/Logging/LogIO.h> 24 24 #include <casa/Containers/Record.h> 25 #include <casa/Containers/RecordField.h> 25 26 26 27 #include <tables/Tables/Table.h> … … 39 40 namespace asap 40 41 { 42 class MSWriterUtils 43 { 44 protected: 45 template<class T> void putField( const String &name, 46 TableRecord &r, 47 T &val ) 48 { 49 RecordFieldPtr<T> rf( r, name ) ; 50 *rf = val ; 51 } 52 template<class T> void defineField( const String &name, 53 TableRecord &r, 54 T &val ) 55 { 56 RecordFieldPtr<T> rf( r, name ) ; 57 rf.define( val ) ; 58 } 59 }; 41 60 42 61 class MSWriter … … 65 84 void fillSource() ; 66 85 void fillWeather() ; 67 void fillSysCal() ; 86 void fillSysCal( std::map< casa::Int,casa::Vector<casa::uInt> > &idrec, 87 std::map< casa::Int,casa::Vector<casa::uInt> > &rowrec ) ; 88 // void fillSysCal( Record &idrec, Record &rowrec ) ; 89 // void fillSysCal() ; 68 90 69 91 // fill empty rows … … 86 108 void queryType( casa::Int type, casa::String &stype, casa::Bool &b, casa::Double &t, Double &l ) ; 87 109 casa::Double getDishDiameter( casa::String antname ) ; 110 void antennaProperty( casa::String &name, casa::String &mount, casa::String &type, casa::Double &diameter ) ; 88 111 89 112 // tool for HPC … … 111 134 casa::LogIO os_ ; 112 135 113 casa::Record tcalIdRec_ ; 114 casa::Record tcalRowRec_ ; 136 // casa::Record tcalIdRec_ ; 137 // casa::Record tcalRowRec_ ; 138 casa::Record srcRec_ ; 115 139 116 140 MSWriter(); -
trunk/src/TableTraverse.cpp
r2289 r2291 251 251 sizes[i] = typeManagers[i]->sizeOf(); 252 252 BaseColumn *baseCol = ROTableColumnBackDoor::baseColumnPtr(cols[i]); 253 PlainColumn *col = dynamic_cast <PlainColumn *>(baseCol); 254 if (col) { 255 const uInt gotRows = col->dataManagerColumn()-> 256 getBlockV(0, nRows, colValues[i]); 257 DebugAssert(gotRows == nRows, AipsError); 258 } else { 259 copyColumnData(colValues[i], typeManagers[i]->sizeOf(), nRows, baseCol); 260 } 253 copyColumnData(colValues[i], typeManagers[i]->sizeOf(), nRows, baseCol); 254 // PlainColumn *col = dynamic_cast <PlainColumn *>(baseCol); 255 // if (col) { 256 // const uInt gotRows = col->dataManagerColumn()-> 257 // getBlockV(0, nRows, colValues[i]); 258 // DebugAssert(gotRows == nRows, AipsError); 259 // } else { 260 // copyColumnData(colValues[i], typeManagers[i]->sizeOf(), nRows, baseCol); 261 // } 261 262 } 262 263
Note:
See TracChangeset
for help on using the changeset viewer.