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