Changeset 2915 for trunk/src/STCalTsys.cpp
- Timestamp:
- 04/02/14 18:38:27 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/STCalTsys.cpp
r2786 r2915 18 18 #include "STSelector.h" 19 19 #include "STCalTsys.h" 20 #include "RowAccumulator.h"21 #include "STIdxIter.h"22 20 #include "STDefs.h" 23 21 #include <atnf/PKSIO/SrcType.h> … … 28 26 namespace asap { 29 27 STCalTsys::STCalTsys(CountedPtr<Scantable> &s, vector<int> &iflist) 30 : STCalibration(s ),28 : STCalibration(s, "TSYS"), 31 29 iflist_(iflist) 32 30 { … … 59 57 } 60 58 61 void STCalTsys::fillCalTable() 59 void STCalTsys::appenddata(uInt scanno, uInt cycleno, 60 uInt beamno, uInt ifno, uInt polno, 61 uInt freqid, Double time, Float elevation, 62 Vector<Float> any_data) 62 63 { 63 RowAccumulator acc(W_TINT); 64 65 vector<string> cols(3); 66 cols[0] = "IFNO"; 67 cols[1] = "POLNO"; 68 cols[2] = "BEAMNO"; 69 STIdxIterAcc iter(scantable_, cols); 70 71 ROScalarColumn<Double> *tcol = new ROScalarColumn<Double>(scantable_->table(), "TIME"); 72 Vector<Double> timeSec = tcol->getColumn() * 86400.0; 73 tcol->attach(scantable_->table(), "INTERVAL"); 74 Vector<Double> intervalSec = tcol->getColumn(); 75 delete tcol; 76 ROScalarColumn<Float> *ecol = new ROScalarColumn<Float>(scantable_->table(), "ELEVATION"); 77 Vector<Float> elevation = ecol->getColumn(); 78 delete ecol; 79 80 ROArrayColumn<Float> specCol(scantable_->table(), "TSYS"); 81 ROArrayColumn<uChar> flagCol(scantable_->table(), "FLAGTRA"); 82 ROScalarColumn<uInt> freqidCol(scantable_->table(), "FREQ_ID"); 83 84 // dummy Tsys: the following process doesn't need Tsys but RowAccumulator 85 // requires to set it with spectral data 86 Vector<Float> tsys(1, 1.0); 87 88 Double timeCen = 0.0; 89 Float elCen = 0.0; 90 uInt count = 0; 91 92 while(!iter.pastEnd()) { 93 Vector<uInt> rows = iter.getRows(SHARE); 94 Vector<uInt> current = iter.current(); 95 //os_ << "current=" << current << LogIO::POST; 96 uInt len = rows.nelements(); 97 if (len == 0) { 98 iter.next(); 99 continue; 100 } 101 else if (len == 1) { 102 STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_)); 103 uInt irow = rows[0]; 104 p->appenddata(0, 0, current[2], current[0], current[1], 105 freqidCol(irow), timeSec[irow], elevation[irow], specCol(irow)); 106 iter.next(); 107 continue; 108 } 109 110 uInt nchan = scantable_->nchan(scantable_->getIF(rows[0])); 111 Vector<uChar> flag(nchan); 112 Vector<Bool> bflag(nchan); 113 Vector<Float> spec(nchan); 114 115 Vector<Double> timeSep(len); 116 for (uInt i = 0; i < len-1; i++) { 117 timeSep[i] = timeSec[rows[i+1]] - timeSec[rows[i]] ; 118 } 119 Double tMedian = median(timeSep(IPosition(1,0), IPosition(1,len-2))); 120 timeSep[len-1] = tMedian * 10000.0 ; // any large value 121 122 uInt irow ; 123 uInt jrow ; 124 for (uInt i = 0; i < len; i++) { 125 //os_ << "start row " << rows[i] << LogIO::POST; 126 irow = rows[i]; 127 jrow = (i < len-1) ? rows[i+1] : rows[i]; 128 // accumulate data 129 flagCol.get(irow, flag); 130 convertArray(bflag, flag); 131 specCol.get(irow, spec); 132 if ( !allEQ(bflag,True) ) 133 acc.add( spec, !bflag, tsys, intervalSec[irow], timeSec[irow] ) ; 134 timeCen += timeSec[irow]; 135 elCen += elevation[irow]; 136 count++; 137 138 // check time gap 139 double gap = 2.0 * timeSep[i] / (intervalSec[jrow] + intervalSec[irow]); 140 if ( gap > 5.0 ) { 141 if ( acc.state() ) { 142 acc.replaceNaN() ; 143 // const Vector<Bool> &msk = acc.getMask(); 144 // convertArray(flag, !msk); 145 // for (uInt k = 0; k < nchan; ++k) { 146 // uChar userFlag = 1 << 7; 147 // if (msk[k]==True) userFlag = 0 << 7; 148 // flag(k) = userFlag; 149 // } 150 timeCen /= (Double)count * 86400.0; // sec->day 151 elCen /= (Float)count; 152 STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_)); 153 p->appenddata(0, 0, current[2], current[0], current[1], 154 freqidCol(irow), timeCen, elCen, acc.getSpectrum()); 155 } 156 acc.reset() ; 157 timeCen = 0.0; 158 elCen = 0.0; 159 count = 0; 160 } 161 } 162 163 iter.next() ; 164 //os_ << "end " << current << LogIO::POST; 165 } 64 STCalTsysTable *p = dynamic_cast<STCalTsysTable *>(&(*applytable_)); 65 p->appenddata(scanno, cycleno, beamno, ifno, polno, 66 freqid, time, elevation, any_data); 166 67 } 167 68
Note: See TracChangeset
for help on using the changeset viewer.