- Timestamp:
- 06/13/12 18:08:22 (13 years ago)
- Location:
- branches/hpc33/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/hpc33/src/STMath.cpp
r2564 r2566 3884 3884 out->attach() ; 3885 3885 out->table().addRow( s->nrow() ) ; 3886 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ;3886 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ; 3887 3887 3888 3888 // process each on scan … … 4034 4034 vector<bool> masks = s->getMask( 0 ) ; 4035 4035 CountedPtr<Scantable> ssig, sref ; 4036 CountedPtr<Scantable> out ; 4036 //CountedPtr<Scantable> out ; 4037 bool insitu = insitu_ ; 4038 insitu_ = False ; 4039 CountedPtr<Scantable> out = getScantable( s, true ) ; 4040 insitu_ = insitu ; 4037 4041 4038 4042 if ( antname.find( "APEX" ) != string::npos ) { 4039 4043 // APEX calibration 4040 4044 // sky scan 4041 STSelector sel = STSelector() ; 4042 types.push_back( SrcType::FLOSKY ) ; 4043 sel.setTypes( types ) ; 4044 s->setSelection( sel ) ; 4045 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ; 4046 CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ; 4047 s->unsetSelection() ; 4048 sel.reset() ; 4049 types.clear() ; 4050 types.push_back( SrcType::FHISKY ) ; 4051 sel.setTypes( types ) ; 4052 s->setSelection( sel ) ; 4053 tmp.clear() ; 4054 tmp.push_back( getScantable( s, false ) ) ; 4055 CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ; 4056 s->unsetSelection() ; 4057 sel.reset() ; 4058 types.clear() ; 4045 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOSKY ) ; 4046 out->attach() ; 4047 CountedPtr<Scantable> askylo = averageWithinSession( out, 4048 masks, 4049 "TINT" ) ; 4050 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHISKY ) ; 4051 out->attach() ; 4052 CountedPtr<Scantable> askyhi = averageWithinSession( out, 4053 masks, 4054 "TINT" ) ; 4059 4055 4060 4056 // hot scan 4061 types.push_back( SrcType::FLOHOT ) ; 4062 sel.setTypes( types ) ; 4063 s->setSelection( sel ) ; 4064 tmp.clear() ; 4065 tmp.push_back( getScantable( s, false ) ) ; 4066 CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ; 4067 s->unsetSelection() ; 4068 sel.reset() ; 4069 types.clear() ; 4070 types.push_back( SrcType::FHIHOT ) ; 4071 sel.setTypes( types ) ; 4072 s->setSelection( sel ) ; 4073 tmp.clear() ; 4074 tmp.push_back( getScantable( s, false ) ) ; 4075 CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ; 4076 s->unsetSelection() ; 4077 sel.reset() ; 4078 types.clear() ; 4057 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOHOT ) ; 4058 out->attach() ; 4059 CountedPtr<Scantable> ahotlo = averageWithinSession( out, 4060 masks, 4061 "TINT" ) ; 4062 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIHOT ) ; 4063 out->attach() ; 4064 CountedPtr<Scantable> ahothi = averageWithinSession( out, 4065 masks, 4066 "TINT" ) ; 4079 4067 4080 4068 // cold scan 4081 4069 CountedPtr<Scantable> acoldlo, acoldhi ; 4082 // types.push_back( SrcType::FLOCOLD ) ; 4083 // sel.setTypes( types ) ; 4084 // s->setSelection( sel ) ; 4085 // tmp.clear() ; 4086 // tmp.push_back( getScantable( s, false ) ) ; 4087 // CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ; 4088 // s->unsetSelection() ; 4089 // sel.reset() ; 4090 // types.clear() ; 4091 // types.push_back( SrcType::FHICOLD ) ; 4092 // sel.setTypes( types ) ; 4093 // s->setSelection( sel ) ; 4094 // tmp.clear() ; 4095 // tmp.push_back( getScantable( s, false ) ) ; 4096 // CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ; 4097 // s->unsetSelection() ; 4098 // sel.reset() ; 4099 // types.clear() ; 4070 // out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOCOLD ) ; 4071 // out->attach() ; 4072 // CountedPtr<Scantable> acoldlo = averageWithinSession( out, 4073 // masks, 4074 // "TINT" ) ; 4075 // out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHICOLD ) ; 4076 // out->attach() ; 4077 // CountedPtr<Scantable> acoldhi = averageWithinSession( out, 4078 // masks, 4079 // "TINT" ) ; 4100 4080 4101 4081 // ref scan 4102 bool insitu = insitu_ ;4103 4082 insitu_ = false ; 4104 4083 sref = getScantable( s, true ) ; 4084 CountedPtr<Scantable> rref = getScantable( s, true ) ; 4105 4085 insitu_ = insitu ; 4106 types.push_back( SrcType::FSLO ) ; 4107 sel.setTypes( types ) ; 4108 s->setSelection( sel ) ; 4109 TableCopy::copyRows( sref->table(), s->table() ) ; 4110 s->unsetSelection() ; 4111 sel.reset() ; 4112 types.clear() ; 4086 rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSLO ) ; 4087 rref->attach() ; 4088 copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ; 4113 4089 4114 4090 // sig scan 4115 4091 insitu_ = false ; 4116 4092 ssig = getScantable( s, true ) ; 4093 CountedPtr<Scantable> rsig = getScantable( s, true ) ; 4117 4094 insitu_ = insitu ; 4118 types.push_back( SrcType::FSHI ) ; 4119 sel.setTypes( types ) ; 4120 s->setSelection( sel ) ; 4121 TableCopy::copyRows( ssig->table(), s->table() ) ; 4122 s->unsetSelection() ; 4123 sel.reset() ; 4124 types.clear() ; 4095 rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSHI ) ; 4096 rsig->attach() ; 4097 copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ; 4125 4098 4126 4099 if ( apexcalmode == 0 ) { 4127 // APEX fs data without off scan 4128 // process each sig and ref scan 4129 ArrayColumn<Float> tsysCollo ; 4130 tsysCollo.attach( ssig->table(), "TSYS" ) ; 4131 ArrayColumn<Float> tsysColhi ; 4132 tsysColhi.attach( sref->table(), "TSYS" ) ; 4133 for ( int i = 0 ; i < ssig->nrow() ; i++ ) { 4134 vector< CountedPtr<Scantable> > sky( 2 ) ; 4135 sky[0] = askylo ; 4136 sky[1] = askyhi ; 4137 vector< CountedPtr<Scantable> > hot( 2 ) ; 4138 hot[0] = ahotlo ; 4139 hot[1] = ahothi ; 4140 vector< CountedPtr<Scantable> > cold( 2 ) ; 4141 //cold[0] = acoldlo ; 4142 //cold[1] = acoldhi ; 4143 vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ; 4144 ssig->setSpectrum( sp, i ) ; 4145 string reftime = ssig->getTime( i ) ; 4146 vector<int> ii( 1, ssig->getIF( i ) ) ; 4147 vector<int> ib( 1, ssig->getBeam( i ) ) ; 4148 vector<int> ip( 1, ssig->getPol( i ) ) ; 4149 sel.setIFs( ii ) ; 4150 sel.setBeams( ib ) ; 4151 sel.setPolarizations( ip ) ; 4152 askylo->setSelection( sel ) ; 4153 vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ; 4154 const Vector<Float> Vtsyslo( sptsys ) ; 4155 tsysCollo.put( i, Vtsyslo ) ; 4156 askylo->unsetSelection() ; 4100 // using STIdxIterAcc 4101 vector<string> cols( 3 ) ; 4102 cols[0] = "BEAMNO" ; 4103 cols[1] = "POLNO" ; 4104 cols[2] = "IFNO" ; 4105 STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ; 4106 STSelector sel ; 4107 vector< CountedPtr<Scantable> > on( 2 ) ; 4108 on[0] = rsig ; 4109 on[1] = rref ; 4110 vector< CountedPtr<Scantable> > sky( 2 ) ; 4111 sky[0] = askylo ; 4112 sky[1] = askyhi ; 4113 vector< CountedPtr<Scantable> > hot( 2 ) ; 4114 hot[0] = ahotlo ; 4115 hot[1] = ahothi ; 4116 vector< CountedPtr<Scantable> > cold( 2 ) ; 4117 while ( !iter->pastEnd() ) { 4118 Vector<uInt> ids = iter->current() ; 4119 stringstream ss ; 4120 ss << "SELECT FROM $1 WHERE " 4121 << "BEAMNO==" << ids[0] << "&&" 4122 << "POLNO==" << ids[1] << "&&" 4123 << "IFNO==" << ids[2] ; 4124 //cout << "TaQL string: " << ss.str() << endl ; 4125 sel.setTaQL( ss.str() ) ; 4126 sky[0]->setSelection( sel ) ; 4127 sky[1]->setSelection( sel ) ; 4128 hot[0]->setSelection( sel ) ; 4129 hot[1]->setSelection( sel ) ; 4130 Vector<uInt> rows = iter->getRows( SHARE ) ; 4131 calibrateAPEXFS( ssig, sref, on, sky, hot, cold, rows ) ; 4132 sky[0]->unsetSelection() ; 4133 sky[1]->unsetSelection() ; 4134 hot[0]->unsetSelection() ; 4135 hot[1]->unsetSelection() ; 4157 4136 sel.reset() ; 4158 sky[0] = askyhi ; 4159 sky[1] = askylo ; 4160 hot[0] = ahothi ; 4161 hot[1] = ahotlo ; 4162 cold[0] = acoldhi ; 4163 cold[1] = acoldlo ; 4164 sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ; 4165 sref->setSpectrum( sp, i ) ; 4166 reftime = sref->getTime( i ) ; 4167 ii[0] = sref->getIF( i ) ; 4168 ib[0] = sref->getBeam( i ) ; 4169 ip[0] = sref->getPol( i ) ; 4170 sel.setIFs( ii ) ; 4171 sel.setBeams( ib ) ; 4172 sel.setPolarizations( ip ) ; 4173 askyhi->setSelection( sel ) ; 4174 sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ; 4175 const Vector<Float> Vtsyshi( sptsys ) ; 4176 tsysColhi.put( i, Vtsyshi ) ; 4177 askyhi->unsetSelection() ; 4178 sel.reset() ; 4179 } 4137 iter->next() ; 4138 } 4139 delete iter ; 4140 4180 4141 } 4181 4142 else if ( apexcalmode == 1 ) { 4182 4143 // APEX fs data with off scan 4183 4144 // off scan 4184 types.push_back( SrcType::FLOOFF ) ; 4185 sel.setTypes( types ) ; 4186 s->setSelection( sel ) ; 4187 tmp.clear() ; 4188 tmp.push_back( getScantable( s, false ) ) ; 4189 CountedPtr<Scantable> aofflo = average( tmp, masks, "TINT", "SCAN" ) ; 4190 s->unsetSelection() ; 4191 sel.reset() ; 4192 types.clear() ; 4193 types.push_back( SrcType::FHIOFF ) ; 4194 sel.setTypes( types ) ; 4195 s->setSelection( sel ) ; 4196 tmp.clear() ; 4197 tmp.push_back( getScantable( s, false ) ) ; 4198 CountedPtr<Scantable> aoffhi = average( tmp, masks, "TINT", "SCAN" ) ; 4199 s->unsetSelection() ; 4200 sel.reset() ; 4201 types.clear() ; 4145 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOOFF ) ; 4146 out->attach() ; 4147 CountedPtr<Scantable> aofflo = averageWithinSession( out, 4148 masks, 4149 "TINT" ) ; 4150 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIOFF ) ; 4151 out->attach() ; 4152 CountedPtr<Scantable> aoffhi = averageWithinSession( out, 4153 masks, 4154 "TINT" ) ; 4202 4155 4203 4156 // process each sig and ref scan 4204 ArrayColumn<Float> tsysCollo ; 4205 tsysCollo.attach( ssig->table(), "TSYS" ) ; 4206 ArrayColumn<Float> tsysColhi ; 4207 tsysColhi.attach( sref->table(), "TSYS" ) ; 4208 for ( int i = 0 ; i < ssig->nrow() ; i++ ) { 4209 vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ; 4210 ssig->setSpectrum( sp, i ) ; 4211 sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ; 4212 string reftime = ssig->getTime( i ) ; 4213 vector<int> ii( 1, ssig->getIF( i ) ) ; 4214 vector<int> ib( 1, ssig->getBeam( i ) ) ; 4215 vector<int> ip( 1, ssig->getPol( i ) ) ; 4216 sel.setIFs( ii ) ; 4217 sel.setBeams( ib ) ; 4218 sel.setPolarizations( ip ) ; 4157 // STSelector sel ; 4158 vector<string> cols( 3 ) ; 4159 cols[0] = "BEAMNO" ; 4160 cols[1] = "POLNO" ; 4161 cols[2] = "IFNO" ; 4162 STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ; 4163 STSelector sel ; 4164 while ( !iter->pastEnd() ) { 4165 Vector<uInt> ids = iter->current() ; 4166 stringstream ss ; 4167 ss << "SELECT FROM $1 WHERE " 4168 << "BEAMNO==" << ids[0] << "&&" 4169 << "POLNO==" << ids[1] << "&&" 4170 << "IFNO==" << ids[2] ; 4171 //cout << "TaQL string: " << ss.str() << endl ; 4172 sel.setTaQL( ss.str() ) ; 4173 aofflo->setSelection( sel ) ; 4174 ahotlo->setSelection( sel ) ; 4219 4175 askylo->setSelection( sel ) ; 4220 vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ; 4221 const Vector<Float> Vtsyslo( sptsys ) ; 4222 tsysCollo.put( i, Vtsyslo ) ; 4176 Vector<uInt> rows = iter->getRows( SHARE ) ; 4177 calibrateCW( ssig, rsig, aofflo, askylo, ahotlo, acoldlo, rows, antname ) ; 4178 aofflo->unsetSelection() ; 4179 ahotlo->unsetSelection() ; 4223 4180 askylo->unsetSelection() ; 4224 4181 sel.reset() ; 4225 sref->setSpectrum( sp, i ) ; 4226 reftime = sref->getTime( i ) ; 4227 ii[0] = sref->getIF( i ) ; 4228 ib[0] = sref->getBeam( i ) ; 4229 ip[0] = sref->getPol( i ) ; 4230 sel.setIFs( ii ) ; 4231 sel.setBeams( ib ) ; 4232 sel.setPolarizations( ip ) ; 4233 askyhi->setSelection( sel ) ; 4234 sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ; 4235 const Vector<Float> Vtsyshi( sptsys ) ; 4236 tsysColhi.put( i, Vtsyshi ) ; 4182 iter->next() ; 4183 } 4184 delete iter ; 4185 iter = new STIdxIterAcc( sref, cols ) ; 4186 while ( !iter->pastEnd() ) { 4187 Vector<uInt> ids = iter->current() ; 4188 stringstream ss ; 4189 ss << "SELECT FROM $1 WHERE " 4190 << "BEAMNO==" << ids[0] << "&&" 4191 << "POLNO==" << ids[1] << "&&" 4192 << "IFNO==" << ids[2] ; 4193 //cout << "TaQL string: " << ss.str() << endl ; 4194 sel.setTaQL( ss.str() ) ; 4195 aoffhi->setSelection( sel ) ; 4196 ahothi->setSelection( sel ) ; 4197 askyhi->setSelection( sel ) ; 4198 Vector<uInt> rows = iter->getRows( SHARE ) ; 4199 calibrateCW( sref, rref, aoffhi, askyhi, ahothi, acoldhi, rows, antname ) ; 4200 aoffhi->unsetSelection() ; 4201 ahothi->unsetSelection() ; 4237 4202 askyhi->unsetSelection() ; 4238 4203 sel.reset() ; 4239 } 4204 iter->next() ; 4205 } 4206 delete iter ; 4240 4207 } 4241 4208 } … … 4243 4210 // non-APEX fs data 4244 4211 // sky scan 4212 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ; 4213 out->attach() ; 4214 CountedPtr<Scantable> asky = averageWithinSession( out, 4215 masks, 4216 "TINT" ) ; 4245 4217 STSelector sel = STSelector() ; 4246 types.push_back( SrcType::SKY ) ; 4247 sel.setTypes( types ) ; 4248 s->setSelection( sel ) ; 4249 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ; 4250 CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ; 4251 s->unsetSelection() ; 4252 sel.reset() ; 4253 types.clear() ; 4254 4218 4255 4219 // hot scan 4256 types.push_back( SrcType::HOT ) ; 4257 sel.setTypes( types ) ; 4258 s->setSelection( sel ) ; 4259 tmp.clear() ; 4260 tmp.push_back( getScantable( s, false ) ) ; 4261 CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ; 4262 s->unsetSelection() ; 4263 sel.reset() ; 4264 types.clear() ; 4220 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ; 4221 out->attach() ; 4222 CountedPtr<Scantable> ahot = averageWithinSession( out, 4223 masks, 4224 "TINT" ) ; 4265 4225 4266 4226 // cold scan 4267 4227 CountedPtr<Scantable> acold ; 4268 // types.push_back( SrcType::COLD ) ; 4269 // sel.setTypes( types ) ; 4270 // s->setSelection( sel ) ; 4271 // tmp.clear() ; 4272 // tmp.push_back( getScantable( s, false ) ) ; 4273 // CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCAN" ) ; 4274 // s->unsetSelection() ; 4275 // sel.reset() ; 4276 // types.clear() ; 4228 // out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ; 4229 // out->attach() ; 4230 // CountedPtr<Scantable> acold = averageWithinSession( out, 4231 // masks, 4232 // "TINT" ) ; 4277 4233 4278 4234 // ref scan … … 4280 4236 insitu_ = false ; 4281 4237 sref = getScantable( s, true ) ; 4238 CountedPtr<Scantable> rref = getScantable( s, true ) ; 4282 4239 insitu_ = insitu ; 4283 types.push_back( SrcType::FSOFF ) ; 4284 sel.setTypes( types ) ; 4285 s->setSelection( sel ) ; 4286 TableCopy::copyRows( sref->table(), s->table() ) ; 4287 s->unsetSelection() ; 4288 sel.reset() ; 4289 types.clear() ; 4240 rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ; 4241 rref->attach() ; 4242 copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ; 4290 4243 4291 4244 // sig scan 4292 4245 insitu_ = false ; 4293 4246 ssig = getScantable( s, true ) ; 4247 CountedPtr<Scantable> rsig = getScantable( s, true ) ; 4294 4248 insitu_ = insitu ; 4295 types.push_back( SrcType::FSON ) ; 4296 sel.setTypes( types ) ; 4297 s->setSelection( sel ) ; 4298 TableCopy::copyRows( ssig->table(), s->table() ) ; 4299 s->unsetSelection() ; 4300 sel.reset() ; 4301 types.clear() ; 4249 rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ; 4250 rsig->attach() ; 4251 copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ; 4302 4252 4303 4253 // process each sig and ref scan 4304 ArrayColumn<Float> tsysColsig ; 4305 tsysColsig.attach( ssig->table(), "TSYS" ) ; 4306 ArrayColumn<Float> tsysColref ; 4307 tsysColref.attach( ssig->table(), "TSYS" ) ; 4308 for ( int i = 0 ; i < ssig->nrow() ; i++ ) { 4309 vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ; 4310 ssig->setSpectrum( sp, i ) ; 4311 string reftime = ssig->getTime( i ) ; 4312 vector<int> ii( 1, ssig->getIF( i ) ) ; 4313 vector<int> ib( 1, ssig->getBeam( i ) ) ; 4314 vector<int> ip( 1, ssig->getPol( i ) ) ; 4315 sel.setIFs( ii ) ; 4316 sel.setBeams( ib ) ; 4317 sel.setPolarizations( ip ) ; 4254 vector<string> cols( 3 ) ; 4255 cols[0] = "BEAMNO" ; 4256 cols[1] = "POLNO" ; 4257 cols[2] = "IFNO" ; 4258 STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ; 4259 while ( !iter->pastEnd() ) { 4260 Vector<uInt> ids = iter->current() ; 4261 stringstream ss ; 4262 ss << "SELECT FROM $1 WHERE " 4263 << "BEAMNO==" << ids[0] << "&&" 4264 << "POLNO==" << ids[1] << "&&" 4265 << "IFNO==" << ids[2] ; 4266 //cout << "TaQL string: " << ss.str() << endl ; 4267 sel.setTaQL( ss.str() ) ; 4268 ahot->setSelection( sel ) ; 4318 4269 asky->setSelection( sel ) ; 4319 vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ; 4320 const Vector<Float> Vtsys( sptsys ) ; 4321 tsysColsig.put( i, Vtsys ) ; 4270 Vector<uInt> rows = iter->getRows( SHARE ) ; 4271 // out should be an exact copy of s except that SPECTRA column is empty 4272 calibrateFS( ssig, sref, rsig, rref, asky, ahot, acold, rows ) ; 4273 ahot->unsetSelection() ; 4322 4274 asky->unsetSelection() ; 4323 4275 sel.reset() ; 4324 sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ; 4325 sref->setSpectrum( sp, i ) ; 4326 tsysColref.put( i, Vtsys ) ; 4327 } 4276 iter->next() ; 4277 } 4278 delete iter ; 4328 4279 } 4329 4280 … … 4331 4282 Table sigtab = ssig->table() ; 4332 4283 Table reftab = sref->table() ; 4333 ScalarColumn<uInt> sigifnoCol ;4334 ScalarColumn<uInt> refifnoCol ;4335 ScalarColumn<uInt> sigfidCol ;4336 4284 ScalarColumn<uInt> reffidCol ; 4337 4285 Int nchan = (Int)ssig->nchan() ; 4338 sigifnoCol.attach( sigtab, "IFNO" ) ;4339 refifnoCol.attach( reftab, "IFNO" ) ;4340 sigfidCol.attach( sigtab, "FREQ_ID" ) ;4341 4286 reffidCol.attach( reftab, "FREQ_ID" ) ; 4342 Vector<uInt> sfids ( sigfidCol.getColumn()) ;4343 Vector<uInt> rfids ( reffidCol.getColumn()) ;4287 Vector<uInt> sfids = ssig->mfreqidCol_.getColumn() ; 4288 Vector<uInt> rfids = sref->mfreqidCol_.getColumn() ; 4344 4289 vector<uInt> sfids_unique ; 4345 4290 vector<uInt> rfids_unique ; … … 4420 4365 } 4421 4366 4422 vector<float> STMath::getSpectrumFromTime( string reftime,4423 CountedPtr<Scantable>& s,4424 string mode )4425 {4426 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;4427 vector<float> sp ;4428 4429 if ( s->nrow() == 0 ) {4430 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;4431 return sp ;4432 }4433 else if ( s->nrow() == 1 ) {4434 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;4435 return s->getSpectrum( 0 ) ;4436 }4437 else {4438 vector<int> idx = getRowIdFromTime( reftime, s ) ;4439 if ( mode == "before" ) {4440 int id = -1 ;4441 if ( idx[0] != -1 ) {4442 id = idx[0] ;4443 }4444 else if ( idx[1] != -1 ) {4445 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;4446 id = idx[1] ;4447 }4448 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4449 sp = s->getSpectrum( id ) ;4450 }4451 else if ( mode == "after" ) {4452 int id = -1 ;4453 if ( idx[1] != -1 ) {4454 id = idx[1] ;4455 }4456 else if ( idx[0] != -1 ) {4457 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;4458 id = idx[1] ;4459 }4460 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4461 sp = s->getSpectrum( id ) ;4462 }4463 else if ( mode == "nearest" ) {4464 int id = -1 ;4465 if ( idx[0] == -1 ) {4466 id = idx[1] ;4467 }4468 else if ( idx[1] == -1 ) {4469 id = idx[0] ;4470 }4471 else if ( idx[0] == idx[1] ) {4472 id = idx[0] ;4473 }4474 else {4475 //double t0 = getMJD( s->getTime( idx[0] ) ) ;4476 //double t1 = getMJD( s->getTime( idx[1] ) ) ;4477 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;4478 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;4479 double tref = getMJD( reftime ) ;4480 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {4481 id = idx[1] ;4482 }4483 else {4484 id = idx[0] ;4485 }4486 }4487 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4488 sp = s->getSpectrum( id ) ;4489 }4490 else if ( mode == "linear" ) {4491 if ( idx[0] == -1 ) {4492 // use after4493 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;4494 int id = idx[1] ;4495 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4496 sp = s->getSpectrum( id ) ;4497 }4498 else if ( idx[1] == -1 ) {4499 // use before4500 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;4501 int id = idx[0] ;4502 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4503 sp = s->getSpectrum( id ) ;4504 }4505 else if ( idx[0] == idx[1] ) {4506 // use before4507 //os << "No need to interporate." << LogIO::POST ;4508 int id = idx[0] ;4509 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4510 sp = s->getSpectrum( id ) ;4511 }4512 else {4513 // do interpolation4514 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;4515 //double t0 = getMJD( s->getTime( idx[0] ) ) ;4516 //double t1 = getMJD( s->getTime( idx[1] ) ) ;4517 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;4518 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;4519 double tref = getMJD( reftime ) ;4520 vector<float> sp0 = s->getSpectrum( idx[0] ) ;4521 vector<float> sp1 = s->getSpectrum( idx[1] ) ;4522 for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {4523 float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;4524 sp.push_back( v ) ;4525 }4526 }4527 }4528 else {4529 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;4530 }4531 return sp ;4532 }4533 }4534 4535 4367 Vector<Float> STMath::getSpectrumFromTime( double reftime, 4536 4368 const Vector<Double> &timeVec, … … 4645 4477 } 4646 4478 4647 double STMath::getMJD( string strtime ) 4648 { 4649 if ( strtime.find("/") == string::npos ) { 4650 // MJD time string 4651 return atof( strtime.c_str() ) ; 4652 } 4653 else { 4654 // string in YYYY/MM/DD/HH:MM:SS format 4655 uInt year = atoi( strtime.substr( 0, 4 ).c_str() ) ; 4656 uInt month = atoi( strtime.substr( 5, 2 ).c_str() ) ; 4657 uInt day = atoi( strtime.substr( 8, 2 ).c_str() ) ; 4658 uInt hour = atoi( strtime.substr( 11, 2 ).c_str() ) ; 4659 uInt minute = atoi( strtime.substr( 14, 2 ).c_str() ) ; 4660 uInt sec = atoi( strtime.substr( 17, 2 ).c_str() ) ; 4661 Time t( year, month, day, hour, minute, sec ) ; 4662 return t.modifiedJulianDay() ; 4663 } 4664 } 4665 4666 vector<int> STMath::getRowIdFromTime( string reftime, CountedPtr<Scantable> &s ) 4667 { 4668 double reft = getMJD( reftime ) ; 4479 vector<int> STMath::getRowIdFromTime( double reftime, const Vector<Double> &t ) 4480 { 4481 // double reft = reftime ; 4669 4482 double dtmin = 1.0e100 ; 4670 4483 double dtmax = -1.0e100 ; 4671 vector<double> dt ;4484 // vector<double> dt ; 4672 4485 int just_before = -1 ; 4673 4486 int just_after = -1 ; 4674 for ( int i = 0 ; i < s->nrow() ; i++ ) { 4675 dt.push_back( getMJD( s->getTime( i ) ) - reft ) ; 4676 } 4487 Vector<Double> dt = t - reftime ; 4677 4488 for ( unsigned int i = 0 ; i < dt.size() ; i++ ) { 4678 4489 if ( dt[i] > 0.0 ) { … … 4700 4511 } 4701 4512 4702 vector<int> v ;4703 v.push_back( just_before ) ;4704 v.push_back( just_after ) ;4705 4706 return v ;4707 }4708 4709 vector<int> STMath::getRowIdFromTime( double reftime, const Vector<Double> &t )4710 {4711 // double reft = reftime ;4712 double dtmin = 1.0e100 ;4713 double dtmax = -1.0e100 ;4714 // vector<double> dt ;4715 int just_before = -1 ;4716 int just_after = -1 ;4717 //cout << setprecision(24) << reft << endl ;4718 // ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;4719 // for ( int i = 0 ; i < s->nrow() ; i++ ) {4720 // cout << setprecision(24) << timeCol(i) << endl ;4721 // //dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;4722 // dt.push_back( timeCol(i) - reft ) ;4723 // }4724 Vector<Double> dt = t - reftime ;4725 for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {4726 if ( dt[i] > 0.0 ) {4727 // after reftime4728 if ( dt[i] < dtmin ) {4729 just_after = i ;4730 dtmin = dt[i] ;4731 }4732 }4733 else if ( dt[i] < 0.0 ) {4734 // before reftime4735 if ( dt[i] > dtmax ) {4736 just_before = i ;4737 dtmax = dt[i] ;4738 }4739 }4740 else {4741 // just a reftime4742 just_before = i ;4743 just_after = i ;4744 dtmax = 0 ;4745 dtmin = 0 ;4746 break ;4747 }4748 }4749 4750 4513 vector<int> v(2) ; 4751 4514 v[0] = just_before ; … … 4753 4516 4754 4517 return v ; 4755 }4756 4757 vector<float> STMath::getTcalFromTime( string reftime,4758 CountedPtr<Scantable>& s,4759 string mode )4760 {4761 LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ;4762 vector<float> tcal ;4763 STTcal tcalTable = s->tcal() ;4764 String time ;4765 Vector<Float> tcalval ;4766 if ( s->nrow() == 0 ) {4767 os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;4768 return tcal ;4769 }4770 else if ( s->nrow() == 1 ) {4771 uInt tcalid = s->getTcalId( 0 ) ;4772 //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;4773 tcalTable.getEntry( time, tcalval, tcalid ) ;4774 tcalval.tovector( tcal ) ;4775 return tcal ;4776 }4777 else {4778 vector<int> idx = getRowIdFromTime( reftime, s ) ;4779 if ( mode == "before" ) {4780 int id = -1 ;4781 if ( idx[0] != -1 ) {4782 id = idx[0] ;4783 }4784 else if ( idx[1] != -1 ) {4785 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;4786 id = idx[1] ;4787 }4788 uInt tcalid = s->getTcalId( id ) ;4789 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4790 tcalTable.getEntry( time, tcalval, tcalid ) ;4791 tcalval.tovector( tcal ) ;4792 }4793 else if ( mode == "after" ) {4794 int id = -1 ;4795 if ( idx[1] != -1 ) {4796 id = idx[1] ;4797 }4798 else if ( idx[0] != -1 ) {4799 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;4800 id = idx[1] ;4801 }4802 uInt tcalid = s->getTcalId( id ) ;4803 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4804 tcalTable.getEntry( time, tcalval, tcalid ) ;4805 tcalval.tovector( tcal ) ;4806 }4807 else if ( mode == "nearest" ) {4808 int id = -1 ;4809 if ( idx[0] == -1 ) {4810 id = idx[1] ;4811 }4812 else if ( idx[1] == -1 ) {4813 id = idx[0] ;4814 }4815 else if ( idx[0] == idx[1] ) {4816 id = idx[0] ;4817 }4818 else {4819 //double t0 = getMJD( s->getTime( idx[0] ) ) ;4820 //double t1 = getMJD( s->getTime( idx[1] ) ) ;4821 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;4822 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;4823 double tref = getMJD( reftime ) ;4824 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {4825 id = idx[1] ;4826 }4827 else {4828 id = idx[0] ;4829 }4830 }4831 uInt tcalid = s->getTcalId( id ) ;4832 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4833 tcalTable.getEntry( time, tcalval, tcalid ) ;4834 tcalval.tovector( tcal ) ;4835 }4836 else if ( mode == "linear" ) {4837 if ( idx[0] == -1 ) {4838 // use after4839 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;4840 int id = idx[1] ;4841 uInt tcalid = s->getTcalId( id ) ;4842 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4843 tcalTable.getEntry( time, tcalval, tcalid ) ;4844 tcalval.tovector( tcal ) ;4845 }4846 else if ( idx[1] == -1 ) {4847 // use before4848 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;4849 int id = idx[0] ;4850 uInt tcalid = s->getTcalId( id ) ;4851 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4852 tcalTable.getEntry( time, tcalval, tcalid ) ;4853 tcalval.tovector( tcal ) ;4854 }4855 else if ( idx[0] == idx[1] ) {4856 // use before4857 //os << "No need to interporate." << LogIO::POST ;4858 int id = idx[0] ;4859 uInt tcalid = s->getTcalId( id ) ;4860 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4861 tcalTable.getEntry( time, tcalval, tcalid ) ;4862 tcalval.tovector( tcal ) ;4863 }4864 else {4865 // do interpolation4866 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;4867 //double t0 = getMJD( s->getTime( idx[0] ) ) ;4868 //double t1 = getMJD( s->getTime( idx[1] ) ) ;4869 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;4870 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;4871 double tref = getMJD( reftime ) ;4872 vector<float> tcal0 ;4873 vector<float> tcal1 ;4874 uInt tcalid0 = s->getTcalId( idx[0] ) ;4875 uInt tcalid1 = s->getTcalId( idx[1] ) ;4876 tcalTable.getEntry( time, tcalval, tcalid0 ) ;4877 tcalval.tovector( tcal0 ) ;4878 tcalTable.getEntry( time, tcalval, tcalid1 ) ;4879 tcalval.tovector( tcal1 ) ;4880 for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) {4881 float v = ( tcal1[i] - tcal0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tcal0[i] ;4882 tcal.push_back( v ) ;4883 }4884 }4885 }4886 else {4887 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;4888 }4889 return tcal ;4890 }4891 4518 } 4892 4519 … … 5011 4638 } 5012 4639 5013 vector<float> STMath::getTsysFromTime( string reftime,5014 CountedPtr<Scantable>& s,5015 string mode )5016 {5017 LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ;5018 ArrayColumn<Float> tsysCol ;5019 tsysCol.attach( s->table(), "TSYS" ) ;5020 vector<float> tsys ;5021 String time ;5022 Vector<Float> tsysval ;5023 if ( s->nrow() == 0 ) {5024 os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ;5025 return tsys ;5026 }5027 else if ( s->nrow() == 1 ) {5028 //os << "use row " << 0 << LogIO::POST ;5029 tsysval = tsysCol( 0 ) ;5030 tsysval.tovector( tsys ) ;5031 return tsys ;5032 }5033 else {5034 vector<int> idx = getRowIdFromTime( reftime, s ) ;5035 if ( mode == "before" ) {5036 int id = -1 ;5037 if ( idx[0] != -1 ) {5038 id = idx[0] ;5039 }5040 else if ( idx[1] != -1 ) {5041 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;5042 id = idx[1] ;5043 }5044 //os << "use row " << id << LogIO::POST ;5045 tsysval = tsysCol( id ) ;5046 tsysval.tovector( tsys ) ;5047 }5048 else if ( mode == "after" ) {5049 int id = -1 ;5050 if ( idx[1] != -1 ) {5051 id = idx[1] ;5052 }5053 else if ( idx[0] != -1 ) {5054 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;5055 id = idx[1] ;5056 }5057 //os << "use row " << id << LogIO::POST ;5058 tsysval = tsysCol( id ) ;5059 tsysval.tovector( tsys ) ;5060 }5061 else if ( mode == "nearest" ) {5062 int id = -1 ;5063 if ( idx[0] == -1 ) {5064 id = idx[1] ;5065 }5066 else if ( idx[1] == -1 ) {5067 id = idx[0] ;5068 }5069 else if ( idx[0] == idx[1] ) {5070 id = idx[0] ;5071 }5072 else {5073 //double t0 = getMJD( s->getTime( idx[0] ) ) ;5074 //double t1 = getMJD( s->getTime( idx[1] ) ) ;5075 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;5076 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;5077 double tref = getMJD( reftime ) ;5078 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {5079 id = idx[1] ;5080 }5081 else {5082 id = idx[0] ;5083 }5084 }5085 //os << "use row " << id << LogIO::POST ;5086 tsysval = tsysCol( id ) ;5087 tsysval.tovector( tsys ) ;5088 }5089 else if ( mode == "linear" ) {5090 if ( idx[0] == -1 ) {5091 // use after5092 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;5093 int id = idx[1] ;5094 //os << "use row " << id << LogIO::POST ;5095 tsysval = tsysCol( id ) ;5096 tsysval.tovector( tsys ) ;5097 }5098 else if ( idx[1] == -1 ) {5099 // use before5100 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;5101 int id = idx[0] ;5102 //os << "use row " << id << LogIO::POST ;5103 tsysval = tsysCol( id ) ;5104 tsysval.tovector( tsys ) ;5105 }5106 else if ( idx[0] == idx[1] ) {5107 // use before5108 //os << "No need to interporate." << LogIO::POST ;5109 int id = idx[0] ;5110 //os << "use row " << id << LogIO::POST ;5111 tsysval = tsysCol( id ) ;5112 tsysval.tovector( tsys ) ;5113 }5114 else {5115 // do interpolation5116 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;5117 //double t0 = getMJD( s->getTime( idx[0] ) ) ;5118 //double t1 = getMJD( s->getTime( idx[1] ) ) ;5119 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;5120 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;5121 double tref = getMJD( reftime ) ;5122 vector<float> tsys0 ;5123 vector<float> tsys1 ;5124 tsysval = tsysCol( idx[0] ) ;5125 tsysval.tovector( tsys0 ) ;5126 tsysval = tsysCol( idx[1] ) ;5127 tsysval.tovector( tsys1 ) ;5128 for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {5129 float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ;5130 tsys.push_back( v ) ;5131 }5132 }5133 }5134 else {5135 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;5136 }5137 return tsys ;5138 }5139 }5140 5141 4640 Vector<Float> STMath::getTsysFromTime( double reftime, 5142 4641 const Vector<Double> &timeVec, … … 5272 4771 Vector<Double> timeHot = timeCol.getColumn() ; 5273 4772 timeCol.attach( on->table(), "TIME" ) ; 5274 ROArrayColumn<Float> tsysCol( off->table(), "SPECTRA" ) ; 5275 Matrix<Float> offspectra = tsysCol.getColumn() ; 5276 tsysCol.attach( sky->table(), "SPECTRA" ) ; 5277 Matrix<Float> skyspectra = tsysCol.getColumn() ; 5278 tsysCol.attach( hot->table(), "SPECTRA" ) ; 5279 Matrix<Float> hotspectra = tsysCol.getColumn() ; 5280 tsysCol.attach( on->table(), "TSYS" ) ; 5281 //ROArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ; 4773 ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ; 4774 Matrix<Float> offspectra = arrayFloatCol.getColumn() ; 4775 arrayFloatCol.attach( sky->table(), "SPECTRA" ) ; 4776 Matrix<Float> skyspectra = arrayFloatCol.getColumn() ; 4777 arrayFloatCol.attach( hot->table(), "SPECTRA" ) ; 4778 Matrix<Float> hotspectra = arrayFloatCol.getColumn() ; 5282 4779 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ; 5283 4780 // I know that the data is contiguous … … 5314 4811 } 5315 4812 5316 vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,5317 CountedPtr<Scantable>& off,5318 CountedPtr<Scantable>& sky,5319 CountedPtr<Scantable>& hot,5320 CountedPtr<Scantable>& cold,5321 int index,5322 string antname )5323 {5324 (void) cold; //currently unused5325 string reftime = on->getTime( index ) ;5326 vector<int> ii( 1, on->getIF( index ) ) ;5327 vector<int> ib( 1, on->getBeam( index ) ) ;5328 vector<int> ip( 1, on->getPol( index ) ) ;5329 STSelector sel = STSelector() ;5330 sel.setIFs( ii ) ;5331 sel.setBeams( ib ) ;5332 sel.setPolarizations( ip ) ;5333 sky->setSelection( sel ) ;5334 hot->setSelection( sel ) ;5335 //cold->setSelection( sel ) ;5336 off->setSelection( sel ) ;5337 vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;5338 vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;5339 //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;5340 vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;5341 vector<float> spec = on->getSpectrum( index ) ;5342 vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;5343 vector<float> sp( tcal.size() ) ;5344 if ( antname.find( "APEX" ) != string::npos ) {5345 // using gain array5346 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {5347 float v = ( ( spec[j] - spoff[j] ) / spoff[j] )5348 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;5349 sp[j] = v ;5350 }5351 }5352 else {5353 // Chopper-Wheel calibration (Ulich & Haas 1976)5354 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {5355 float v = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;5356 sp[j] = v ;5357 }5358 }5359 sel.reset() ;5360 sky->unsetSelection() ;5361 hot->unsetSelection() ;5362 //cold->unsetSelection() ;5363 off->unsetSelection() ;5364 5365 return sp ;5366 }5367 5368 4813 void STMath::calibrateALMA( CountedPtr<Scantable>& out, 5369 4814 const CountedPtr<Scantable>& on, … … 5380 4825 Vector<Double> timeVec = timeCol.getColumn() ; 5381 4826 timeCol.attach( on->table(), "TIME" ) ; 5382 ROArrayColumn<Float> tsysCol( off->table(), "SPECTRA" ) ; 5383 Matrix<Float> offspectra = tsysCol.getColumn() ; 5384 tsysCol.attach( on->table(), "TSYS" ) ; 5385 //ROArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ; 4827 ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ; 4828 Matrix<Float> offspectra = arrayFloatCol.getColumn() ; 5386 4829 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ; 5387 4830 // I know that the data is contiguous … … 5394 4837 //Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ; 5395 4838 Vector<Float> spec = on->specCol_( *p ) ; 5396 Vector<Float> tsys = tsysCol( *p ) ;4839 Vector<Float> tsys = on->tsysCol_( *p ) ; 5397 4840 // ALMA Calibration 5398 4841 // … … 5413 4856 } 5414 4857 5415 vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig, 5416 CountedPtr<Scantable>& ref, 5417 CountedPtr<Scantable>& sky, 5418 CountedPtr<Scantable>& hot, 5419 CountedPtr<Scantable>& cold, 5420 int index ) 5421 { 5422 (void) cold; //currently unused 5423 string reftime = sig->getTime( index ) ; 5424 vector<int> ii( 1, sig->getIF( index ) ) ; 5425 vector<int> ib( 1, sig->getBeam( index ) ) ; 5426 vector<int> ip( 1, sig->getPol( index ) ) ; 5427 vector<int> ic( 1, sig->getScan( index ) ) ; 5428 STSelector sel = STSelector() ; 5429 sel.setIFs( ii ) ; 5430 sel.setBeams( ib ) ; 5431 sel.setPolarizations( ip ) ; 5432 sky->setSelection( sel ) ; 5433 hot->setSelection( sel ) ; 5434 //cold->setSelection( sel ) ; 5435 vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ; 5436 vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ; 5437 //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ; 5438 vector<float> spref = ref->getSpectrum( index ) ; 5439 vector<float> spsig = sig->getSpectrum( index ) ; 5440 vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ; 5441 vector<float> sp( tcal.size() ) ; 5442 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) { 5443 float v = tcal[j] * spsky[j] / ( sphot[j] - spsky[j] ) * ( spsig[j] - spref[j] ) / spref[j] ; 5444 sp[j] = v ; 5445 } 5446 sel.reset() ; 5447 sky->unsetSelection() ; 5448 hot->unsetSelection() ; 5449 //cold->unsetSelection() ; 5450 5451 return sp ; 5452 } 5453 5454 vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig, 5455 CountedPtr<Scantable>& ref, 5456 vector< CountedPtr<Scantable> >& sky, 5457 vector< CountedPtr<Scantable> >& hot, 5458 vector< CountedPtr<Scantable> >& cold, 5459 int index ) 5460 { 5461 (void) cold; //currently unused 5462 string reftime = sig->getTime( index ) ; 5463 vector<int> ii( 1, sig->getIF( index ) ) ; 5464 vector<int> ib( 1, sig->getBeam( index ) ) ; 5465 vector<int> ip( 1, sig->getPol( index ) ) ; 5466 vector<int> ic( 1, sig->getScan( index ) ) ; 5467 STSelector sel = STSelector() ; 5468 sel.setIFs( ii ) ; 5469 sel.setBeams( ib ) ; 5470 sel.setPolarizations( ip ) ; 5471 sky[0]->setSelection( sel ) ; 5472 hot[0]->setSelection( sel ) ; 5473 //cold[0]->setSelection( sel ) ; 5474 vector<float> spskys = getSpectrumFromTime( reftime, sky[0], "linear" ) ; 5475 vector<float> sphots = getSpectrumFromTime( reftime, hot[0], "linear" ) ; 5476 //vector<float> spcolds = getSpectrumFromTime( reftime, cold[0], "linear" ) ; 5477 vector<float> tcals = getTcalFromTime( reftime, sky[0], "linear" ) ; 5478 sel.reset() ; 5479 ii[0] = ref->getIF( index ) ; 5480 sel.setIFs( ii ) ; 5481 sel.setBeams( ib ) ; 5482 sel.setPolarizations( ip ) ; 5483 sky[1]->setSelection( sel ) ; 5484 hot[1]->setSelection( sel ) ; 5485 //cold[1]->setSelection( sel ) ; 5486 vector<float> spskyr = getSpectrumFromTime( reftime, sky[1], "linear" ) ; 5487 vector<float> sphotr = getSpectrumFromTime( reftime, hot[1], "linear" ) ; 5488 //vector<float> spcoldr = getSpectrumFromTime( reftime, cold[1], "linear" ) ; 5489 vector<float> tcalr = getTcalFromTime( reftime, sky[1], "linear" ) ; 5490 vector<float> spref = ref->getSpectrum( index ) ; 5491 vector<float> spsig = sig->getSpectrum( index ) ; 5492 vector<float> sp( tcals.size() ) ; 5493 for ( unsigned int j = 0 ; j < tcals.size() ; j++ ) { 5494 float v = tcals[j] * spsig[j] / ( sphots[j] - spskys[j] ) - tcalr[j] * spref[j] / ( sphotr[j] - spskyr[j] ) ; 5495 sp[j] = v ; 5496 } 5497 sel.reset() ; 5498 sky[0]->unsetSelection() ; 5499 hot[0]->unsetSelection() ; 5500 //cold[0]->unsetSelection() ; 5501 sky[1]->unsetSelection() ; 5502 hot[1]->unsetSelection() ; 5503 //cold[1]->unsetSelection() ; 5504 5505 return sp ; 4858 void STMath::calibrateAPEXFS( CountedPtr<Scantable> &sig, 4859 CountedPtr<Scantable> &ref, 4860 const vector< CountedPtr<Scantable> >& on, 4861 const vector< CountedPtr<Scantable> >& sky, 4862 const vector< CountedPtr<Scantable> >& hot, 4863 const vector< CountedPtr<Scantable> >& cold, 4864 const Vector<uInt> &rows ) 4865 { 4866 // if rows is empty, just return 4867 if ( rows.nelements() == 0 ) 4868 return ; 4869 ROScalarColumn<Double> timeCol( sky[0]->table(), "TIME" ) ; 4870 Vector<Double> timeSkyS = timeCol.getColumn() ; 4871 timeCol.attach( sky[1]->table(), "TIME" ) ; 4872 Vector<Double> timeSkyR = timeCol.getColumn() ; 4873 timeCol.attach( hot[0]->table(), "TIME" ) ; 4874 Vector<Double> timeHotS = timeCol.getColumn() ; 4875 timeCol.attach( hot[1]->table(), "TIME" ) ; 4876 Vector<Double> timeHotR = timeCol.getColumn() ; 4877 timeCol.attach( sig->table(), "TIME" ) ; 4878 ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ; 4879 ROArrayColumn<Float> arrayFloatCol( sky[0]->table(), "SPECTRA" ) ; 4880 Matrix<Float> skyspectraS = arrayFloatCol.getColumn() ; 4881 arrayFloatCol.attach( sky[1]->table(), "SPECTRA" ) ; 4882 Matrix<Float> skyspectraR = arrayFloatCol.getColumn() ; 4883 arrayFloatCol.attach( hot[0]->table(), "SPECTRA" ) ; 4884 Matrix<Float> hotspectraS = arrayFloatCol.getColumn() ; 4885 arrayFloatCol.attach( hot[1]->table(), "SPECTRA" ) ; 4886 Matrix<Float> hotspectraR = arrayFloatCol.getColumn() ; 4887 unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ; 4888 Vector<Float> spec( spsize ) ; 4889 // I know that the data is contiguous 4890 const uInt *p = rows.data() ; 4891 vector<int> ids( 2 ) ; 4892 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) { 4893 double reftime = timeCol.asdouble(*p) ; 4894 ids = getRowIdFromTime( reftime, timeSkyS ) ; 4895 Vector<Float> spskyS = getSpectrumFromTime( reftime, timeSkyS, ids, skyspectraS, "linear" ) ; 4896 Vector<Float> tcalS = getTcalFromTime( reftime, timeSkyS, ids, sky[0], "linear" ) ; 4897 Vector<Float> tsysS = getTsysFromTime( reftime, timeSkyS, ids, sky[0], "linear" ) ; 4898 ids = getRowIdFromTime( reftime, timeHotS ) ; 4899 Vector<Float> sphotS = getSpectrumFromTime( reftime, timeHotS, ids, hotspectraS ) ; 4900 reftime = timeCol2.asdouble(*p) ; 4901 ids = getRowIdFromTime( reftime, timeSkyR ) ; 4902 Vector<Float> spskyR = getSpectrumFromTime( reftime, timeSkyR, ids, skyspectraR, "linear" ) ; 4903 Vector<Float> tcalR = getTcalFromTime( reftime, timeSkyR, ids, sky[1], "linear" ) ; 4904 Vector<Float> tsysR = getTsysFromTime( reftime, timeSkyR, ids, sky[1], "linear" ) ; 4905 ids = getRowIdFromTime( reftime, timeHotR ) ; 4906 Vector<Float> sphotR = getSpectrumFromTime( reftime, timeHotR, ids, hotspectraR ) ; 4907 Vector<Float> spsig = on[0]->specCol_( *p ) ; 4908 Vector<Float> spref = on[1]->specCol_( *p ) ; 4909 for ( unsigned int j = 0 ; j < spsize ; j++ ) { 4910 spec[j] = tcalS[j] * spsig[j] / ( sphotS[j] - spskyS[j] ) 4911 - tcalR[j] * spref[j] / ( sphotR[j] - spskyR[j] ) ; 4912 } 4913 sig->specCol_.put( *p, spec ) ; 4914 sig->tsysCol_.put( *p, tsysS ) ; 4915 spec *= (Float)-1.0 ; 4916 ref->specCol_.put( *p, spec ) ; 4917 ref->tsysCol_.put( *p, tsysR ) ; 4918 p++ ; 4919 } 4920 } 4921 4922 void STMath::calibrateFS( CountedPtr<Scantable> &sig, 4923 CountedPtr<Scantable> &ref, 4924 const CountedPtr<Scantable>& rsig, 4925 const CountedPtr<Scantable>& rref, 4926 const CountedPtr<Scantable>& sky, 4927 const CountedPtr<Scantable>& hot, 4928 const CountedPtr<Scantable>& cold, 4929 const Vector<uInt> &rows ) 4930 { 4931 // if rows is empty, just return 4932 if ( rows.nelements() == 0 ) 4933 return ; 4934 ROScalarColumn<Double> timeCol( sky->table(), "TIME" ) ; 4935 Vector<Double> timeSky = timeCol.getColumn() ; 4936 timeCol.attach( hot->table(), "TIME" ) ; 4937 Vector<Double> timeHot = timeCol.getColumn() ; 4938 timeCol.attach( sig->table(), "TIME" ) ; 4939 ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ; 4940 ROArrayColumn<Float> arrayFloatCol( sky->table(), "SPECTRA" ) ; 4941 Matrix<Float> skyspectra = arrayFloatCol.getColumn() ; 4942 arrayFloatCol.attach( hot->table(), "SPECTRA" ) ; 4943 Matrix<Float> hotspectra = arrayFloatCol.getColumn() ; 4944 unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ; 4945 Vector<Float> spec( spsize ) ; 4946 // I know that the data is contiguous 4947 const uInt *p = rows.data() ; 4948 vector<int> ids( 2 ) ; 4949 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) { 4950 double reftime = timeCol.asdouble(*p) ; 4951 ids = getRowIdFromTime( reftime, timeSky ) ; 4952 Vector<Float> spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ; 4953 Vector<Float> tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ; 4954 Vector<Float> tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ; 4955 ids = getRowIdFromTime( reftime, timeHot ) ; 4956 Vector<Float> sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra ) ; 4957 Vector<Float> spsig = rsig->specCol_( *p ) ; 4958 Vector<Float> spref = rref->specCol_( *p ) ; 4959 // using gain array 4960 for ( unsigned int j = 0 ; j < spsize ; j++ ) { 4961 spec[j] = ( ( spsig[j] - spref[j] ) / spref[j] ) 4962 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ; 4963 } 4964 sig->specCol_.put( *p, spec ) ; 4965 sig->tsysCol_.put( *p, tsys ) ; 4966 4967 reftime = timeCol2.asdouble(*p) ; 4968 spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ; 4969 tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ; 4970 tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ; 4971 ids = getRowIdFromTime( reftime, timeHot ) ; 4972 sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra ) ; 4973 // using gain array 4974 for ( unsigned int j = 0 ; j < spsize ; j++ ) { 4975 spec[j] = ( ( spref[j] - spsig[j] ) / spsig[j] ) 4976 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ; 4977 } 4978 ref->specCol_.put( *p, spec ) ; 4979 ref->tsysCol_.put( *p, tsys ) ; 4980 p++ ; 4981 } 5506 4982 } 5507 4983 -
branches/hpc33/src/STMath.h
r2563 r2566 402 402 flagsFromMA(const casa::MaskedArray<casa::Float>& ma); 403 403 404 vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;405 404 casa::Vector<casa::Float> getSpectrumFromTime( double reftime, const casa::Vector<casa::Double> &v, const vector<int> &idx, const casa::Matrix<casa::Float>& sp, string mode = "before" ) ; 406 vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;407 405 casa::Vector<casa::Float> getTcalFromTime( double reftime, const casa::Vector<casa::Double> &v, const vector<int> &idx, const casa::CountedPtr<Scantable>& s, string mode="before" ) ; 408 vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;409 406 casa::Vector<casa::Float> getTsysFromTime( double reftime, const casa::Vector<casa::Double> &v, const vector<int> &idx, const casa::CountedPtr<Scantable>& s, string mode="before" ) ; 410 vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;411 407 vector<int> getRowIdFromTime( double reftime, const casa::Vector<casa::Double>& t ) ; 412 408 … … 420 416 const casa::Vector<casa::uInt> &rows, 421 417 const casa::String &antname ) ; 422 vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on, 423 casa::CountedPtr<Scantable>& off, 424 casa::CountedPtr<Scantable>& sky, 425 casa::CountedPtr<Scantable>& hot, 426 casa::CountedPtr<Scantable>& cold, 427 int index, 428 string antname ) ; 418 429 419 // Tsys * (ON-OFF)/OFF 430 420 void calibrateALMA( casa::CountedPtr<Scantable>& out, … … 432 422 const casa::CountedPtr<Scantable>& off, 433 423 const casa::Vector<casa::uInt>& rows ) ; 434 vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig, 435 casa::CountedPtr<Scantable>& ref, 436 casa::CountedPtr<Scantable>& sky, 437 casa::CountedPtr<Scantable>& hot, 438 casa::CountedPtr<Scantable>& cold, 439 int index ) ; 440 vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig, 441 casa::CountedPtr<Scantable>& ref, 442 vector< casa::CountedPtr<Scantable> >& sky, 443 vector< casa::CountedPtr<Scantable> >& hot, 444 vector< casa::CountedPtr<Scantable> >& cold, 445 int index ) ; 446 double getMJD( string strtime ) ; 424 425 // Frequency switching 426 void calibrateFS( casa::CountedPtr<Scantable> &sig, 427 casa::CountedPtr<Scantable> &ref, 428 const casa::CountedPtr<Scantable> &rsig, 429 const casa::CountedPtr<Scantable> &rref, 430 const casa::CountedPtr<Scantable> &sky, 431 const casa::CountedPtr<Scantable> &hot, 432 const casa::CountedPtr<Scantable> &cold, 433 const casa::Vector<casa::uInt> &rows ) ; 434 void calibrateAPEXFS( casa::CountedPtr<Scantable> &sig, 435 casa::CountedPtr<Scantable> &ref, 436 const vector< casa::CountedPtr<Scantable> > &on, 437 const vector< casa::CountedPtr<Scantable> > &sky, 438 const vector< casa::CountedPtr<Scantable> > &hot, 439 const vector< casa::CountedPtr<Scantable> > &cold, 440 const casa::Vector<casa::uInt> &rows ) ; 447 441 void copyRows( casa::Table &out, 448 442 const casa::Table &in,
Note:
See TracChangeset
for help on using the changeset viewer.