Changeset 2915 for trunk/src/STCalibration.cpp
- Timestamp:
- 04/02/14 18:38:27 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/STCalibration.cpp
r2786 r2915 13 13 #include "STCalibration.h" 14 14 #include "STCalSkyTable.h" 15 #include "RowAccumulator.h" 16 #include "STIdxIter.h" 15 17 16 18 using namespace casa; 17 19 18 20 namespace asap { 19 STCalibration::STCalibration(CountedPtr<Scantable> &s) 20 : scantable_(s) 21 STCalibration::STCalibration(CountedPtr<Scantable> &s, const String target_column) 22 : scantable_(s), 23 target_column_(target_column) 21 24 { 22 25 } … … 33 36 } 34 37 38 void STCalibration::fillCalTable() 39 { 40 RowAccumulator acc(W_TINT); 41 42 vector<string> cols(3); 43 cols[0] = "IFNO"; 44 cols[1] = "POLNO"; 45 cols[2] = "BEAMNO"; 46 STIdxIterAcc iter(scantable_, cols); 47 48 ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(scantable_->table(), "TIME"); 49 Vector<Double> timeSec = tcol->getColumn() * 86400.0; 50 tcol->attach(scantable_->table(), "INTERVAL"); 51 Vector<Double> intervalSec = tcol->getColumn(); 52 delete tcol; 53 ROScalarColumn<Float> *ecol = new ROScalarColumn<Float>(scantable_->table(), "ELEVATION"); 54 Vector<Float> elevation = ecol->getColumn(); 55 delete ecol; 56 57 ROArrayColumn<Float> specCol(scantable_->table(), target_column_); 58 ROArrayColumn<uChar> flagCol(scantable_->table(), "FLAGTRA"); 59 ROScalarColumn<uInt> freqidCol(scantable_->table(), "FREQ_ID"); 60 61 // dummy Tsys: the following process doesn't need Tsys but RowAccumulator 62 // requires to set it with spectral data 63 Vector<Float> tsys(1, 1.0); 64 65 Double timeCen = 0.0; 66 Float elCen = 0.0; 67 uInt count = 0; 68 69 while(!iter.pastEnd()) { 70 Vector<uInt> rows = iter.getRows(SHARE); 71 Vector<uInt> current = iter.current(); 72 //os_ << "current=" << current << LogIO::POST; 73 uInt len = rows.nelements(); 74 if (len == 0) { 75 iter.next(); 76 continue; 77 } 78 else if (len == 1) { 79 uInt irow = rows[0]; 80 appenddata(0, 0, current[2], current[0], current[1], 81 freqidCol(irow), timeSec[irow], elevation[irow], specCol(irow)); 82 iter.next(); 83 continue; 84 } 85 86 uInt nchan = scantable_->nchan(scantable_->getIF(rows[0])); 87 Vector<uChar> flag(nchan); 88 Vector<Bool> bflag(nchan); 89 Vector<Float> spec(nchan); 90 91 Vector<Double> timeSep(len); 92 for (uInt i = 0; i < len-1; i++) { 93 timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ; 94 } 95 Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2))); 96 timeSep[len-1] = tMedian * 10000.0 ; // any large value 97 98 uInt irow ; 99 uInt jrow ; 100 for (uInt i = 0; i < len; i++) { 101 //os_ << "start row " << rows[i] << LogIO::POST; 102 irow = rows[i]; 103 jrow = (i < len-1) ? rows[i+1] : rows[i]; 104 // accumulate data 105 flagCol.get(irow, flag); 106 convertArray(bflag, flag); 107 specCol.get(irow, spec); 108 if ( !allEQ(bflag,True) ) 109 acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ; 110 timeCen += timeSec[irow]; 111 elCen += elevation[irow]; 112 count++; 113 114 // check time gap 115 double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]); 116 if ( gap > 1.1 ) { 117 if ( acc.state() ) { 118 acc.replaceNaN() ; 119 // const Vector<Bool> &msk = acc.getMask(); 120 // convertArray(flag, !msk); 121 // for (uInt k = 0; k < nchan; ++k) { 122 // uChar userFlag = 1 << 7; 123 // if (msk[k]==True) userFlag = 0 << 7; 124 // flag(k) = userFlag; 125 // } 126 timeCen /= (Double)count * 86400.0; // sec->day 127 elCen /= (Float)count; 128 appenddata(0, 0, current[2], current[0], current[1], 129 freqidCol(irow), timeCen, elCen, acc.getSpectrum()); 130 } 131 acc.reset() ; 132 timeCen = 0.0; 133 elCen = 0.0; 134 count = 0; 135 } 136 } 137 138 iter.next() ; 139 //os_ << "end " << current << LogIO::POST; 140 } 35 141 } 142 143 }
Note: See TracChangeset
for help on using the changeset viewer.