Changeset 2432 for branches/hpc33
- Timestamp:
- 03/14/12 19:20:33 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/hpc33/src/STMath.cpp
r2345 r2432 57 57 using namespace asap; 58 58 59 // 2012/02/17 TN 60 // Since STGrid is implemented, average doesn't consider direction 61 // when accumulating 59 62 // tolerance for direction comparison (rad) 60 #define TOL_OTF 1.0e-1561 #define TOL_POINT 2.9088821e-4 // 1 arcmin63 // #define TOL_OTF 1.0e-15 64 // #define TOL_POINT 2.9088821e-4 // 1 arcmin 62 65 63 66 STMath::STMath(bool insitu) : … … 83 86 WeightType wtype = stringToWeight(weight); 84 87 88 // 2012/02/17 TN 89 // Since STGrid is implemented, average doesn't consider direction 90 // when accumulating 85 91 // check if OTF observation 86 String obstype = in[0]->getHeader().obstype ;87 Double tol = 0.0 ;88 if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {89 tol = TOL_OTF ;90 }91 else {92 tol = TOL_POINT ;93 }92 // String obstype = in[0]->getHeader().obstype ; 93 // Double tol = 0.0 ; 94 // if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) { 95 // tol = TOL_OTF ; 96 // } 97 // else { 98 // tol = TOL_POINT ; 99 // } 94 100 95 101 // output … … 142 148 while (!iter.pastEnd()) { 143 149 Table subt = iter.table(); 144 // // copy the first row of this selection into the new table 145 // tout.addRow(); 146 // TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 147 // // re-index to 0 148 // if ( avmode != "SCAN" && avmode != "SOURCE" ) { 149 // scanColOut.put(outrowCount, uInt(0)); 150 // copy the first row of this selection into the new table 151 tout.addRow(); 152 TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 153 // re-index to 0 154 if ( avmode != "SCAN" && avmode != "SOURCE" ) { 155 scanColOut.put(outrowCount, uInt(0)); 156 } 157 ++outrowCount; 158 // 2012/02/17 TN 159 // Since STGrid is implemented, average doesn't consider direction 160 // when accumulating 161 // MDirection::ScalarColumn dircol ; 162 // dircol.attach( subt, "DIRECTION" ) ; 163 // Int length = subt.nrow() ; 164 // vector< Vector<Double> > dirs ; 165 // vector<int> indexes ; 166 // for ( Int i = 0 ; i < length ; i++ ) { 167 // Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ; 168 // //os << << count++ << ": " ; 169 // //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ; 170 // bool adddir = true ; 171 // for ( uInt j = 0 ; j < dirs.size() ; j++ ) { 172 // //if ( allTrue( t == dirs[j] ) ) { 173 // Double dx = t[0] - dirs[j][0] ; 174 // Double dy = t[1] - dirs[j][1] ; 175 // Double dd = sqrt( dx * dx + dy * dy ) ; 176 // //if ( allNearAbs( t, dirs[j], tol ) ) { 177 // if ( dd <= tol ) { 178 // adddir = false ; 179 // break ; 180 // } 181 // } 182 // if ( adddir ) { 183 // dirs.push_back( t ) ; 184 // indexes.push_back( i ) ; 185 // } 150 186 // } 151 // ++outrowCount; 152 MDirection::ScalarColumn dircol ; 153 dircol.attach( subt, "DIRECTION" ) ; 154 Int length = subt.nrow() ; 155 vector< Vector<Double> > dirs ; 156 vector<int> indexes ; 157 for ( Int i = 0 ; i < length ; i++ ) { 158 Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ; 159 //os << << count++ << ": " ; 160 //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ; 161 bool adddir = true ; 162 for ( uInt j = 0 ; j < dirs.size() ; j++ ) { 163 //if ( allTrue( t == dirs[j] ) ) { 164 Double dx = t[0] - dirs[j][0] ; 165 Double dy = t[1] - dirs[j][1] ; 166 Double dd = sqrt( dx * dx + dy * dy ) ; 167 //if ( allNearAbs( t, dirs[j], tol ) ) { 168 if ( dd <= tol ) { 169 adddir = false ; 170 break ; 171 } 172 } 173 if ( adddir ) { 174 dirs.push_back( t ) ; 175 indexes.push_back( i ) ; 176 } 177 } 178 uInt rowNum = dirs.size() ; 179 tout.addRow( rowNum ) ; 180 for ( uInt i = 0 ; i < rowNum ; i++ ) { 181 TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ; 182 // re-index to 0 183 if ( avmode != "SCAN" && avmode != "SOURCE" ) { 184 scanColOut.put(outrowCount+i, uInt(0)); 185 } 186 } 187 outrowCount += rowNum ; 187 // uInt rowNum = dirs.size() ; 188 // tout.addRow( rowNum ) ; 189 // for ( uInt i = 0 ; i < rowNum ; i++ ) { 190 // TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ; 191 // // re-index to 0 192 // if ( avmode != "SCAN" && avmode != "SOURCE" ) { 193 // scanColOut.put(outrowCount+i, uInt(0)); 194 // } 195 // } 196 // outrowCount += rowNum ; 188 197 ++iter; 189 198 } … … 218 227 } 219 228 220 vector<uInt> removeRows ; 221 uInt nrsubt = subt.nrow() ; 222 for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) { 223 //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) { 224 Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ; 225 Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ; 226 double dx = x0[0] - x1[0]; 227 double dy = x0[1] - x1[1]; 228 Double dd = sqrt( dx * dx + dy * dy ) ; 229 //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) { 230 if ( dd > tol ) { 231 removeRows.push_back( irow ) ; 232 } 233 } 234 if ( removeRows.size() != 0 ) { 235 subt.removeRow( removeRows ) ; 236 } 229 // 2012/02/17 TN 230 // Since STGrid is implemented, average doesn't consider direction 231 // when accumulating 232 // vector<uInt> removeRows ; 233 // uInt nrsubt = subt.nrow() ; 234 // for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) { 235 // //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) { 236 // Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ; 237 // Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ; 238 // double dx = x0[0] - x1[0]; 239 // double dy = x0[1] - x1[1]; 240 // Double dd = sqrt( dx * dx + dy * dy ) ; 241 // //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) { 242 // if ( dd > tol ) { 243 // removeRows.push_back( irow ) ; 244 // } 245 // } 246 // if ( removeRows.size() != 0 ) { 247 // subt.removeRow( removeRows ) ; 248 // } 237 249 238 if ( nrsubt == removeRows.size() )239 throw(AipsError("Averaging data is empty.")) ;250 // if ( nrsubt == removeRows.size() ) 251 // throw(AipsError("Averaging data is empty.")) ; 240 252 241 253 specCol.attach(subt,"SPECTRA"); … … 321 333 { 322 334 (void) mode; // currently unused 335 // 2012/02/17 TN 336 // Since STGrid is implemented, average doesn't consider direction 337 // when accumulating 323 338 // check if OTF observation 324 String obstype = in->getHeader().obstype ;325 Double tol = 0.0 ;326 if ( obstype.find( "OTF" ) != String::npos ) {327 tol = TOL_OTF ;328 }329 else {330 tol = TOL_POINT ;331 }339 // String obstype = in->getHeader().obstype ; 340 // Double tol = 0.0 ; 341 // if ( obstype.find( "OTF" ) != String::npos ) { 342 // tol = TOL_OTF ; 343 // } 344 // else { 345 // tol = TOL_POINT ; 346 // } 332 347 333 348 // clone as this is non insitu … … 362 377 flagCol.attach(subt,"FLAGTRA"); 363 378 tsysCol.attach(subt,"TSYS"); 364 // tout.addRow(); 365 // TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 366 // if ( avmode != "SCAN") { 367 // scanColOut.put(outrowCount, uInt(0)); 379 380 tout.addRow(); 381 TableCopy::copyRows(tout, subt, outrowCount, 0, 1); 382 if ( avmode != "SCAN") { 383 scanColOut.put(outrowCount, uInt(0)); 384 } 385 Vector<Float> tmp; 386 specCol.get(0, tmp); 387 uInt nchan = tmp.nelements(); 388 // have to do channel by channel here as MaskedArrMath 389 // doesn't have partialMedians 390 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0))); 391 Vector<Float> outspec(nchan); 392 Vector<uChar> outflag(nchan,0); 393 Vector<Float> outtsys(1);/// @fixme when tsys is channel based 394 for (uInt i=0; i<nchan; ++i) { 395 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i))); 396 MaskedArray<Float> ma = maskedArray(specs,flags); 397 outspec[i] = median(ma); 398 if ( allEQ(ma.getMask(), False) ) 399 outflag[i] = userflag;// flag data 400 } 401 outtsys[0] = median(tsysCol.getColumn()); 402 specColOut.put(outrowCount, outspec); 403 flagColOut.put(outrowCount, outflag); 404 tsysColOut.put(outrowCount, outtsys); 405 Double intsum = sum(intCol.getColumn()); 406 intColOut.put(outrowCount, intsum); 407 ++outrowCount; 408 ++iter; 409 410 // 2012/02/17 TN 411 // Since STGrid is implemented, average doesn't consider direction 412 // when accumulating 413 // MDirection::ScalarColumn dircol ; 414 // dircol.attach( subt, "DIRECTION" ) ; 415 // Int length = subt.nrow() ; 416 // vector< Vector<Double> > dirs ; 417 // vector<int> indexes ; 418 // // Handle MX mode averaging 419 // if (in->nbeam() > 1 ) { 420 // length = 1; 368 421 // } 369 // Vector<Float> tmp; 370 // specCol.get(0, tmp); 371 // uInt nchan = tmp.nelements(); 372 // // have to do channel by channel here as MaskedArrMath 373 // // doesn't have partialMedians 374 // Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0))); 375 // Vector<Float> outspec(nchan); 376 // Vector<uChar> outflag(nchan,0); 377 // Vector<Float> outtsys(1);/// @fixme when tsys is channel based 378 // for (uInt i=0; i<nchan; ++i) { 379 // Vector<Float> specs = specCol.getColumn(Slicer(Slice(i))); 380 // MaskedArray<Float> ma = maskedArray(specs,flags); 381 // outspec[i] = median(ma); 382 // if ( allEQ(ma.getMask(), False) ) 383 // outflag[i] = userflag;// flag data 422 // for ( Int i = 0 ; i < length ; i++ ) { 423 // Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ; 424 // bool adddir = true ; 425 // for ( uInt j = 0 ; j < dirs.size() ; j++ ) { 426 // //if ( allTrue( t == dirs[j] ) ) { 427 // Double dx = t[0] - dirs[j][0] ; 428 // Double dy = t[1] - dirs[j][1] ; 429 // Double dd = sqrt( dx * dx + dy * dy ) ; 430 // //if ( allNearAbs( t, dirs[j], tol ) ) { 431 // if ( dd <= tol ) { 432 // adddir = false ; 433 // break ; 434 // } 435 // } 436 // if ( adddir ) { 437 // dirs.push_back( t ) ; 438 // indexes.push_back( i ) ; 439 // } 384 440 // } 385 // outtsys[0] = median(tsysCol.getColumn()); 386 // specColOut.put(outrowCount, outspec); 387 // flagColOut.put(outrowCount, outflag); 388 // tsysColOut.put(outrowCount, outtsys); 389 // Double intsum = sum(intCol.getColumn()); 390 // intColOut.put(outrowCount, intsum); 391 // ++outrowCount; 441 // uInt rowNum = dirs.size() ; 442 // tout.addRow( rowNum ); 443 // for ( uInt i = 0 ; i < rowNum ; i++ ) { 444 // TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ; 445 // // Handle MX mode averaging 446 // if ( avmode != "SCAN") { 447 // scanColOut.put(outrowCount+i, uInt(0)); 448 // } 449 // } 450 // MDirection::ScalarColumn dircolOut ; 451 // dircolOut.attach( tout, "DIRECTION" ) ; 452 // for ( uInt irow = 0 ; irow < rowNum ; irow++ ) { 453 // Vector<Double> t = \ 454 // dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ; 455 // Vector<Float> tmp; 456 // specCol.get(0, tmp); 457 // uInt nchan = tmp.nelements(); 458 // // have to do channel by channel here as MaskedArrMath 459 // // doesn't have partialMedians 460 // Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0))); 461 // // mask spectra for different DIRECTION 462 // for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) { 463 // Vector<Double> direction = \ 464 // dircol(jrow).getAngle(Unit(String("rad"))).getValue() ; 465 // //if ( t[0] != direction[0] || t[1] != direction[1] ) { 466 // Double dx = t[0] - direction[0]; 467 // Double dy = t[1] - direction[1]; 468 // Double dd = sqrt(dx*dx + dy*dy); 469 // //if ( !allNearAbs( t, direction, tol ) ) { 470 // if ( dd > tol && in->nbeam() < 2 ) { 471 // flags[jrow] = userflag ; 472 // } 473 // } 474 // Vector<Float> outspec(nchan); 475 // Vector<uChar> outflag(nchan,0); 476 // Vector<Float> outtsys(1);/// @fixme when tsys is channel based 477 // for (uInt i=0; i<nchan; ++i) { 478 // Vector<Float> specs = specCol.getColumn(Slicer(Slice(i))); 479 // MaskedArray<Float> ma = maskedArray(specs,flags); 480 // outspec[i] = median(ma); 481 // if ( allEQ(ma.getMask(), False) ) 482 // outflag[i] = userflag;// flag data 483 // } 484 // outtsys[0] = median(tsysCol.getColumn()); 485 // specColOut.put(outrowCount+irow, outspec); 486 // flagColOut.put(outrowCount+irow, outflag); 487 // tsysColOut.put(outrowCount+irow, outtsys); 488 // Vector<Double> integ = intCol.getColumn() ; 489 // MaskedArray<Double> mi = maskedArray( integ, flags ) ; 490 // Double intsum = sum(mi); 491 // intColOut.put(outrowCount+irow, intsum); 492 // } 493 // outrowCount += rowNum ; 392 494 // ++iter; 393 MDirection::ScalarColumn dircol ;394 dircol.attach( subt, "DIRECTION" ) ;395 Int length = subt.nrow() ;396 vector< Vector<Double> > dirs ;397 vector<int> indexes ;398 // Handle MX mode averaging399 if (in->nbeam() > 1 ) {400 length = 1;401 }402 for ( Int i = 0 ; i < length ; i++ ) {403 Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;404 bool adddir = true ;405 for ( uInt j = 0 ; j < dirs.size() ; j++ ) {406 //if ( allTrue( t == dirs[j] ) ) {407 Double dx = t[0] - dirs[j][0] ;408 Double dy = t[1] - dirs[j][1] ;409 Double dd = sqrt( dx * dx + dy * dy ) ;410 //if ( allNearAbs( t, dirs[j], tol ) ) {411 if ( dd <= tol ) {412 adddir = false ;413 break ;414 }415 }416 if ( adddir ) {417 dirs.push_back( t ) ;418 indexes.push_back( i ) ;419 }420 }421 uInt rowNum = dirs.size() ;422 tout.addRow( rowNum );423 for ( uInt i = 0 ; i < rowNum ; i++ ) {424 TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;425 // Handle MX mode averaging426 if ( avmode != "SCAN") {427 scanColOut.put(outrowCount+i, uInt(0));428 }429 }430 MDirection::ScalarColumn dircolOut ;431 dircolOut.attach( tout, "DIRECTION" ) ;432 for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {433 Vector<Double> t = \434 dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;435 Vector<Float> tmp;436 specCol.get(0, tmp);437 uInt nchan = tmp.nelements();438 // have to do channel by channel here as MaskedArrMath439 // doesn't have partialMedians440 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));441 // mask spectra for different DIRECTION442 for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {443 Vector<Double> direction = \444 dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;445 //if ( t[0] != direction[0] || t[1] != direction[1] ) {446 Double dx = t[0] - direction[0];447 Double dy = t[1] - direction[1];448 Double dd = sqrt(dx*dx + dy*dy);449 //if ( !allNearAbs( t, direction, tol ) ) {450 if ( dd > tol && in->nbeam() < 2 ) {451 flags[jrow] = userflag ;452 }453 }454 Vector<Float> outspec(nchan);455 Vector<uChar> outflag(nchan,0);456 Vector<Float> outtsys(1);/// @fixme when tsys is channel based457 for (uInt i=0; i<nchan; ++i) {458 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));459 MaskedArray<Float> ma = maskedArray(specs,flags);460 outspec[i] = median(ma);461 if ( allEQ(ma.getMask(), False) )462 outflag[i] = userflag;// flag data463 }464 outtsys[0] = median(tsysCol.getColumn());465 specColOut.put(outrowCount+irow, outspec);466 flagColOut.put(outrowCount+irow, outflag);467 tsysColOut.put(outrowCount+irow, outtsys);468 Vector<Double> integ = intCol.getColumn() ;469 MaskedArray<Double> mi = maskedArray( integ, flags ) ;470 Double intsum = sum(mi);471 intColOut.put(outrowCount+irow, intsum);472 }473 outrowCount += rowNum ;474 ++iter;475 495 } 476 496 return out; … … 2954 2974 "Use merge first.")); 2955 2975 2976 // 2012/02/17 TN 2977 // Since STGrid is implemented, average doesn't consider direction 2978 // when accumulating 2956 2979 // check if OTF observation 2957 String obstype = in[0]->getHeader().obstype ;2958 Double tol = 0.0 ;2959 if ( obstype.find( "OTF" ) != String::npos ) {2960 tol = TOL_OTF ;2961 }2962 else {2963 tol = TOL_POINT ;2964 }2980 // String obstype = in[0]->getHeader().obstype ; 2981 // Double tol = 0.0 ; 2982 // if ( obstype.find( "OTF" ) != String::npos ) { 2983 // tol = TOL_OTF ; 2984 // } 2985 // else { 2986 // tol = TOL_POINT ; 2987 // } 2965 2988 2966 2989 CountedPtr<Scantable> out ; // processed result … … 3557 3580 ScalarColumn<uInt> freqidColOut ; 3558 3581 freqidColOut.attach( out->table(), "FREQ_ID" ) ; 3559 MDirection::ScalarColumn dirColOut ;3560 dirColOut.attach( out->table(), "DIRECTION" ) ;3582 // MDirection::ScalarColumn dirColOut ; 3583 // dirColOut.attach( out->table(), "DIRECTION" ) ; 3561 3584 Table &tab = tmpout->table() ; 3562 3585 Block<String> cols(1); … … 3575 3598 ScalarColumn<uInt> polnos ; 3576 3599 polnos.attach( iter.table(), "POLNO" ) ; 3577 MDirection::ScalarColumn dircol ;3578 dircol.attach( iter.table(), "DIRECTION" ) ;3600 // MDirection::ScalarColumn dircol ; 3601 // dircol.attach( iter.table(), "DIRECTION" ) ; 3579 3602 uInt polno = polnos( 0 ) ; 3580 3603 //os << "POLNO iteration: " << polno << LogIO::POST ; … … 3635 3658 if ( polout == polno ) { 3636 3659 uInt ifout = ifnoColOut( irow ) ; 3637 Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;3660 // Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ; 3638 3661 uInt igrp ; 3639 3662 for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) { … … 3646 3669 for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) { 3647 3670 uInt ifno = ifnoCol( jrow ) ; 3648 Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ; 3649 //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction ) ) { 3650 Double dx = tdir[0] - direction[0] ; 3651 Double dy = tdir[1] - direction[1] ; 3652 Double dd = sqrt( dx * dx + dy * dy ) ; 3671 // 2012/02/17 TN 3672 // Since STGrid is implemented, average doesn't consider direction 3673 // when accumulating 3674 // Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ; 3675 // //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction ) ) { 3676 // Double dx = tdir[0] - direction[0] ; 3677 // Double dy = tdir[1] - direction[1] ; 3678 // Double dd = sqrt( dx * dx + dy * dy ) ; 3653 3679 //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) { 3654 if ( ifno == freqgrp[igrp][imem] && dd <= tol ) { 3680 // if ( ifno == freqgrp[igrp][imem] && dd <= tol ) { 3681 if ( ifno == freqgrp[igrp][imem] ) { 3655 3682 Vector<Float> spec = specCols( jrow ) ; 3656 3683 Vector<uChar> flag = flagCols( jrow ) ;
Note:
See TracChangeset
for help on using the changeset viewer.