- Timestamp:
- 06/12/12 18:45:52 (13 years ago)
- Location:
- branches/hpc33/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/hpc33/src/STMath.cpp
r2561 r2563 3843 3843 vector<int> types ; 3844 3844 3845 // save original table selection 3846 Table torg = s->table_ ; 3847 3845 3848 // sky scan 3846 STSelector sel = STSelector();3847 types.push_back( SrcType::SKY );3848 sel.setTypes( types ) ;3849 s->setSelection( sel) ;3850 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) );3851 CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN") ;3852 s->unsetSelection() ;3853 sel.reset() ;3854 types.clear() ;3855 3849 bool insitu = insitu_ ; 3850 insitu_ = false ; 3851 // share calibration scans before average with out 3852 CountedPtr<Scantable> out = getScantable( s, true ) ; 3853 insitu_ = insitu ; 3854 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ; 3855 out->attach() ; 3856 CountedPtr<Scantable> asky = averageWithinSession( out, 3857 masks, 3858 "TINT" ) ; 3856 3859 // hot scan 3857 types.push_back( SrcType::HOT ) ; 3858 sel.setTypes( types ) ; 3859 s->setSelection( sel ) ; 3860 tmp.clear() ; 3861 tmp.push_back( getScantable( s, false ) ) ; 3862 CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ; 3863 s->unsetSelection() ; 3864 sel.reset() ; 3865 types.clear() ; 3866 3860 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ; 3861 out->attach() ; 3862 CountedPtr<Scantable> ahot = averageWithinSession( out, 3863 masks, 3864 "TINT" ) ; 3867 3865 // cold scan 3868 3866 CountedPtr<Scantable> acold ; 3869 // types.push_back( SrcType::COLD ) ; 3870 // sel.setTypes( types ) ; 3871 // s->setSelection( sel ) ; 3872 // tmp.clear() ; 3873 // tmp.push_back( getScantable( s, false ) ) ; 3874 // CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCNAN" ) ; 3875 // s->unsetSelection() ; 3876 // sel.reset() ; 3877 // types.clear() ; 3867 // out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ; 3868 // out->attach() ; 3869 // CountedPtr<Scantable> acold = averageWithinSession( out, 3870 // masks, 3871 // "TINT" ) ; 3878 3872 3879 3873 // off scan 3880 types.push_back( SrcType::PSOFF ) ; 3881 sel.setTypes( types ) ; 3882 s->setSelection( sel ) ; 3883 tmp.clear() ; 3884 tmp.push_back( getScantable( s, false ) ) ; 3885 CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ; 3886 s->unsetSelection() ; 3887 sel.reset() ; 3888 types.clear() ; 3874 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ; 3875 out->attach() ; 3876 CountedPtr<Scantable> aoff = averageWithinSession( out, 3877 masks, 3878 "TINT" ) ; 3889 3879 3890 3880 // on scan 3891 bool insitu = insitu_ ; 3892 insitu_ = false ; 3893 CountedPtr<Scantable> out = getScantable( s, true ) ; 3894 insitu_ = insitu ; 3895 types.push_back( SrcType::PSON ) ; 3896 sel.setTypes( types ) ; 3897 s->setSelection( sel ) ; 3898 TableCopy::copyRows( out->table(), s->table() ) ; 3899 s->unsetSelection() ; 3900 sel.reset() ; 3901 types.clear() ; 3881 s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ; 3882 s->attach() ; 3883 out->table_ = out->originalTable_ ; 3884 out->attach() ; 3885 out->table().addRow( s->nrow() ) ; 3886 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ; 3902 3887 3903 3888 // process each on scan 3904 ArrayColumn<Float> tsysCol ; 3905 tsysCol.attach( out->table(), "TSYS" ) ; 3906 for ( int i = 0 ; i < out->nrow() ; i++ ) { 3907 vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ; 3908 out->setSpectrum( sp, i ) ; 3909 string reftime = out->getTime( i ) ; 3910 vector<int> ii( 1, out->getIF( i ) ) ; 3911 vector<int> ib( 1, out->getBeam( i ) ) ; 3912 vector<int> ip( 1, out->getPol( i ) ) ; 3913 sel.setIFs( ii ) ; 3914 sel.setBeams( ib ) ; 3915 sel.setPolarizations( ip ) ; 3916 asky->setSelection( sel ) ; 3917 vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ; 3918 const Vector<Float> Vtsys( sptsys ) ; 3919 tsysCol.put( i, Vtsys ) ; 3889 STSelector sel ; 3890 vector<string> cols( 3 ) ; 3891 cols[0] = "BEAMNO" ; 3892 cols[1] = "POLNO" ; 3893 cols[2] = "IFNO" ; 3894 STIdxIter *iter = new STIdxIterAcc( out, cols ) ; 3895 while ( !iter->pastEnd() ) { 3896 Vector<uInt> ids = iter->current() ; 3897 stringstream ss ; 3898 ss << "SELECT FROM $1 WHERE " 3899 << "BEAMNO==" << ids[0] << "&&" 3900 << "POLNO==" << ids[1] << "&&" 3901 << "IFNO==" << ids[2] ; 3902 //cout << "TaQL string: " << ss.str() << endl ; 3903 sel.setTaQL( ss.str() ) ; 3904 aoff->setSelection( sel ) ; 3905 ahot->setSelection( sel ) ; 3906 asky->setSelection( sel ) ; 3907 Vector<uInt> rows = iter->getRows( SHARE ) ; 3908 // out should be an exact copy of s except that SPECTRA column is empty 3909 calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ; 3910 aoff->unsetSelection() ; 3911 ahot->unsetSelection() ; 3920 3912 asky->unsetSelection() ; 3921 3913 sel.reset() ; 3922 } 3914 iter->next() ; 3915 } 3916 delete iter ; 3917 s->table_ = torg ; 3918 s->attach() ; 3923 3919 3924 3920 // flux unit … … 4538 4534 4539 4535 Vector<Float> STMath::getSpectrumFromTime( double reftime, 4540 Vector<Double> &timeVec, 4541 const CountedPtr<Scantable>& s, 4536 const Vector<Double> &timeVec, 4537 const vector<int> &idx, 4538 const Matrix<Float>& spectra, 4542 4539 string mode ) 4543 4540 { 4544 4541 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ; 4545 4542 Vector<Float> sp ; 4546 4547 if ( s->nrow() == 0 ) { 4543 uInt ncol = spectra.ncolumn() ; 4544 4545 if ( ncol == 0 ) { 4548 4546 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ; 4549 4547 return sp ; 4550 4548 } 4551 else if ( s->nrow()== 1 ) {4549 else if ( ncol == 1 ) { 4552 4550 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ; 4553 s ->specCol_.get( 0, sp) ;4551 sp.reference( spectra.column( 0 ) ) ; 4554 4552 return sp ; 4555 4553 } 4556 4554 else { 4557 vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;4558 4555 if ( mode == "before" ) { 4559 4556 int id = -1 ; … … 4566 4563 } 4567 4564 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4568 s->specCol_.get( id, sp ) ; 4569 //sp = s->getSpectrum( id ) ; 4565 sp.reference( spectra.column( id ) ) ; 4570 4566 } 4571 4567 else if ( mode == "after" ) { … … 4579 4575 } 4580 4576 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4581 s->specCol_.get( id, sp ) ; 4582 //sp = s->getSpectrum( id ) ; 4577 sp.reference( spectra.column( id ) ) ; 4583 4578 } 4584 4579 else if ( mode == "nearest" ) { … … 4605 4600 } 4606 4601 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4607 s->specCol_.get( id, sp ) ; 4608 //sp = s->getSpectrum( id ) ; 4602 sp.reference( spectra.column( id ) ) ; 4609 4603 } 4610 4604 else if ( mode == "linear" ) { … … 4614 4608 int id = idx[1] ; 4615 4609 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4616 s->specCol_.get( id, sp ) ; 4617 //sp = s->getSpectrum( id ) ; 4610 sp.reference( spectra.column( id ) ) ; 4618 4611 } 4619 4612 else if ( idx[1] == -1 ) { … … 4622 4615 int id = idx[0] ; 4623 4616 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4624 s->specCol_.get( id, sp ) ; 4625 //sp = s->getSpectrum( id ) ; 4617 sp.reference( spectra.column( id ) ) ; 4626 4618 } 4627 4619 else if ( idx[0] == idx[1] ) { … … 4630 4622 int id = idx[0] ; 4631 4623 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 4632 s->specCol_.get( id, sp ) ; 4633 //sp = s->getSpectrum( id ) ; 4624 sp.reference( spectra.column( id ) ) ; 4634 4625 } 4635 4626 else { … … 4639 4630 double t1 = timeVec[idx[1]] ; 4640 4631 double tref = reftime ; 4641 s->specCol_.get( idx[0], sp ) ;4642 //sp = s->getSpectrum( idx[0] ) ;4643 //vector<float> sp1 = s->getSpectrum( idx[1] ) ;4644 Vector<Float> sp1 = s->specCol_( idx[1] ) ;4645 double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;4646 for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {4647 sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ;4648 }4649 }4650 }4651 else {4652 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;4653 }4654 return sp ;4655 }4656 }4657 4658 Vector<Float> STMath::getSpectrumFromTime( double reftime,4659 const Vector<Double> &timeVec,4660 const Matrix<Float>& spectra,4661 string mode )4662 {4663 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;4664 Vector<Float> sp ;4665 uInt ncol = spectra.ncolumn() ;4666 4667 if ( ncol == 0 ) {4668 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;4669 return sp ;4670 }4671 else if ( ncol == 1 ) {4672 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;4673 //s->specCol_.get( 0, sp ) ;4674 sp.reference( spectra.column( 0 ) ) ;4675 return sp ;4676 }4677 else {4678 vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;4679 if ( mode == "before" ) {4680 int id = -1 ;4681 if ( idx[0] != -1 ) {4682 id = idx[0] ;4683 }4684 else if ( idx[1] != -1 ) {4685 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;4686 id = idx[1] ;4687 }4688 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4689 //s->specCol_.get( id, sp ) ;4690 //sp = s->getSpectrum( id ) ;4691 sp.reference( spectra.column( id ) ) ;4692 }4693 else if ( mode == "after" ) {4694 int id = -1 ;4695 if ( idx[1] != -1 ) {4696 id = idx[1] ;4697 }4698 else if ( idx[0] != -1 ) {4699 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;4700 id = idx[1] ;4701 }4702 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4703 //s->specCol_.get( id, sp ) ;4704 //sp = s->getSpectrum( id ) ;4705 sp.reference( spectra.column( id ) ) ;4706 }4707 else if ( mode == "nearest" ) {4708 int id = -1 ;4709 if ( idx[0] == -1 ) {4710 id = idx[1] ;4711 }4712 else if ( idx[1] == -1 ) {4713 id = idx[0] ;4714 }4715 else if ( idx[0] == idx[1] ) {4716 id = idx[0] ;4717 }4718 else {4719 double t0 = timeVec[idx[0]] ;4720 double t1 = timeVec[idx[1]] ;4721 double tref = reftime ;4722 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {4723 id = idx[1] ;4724 }4725 else {4726 id = idx[0] ;4727 }4728 }4729 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4730 //s->specCol_.get( id, sp ) ;4731 //sp = s->getSpectrum( id ) ;4732 sp.reference( spectra.column( id ) ) ;4733 }4734 else if ( mode == "linear" ) {4735 if ( idx[0] == -1 ) {4736 // use after4737 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;4738 int id = idx[1] ;4739 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4740 //s->specCol_.get( id, sp ) ;4741 //sp = s->getSpectrum( id ) ;4742 sp.reference( spectra.column( id ) ) ;4743 }4744 else if ( idx[1] == -1 ) {4745 // use before4746 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;4747 int id = idx[0] ;4748 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4749 //s->specCol_.get( id, sp ) ;4750 //sp = s->getSpectrum( id ) ;4751 sp.reference( spectra.column( id ) ) ;4752 }4753 else if ( idx[0] == idx[1] ) {4754 // use before4755 //os << "No need to interporate." << LogIO::POST ;4756 int id = idx[0] ;4757 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4758 //s->specCol_.get( id, sp ) ;4759 //sp = s->getSpectrum( id ) ;4760 sp.reference( spectra.column( id ) ) ;4761 }4762 else {4763 // do interpolation4764 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;4765 double t0 = timeVec[idx[0]] ;4766 double t1 = timeVec[idx[1]] ;4767 double tref = reftime ;4768 //s->specCol_.get( idx[0], sp ) ;4769 //sp = s->getSpectrum( idx[0] ) ;4770 //vector<float> sp1 = s->getSpectrum( idx[1] ) ;4771 4632 sp = spectra.column( idx[0] ).copy() ; 4772 //Vector<Float> sp1 = s->specCol_( idx[1] ) ;4773 4633 Vector<Float> sp1( spectra.column( idx[1] ) ) ; 4774 4634 double tfactor = ( tref - t0 ) / ( t1 - t0 ) ; … … 4782 4642 } 4783 4643 return sp ; 4784 }4785 }4786 4787 void STMath::getSpectrumFromTime2( Vector<Float> &sp,4788 double reftime,4789 Vector<Double> &timeVec,4790 const CountedPtr<Scantable>& s,4791 string mode )4792 {4793 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;4794 //Vector<Float> sp ;4795 4796 if ( s->nrow() == 0 ) {4797 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;4798 //return sp ;4799 }4800 else if ( s->nrow() == 1 ) {4801 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;4802 //return s->getSpectrum( 0 ) ;4803 s->specCol_.get( 0, sp ) ;4804 //return sp ;4805 }4806 else {4807 vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;4808 if ( mode == "before" ) {4809 int id = -1 ;4810 if ( idx[0] != -1 ) {4811 id = idx[0] ;4812 }4813 else if ( idx[1] != -1 ) {4814 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;4815 id = idx[1] ;4816 }4817 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4818 //sp = s->getSpectrum( id ) ;4819 s->specCol_.get( id, sp ) ;4820 }4821 else if ( mode == "after" ) {4822 int id = -1 ;4823 if ( idx[1] != -1 ) {4824 id = idx[1] ;4825 }4826 else if ( idx[0] != -1 ) {4827 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;4828 id = idx[1] ;4829 }4830 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4831 //sp = s->getSpectrum( id ) ;4832 s->specCol_.get( id, sp ) ;4833 }4834 else if ( mode == "nearest" ) {4835 int id = -1 ;4836 if ( idx[0] == -1 ) {4837 id = idx[1] ;4838 }4839 else if ( idx[1] == -1 ) {4840 id = idx[0] ;4841 }4842 else if ( idx[0] == idx[1] ) {4843 id = idx[0] ;4844 }4845 else {4846 double t0 = timeVec[idx[0]] ;4847 double t1 = timeVec[idx[1]] ;4848 double tref = reftime ;4849 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {4850 id = idx[1] ;4851 }4852 else {4853 id = idx[0] ;4854 }4855 }4856 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4857 //sp = s->getSpectrum( id ) ;4858 s->specCol_.get( id, sp ) ;4859 }4860 else if ( mode == "linear" ) {4861 if ( idx[0] == -1 ) {4862 // use after4863 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;4864 int id = idx[1] ;4865 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4866 //sp = s->getSpectrum( id ) ;4867 s->specCol_.get( id, sp ) ;4868 }4869 else if ( idx[1] == -1 ) {4870 // use before4871 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;4872 int id = idx[0] ;4873 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4874 //sp = s->getSpectrum( id ) ;4875 s->specCol_.get( id, sp ) ;4876 }4877 else if ( idx[0] == idx[1] ) {4878 // use before4879 //os << "No need to interporate." << LogIO::POST ;4880 int id = idx[0] ;4881 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4882 //sp = s->getSpectrum( id ) ;4883 s->specCol_.get( id, sp ) ;4884 }4885 else {4886 // do interpolation4887 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;4888 //double t0 = timeVec[idx[0]] ;4889 //double t1 = timeVec[idx[1]] ;4890 //double tref = reftime ;4891 //sp = s->getSpectrum( idx[0] ) ;4892 s->specCol_.get( idx[0], sp ) ;4893 //vector<float> sp1 = s->getSpectrum( idx[1] ) ;4894 Vector<Float> sp1 = s->specCol_( idx[1] ) ;4895 //double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;4896 double tfactor = ( reftime - timeVec[idx[0]] ) / ( timeVec[idx[1]] - timeVec[idx[0]] ) ;4897 for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {4898 sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ;4899 }4900 }4901 }4902 else {4903 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;4904 }4905 //return sp ;4906 4644 } 4907 4645 } … … 5153 4891 } 5154 4892 4893 Vector<Float> STMath::getTcalFromTime( double reftime, 4894 const Vector<Double> &timeVec, 4895 const vector<int> &idx, 4896 const CountedPtr<Scantable>& s, 4897 string mode ) 4898 { 4899 LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ; 4900 STTcal tcalTable = s->tcal() ; 4901 String time ; 4902 Vector<Float> tcalval ; 4903 if ( s->nrow() == 0 ) { 4904 os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ; 4905 return tcalval ; 4906 } 4907 else if ( s->nrow() == 1 ) { 4908 uInt tcalid = s->getTcalId( 0 ) ; 4909 //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4910 tcalTable.getEntry( time, tcalval, tcalid ) ; 4911 return tcalval ; 4912 } 4913 else { 4914 if ( mode == "before" ) { 4915 int id = -1 ; 4916 if ( idx[0] != -1 ) { 4917 id = idx[0] ; 4918 } 4919 else if ( idx[1] != -1 ) { 4920 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ; 4921 id = idx[1] ; 4922 } 4923 uInt tcalid = s->getTcalId( id ) ; 4924 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4925 tcalTable.getEntry( time, tcalval, tcalid ) ; 4926 } 4927 else if ( mode == "after" ) { 4928 int id = -1 ; 4929 if ( idx[1] != -1 ) { 4930 id = idx[1] ; 4931 } 4932 else if ( idx[0] != -1 ) { 4933 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ; 4934 id = idx[1] ; 4935 } 4936 uInt tcalid = s->getTcalId( id ) ; 4937 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4938 tcalTable.getEntry( time, tcalval, tcalid ) ; 4939 } 4940 else if ( mode == "nearest" ) { 4941 int id = -1 ; 4942 if ( idx[0] == -1 ) { 4943 id = idx[1] ; 4944 } 4945 else if ( idx[1] == -1 ) { 4946 id = idx[0] ; 4947 } 4948 else if ( idx[0] == idx[1] ) { 4949 id = idx[0] ; 4950 } 4951 else { 4952 double t0 = timeVec[idx[0]] ; 4953 double t1 = timeVec[idx[1]] ; 4954 if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) { 4955 id = idx[1] ; 4956 } 4957 else { 4958 id = idx[0] ; 4959 } 4960 } 4961 uInt tcalid = s->getTcalId( id ) ; 4962 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4963 tcalTable.getEntry( time, tcalval, tcalid ) ; 4964 } 4965 else if ( mode == "linear" ) { 4966 if ( idx[0] == -1 ) { 4967 // use after 4968 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ; 4969 int id = idx[1] ; 4970 uInt tcalid = s->getTcalId( id ) ; 4971 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4972 tcalTable.getEntry( time, tcalval, tcalid ) ; 4973 } 4974 else if ( idx[1] == -1 ) { 4975 // use before 4976 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ; 4977 int id = idx[0] ; 4978 uInt tcalid = s->getTcalId( id ) ; 4979 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4980 tcalTable.getEntry( time, tcalval, tcalid ) ; 4981 } 4982 else if ( idx[0] == idx[1] ) { 4983 // use before 4984 //os << "No need to interporate." << LogIO::POST ; 4985 int id = idx[0] ; 4986 uInt tcalid = s->getTcalId( id ) ; 4987 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4988 tcalTable.getEntry( time, tcalval, tcalid ) ; 4989 } 4990 else { 4991 // do interpolation 4992 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ; 4993 double t0 = timeVec[idx[0]] ; 4994 double t1 = timeVec[idx[1]] ; 4995 Vector<Float> tcal0 ; 4996 uInt tcalid0 = s->getTcalId( idx[0] ) ; 4997 uInt tcalid1 = s->getTcalId( idx[1] ) ; 4998 tcalTable.getEntry( time, tcal0, tcalid0 ) ; 4999 tcalTable.getEntry( time, tcalval, tcalid1 ) ; 5000 double tfactor = (reftime - t0) / (t1 - t0) ; 5001 for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) { 5002 tcalval[i] = ( tcalval[i] - tcal0[i] ) * tfactor + tcal0[i] ; 5003 } 5004 } 5005 } 5006 else { 5007 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ; 5008 } 5009 return tcalval ; 5010 } 5011 } 5012 5155 5013 vector<float> STMath::getTsysFromTime( string reftime, 5156 5014 CountedPtr<Scantable>& s, … … 5267 5125 tsysval.tovector( tsys0 ) ; 5268 5126 tsysval = tsysCol( idx[1] ) ; 5269 tsysval.tovector( tsys1 ) ; 5127 tsysval.tovector( tsys1 ) ; 5270 5128 for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) { 5271 5129 float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ; … … 5278 5136 } 5279 5137 return tsys ; 5138 } 5139 } 5140 5141 Vector<Float> STMath::getTsysFromTime( double reftime, 5142 const Vector<Double> &timeVec, 5143 const vector<int> &idx, 5144 const CountedPtr<Scantable> &s, 5145 string mode ) 5146 { 5147 LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ; 5148 ArrayColumn<Float> tsysCol ; 5149 tsysCol.attach( s->table(), "TSYS" ) ; 5150 Vector<Float> tsysval ; 5151 if ( s->nrow() == 0 ) { 5152 os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ; 5153 return tsysval ; 5154 } 5155 else if ( s->nrow() == 1 ) { 5156 //os << "use row " << 0 << LogIO::POST ; 5157 tsysval = tsysCol( 0 ) ; 5158 return tsysval ; 5159 } 5160 else { 5161 if ( mode == "before" ) { 5162 int id = -1 ; 5163 if ( idx[0] != -1 ) { 5164 id = idx[0] ; 5165 } 5166 else if ( idx[1] != -1 ) { 5167 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ; 5168 id = idx[1] ; 5169 } 5170 //os << "use row " << id << LogIO::POST ; 5171 tsysval = tsysCol( id ) ; 5172 } 5173 else if ( mode == "after" ) { 5174 int id = -1 ; 5175 if ( idx[1] != -1 ) { 5176 id = idx[1] ; 5177 } 5178 else if ( idx[0] != -1 ) { 5179 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ; 5180 id = idx[1] ; 5181 } 5182 //os << "use row " << id << LogIO::POST ; 5183 tsysval = tsysCol( id ) ; 5184 } 5185 else if ( mode == "nearest" ) { 5186 int id = -1 ; 5187 if ( idx[0] == -1 ) { 5188 id = idx[1] ; 5189 } 5190 else if ( idx[1] == -1 ) { 5191 id = idx[0] ; 5192 } 5193 else if ( idx[0] == idx[1] ) { 5194 id = idx[0] ; 5195 } 5196 else { 5197 double t0 = timeVec[idx[0]] ; 5198 double t1 = timeVec[idx[1]] ; 5199 if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) { 5200 id = idx[1] ; 5201 } 5202 else { 5203 id = idx[0] ; 5204 } 5205 } 5206 //os << "use row " << id << LogIO::POST ; 5207 tsysval = tsysCol( id ) ; 5208 } 5209 else if ( mode == "linear" ) { 5210 if ( idx[0] == -1 ) { 5211 // use after 5212 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ; 5213 int id = idx[1] ; 5214 //os << "use row " << id << LogIO::POST ; 5215 tsysval = tsysCol( id ) ; 5216 } 5217 else if ( idx[1] == -1 ) { 5218 // use before 5219 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ; 5220 int id = idx[0] ; 5221 //os << "use row " << id << LogIO::POST ; 5222 tsysval = tsysCol( id ) ; 5223 } 5224 else if ( idx[0] == idx[1] ) { 5225 // use before 5226 //os << "No need to interporate." << LogIO::POST ; 5227 int id = idx[0] ; 5228 //os << "use row " << id << LogIO::POST ; 5229 tsysval = tsysCol( id ) ; 5230 } 5231 else { 5232 // do interpolation 5233 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ; 5234 double t0 = timeVec[idx[0]] ; 5235 double t1 = timeVec[idx[1]] ; 5236 Vector<Float> tsys0 ; 5237 tsys0 = tsysCol( idx[0] ) ; 5238 tsysval = tsysCol( idx[1] ) ; 5239 double tfactor = (reftime - t0) / (t1 - t0) ; 5240 for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) { 5241 tsysval[i] = ( tsysval[i] - tsys0[i] ) * tfactor + tsys0[i] ; 5242 } 5243 } 5244 } 5245 else { 5246 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ; 5247 } 5248 return tsysval ; 5249 } 5250 } 5251 5252 void STMath::calibrateCW( CountedPtr<Scantable> &out, 5253 const CountedPtr<Scantable>& on, 5254 const CountedPtr<Scantable>& off, 5255 const CountedPtr<Scantable>& sky, 5256 const CountedPtr<Scantable>& hot, 5257 const CountedPtr<Scantable>& cold, 5258 const Vector<uInt> &rows, 5259 const String &antname ) 5260 { 5261 // 2012/05/22 TN 5262 // Assume that out has empty SPECTRA column 5263 5264 // if rows is empty, just return 5265 if ( rows.nelements() == 0 ) 5266 return ; 5267 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ; 5268 Vector<Double> timeOff = timeCol.getColumn() ; 5269 timeCol.attach( sky->table(), "TIME" ) ; 5270 Vector<Double> timeSky = timeCol.getColumn() ; 5271 timeCol.attach( hot->table(), "TIME" ) ; 5272 Vector<Double> timeHot = timeCol.getColumn() ; 5273 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" ) ; 5282 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ; 5283 // I know that the data is contiguous 5284 const uInt *p = rows.data() ; 5285 vector<int> ids( 2 ) ; 5286 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) { 5287 double reftime = timeCol.asdouble(*p) ; 5288 ids = getRowIdFromTime( reftime, timeOff ) ; 5289 Vector<Float> spoff = getSpectrumFromTime( reftime, timeOff, ids, offspectra, "linear" ) ; 5290 ids = getRowIdFromTime( reftime, timeSky ) ; 5291 Vector<Float> spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ; 5292 Vector<Float> tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ; 5293 Vector<Float> tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ; 5294 ids = getRowIdFromTime( reftime, timeHot ) ; 5295 Vector<Float> sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra, "linear" ) ; 5296 Vector<Float> spec = on->specCol_( *p ) ; 5297 if ( antname.find( "APEX" ) != String::npos ) { 5298 // using gain array 5299 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) { 5300 spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] ) 5301 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ; 5302 } 5303 } 5304 else { 5305 // Chopper-Wheel calibration (Ulich & Haas 1976) 5306 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) { 5307 spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ; 5308 } 5309 } 5310 out->specCol_.put( *p, spec ) ; 5311 out->tsysCol_.put( *p, tsys ) ; 5312 p++ ; 5280 5313 } 5281 5314 } … … 5333 5366 } 5334 5367 5335 vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,5336 CountedPtr<Scantable>& off,5337 int index )5338 {5339 // string reftime = on->getTime( index ) ;5340 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;5341 Vector<Double> timeVec = timeCol.getColumn() ;5342 //ROTableColumn timeCol( on->table(), "TIME" ) ;5343 timeCol.attach( on->table(), "TIME" ) ;5344 double reftime = timeCol.asdouble(index) ;5345 vector<int> ii( 1, on->getIF( index ) ) ;5346 vector<int> ib( 1, on->getBeam( index ) ) ;5347 vector<int> ip( 1, on->getPol( index ) ) ;5348 STSelector sel = STSelector() ;5349 sel.setIFs( ii ) ;5350 sel.setBeams( ib ) ;5351 sel.setPolarizations( ip ) ;5352 off->setSelection( sel ) ;5353 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;5354 vector<float> spec = on->getSpectrum( index ) ;5355 //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;5356 //vector<float> tsys = on->getTsysVec( index ) ;5357 ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;5358 Vector<Float> tsys = tsysCol( index ) ;5359 vector<float> sp( spec.size() ) ;5360 // ALMA Calibration5361 //5362 // Ta* = Tsys * ( ON - OFF ) / OFF5363 //5364 // 2010/01/07 Takeshi Nakazato5365 unsigned int tsyssize = tsys.nelements() ;5366 unsigned int spsize = sp.size() ;5367 for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {5368 float tscale = 0.0 ;5369 if ( tsyssize == spsize )5370 tscale = tsys[j] ;5371 else5372 tscale = tsys[0] ;5373 float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;5374 sp[j] = v ;5375 }5376 sel.reset() ;5377 off->unsetSelection() ;5378 5379 return sp ;5380 }5381 5382 5368 void STMath::calibrateALMA( CountedPtr<Scantable>& out, 5383 5369 const CountedPtr<Scantable>& on, 5384 5370 const CountedPtr<Scantable>& off, 5385 Vector<uInt>& rows )5371 const Vector<uInt>& rows ) 5386 5372 { 5387 5373 // 2012/05/22 TN … … 5401 5387 // I know that the data is contiguous 5402 5388 const uInt *p = rows.data() ; 5389 vector<int> ids( 2 ) ; 5403 5390 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) { 5404 5391 double reftime = timeCol.asdouble(*p) ; 5405 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, offspectra, "linear" ) ; 5392 ids = getRowIdFromTime( reftime, timeVec ) ; 5393 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, ids, offspectra, "linear" ) ; 5406 5394 //Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ; 5407 5395 Vector<Float> spec = on->specCol_( *p ) ; -
branches/hpc33/src/STMath.h
r2560 r2563 403 403 404 404 vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ; 405 casa::Vector<casa::Float> getSpectrumFromTime( double reftime, casa::Vector<casa::Double> &v, const casa::CountedPtr<Scantable>& s, string mode = "before" ) ; 406 casa::Vector<casa::Float> getSpectrumFromTime( double reftime, const casa::Vector<casa::Double> &v, const casa::Matrix<casa::Float>& sp, string mode = "before" ) ; 407 void getSpectrumFromTime2( casa::Vector<casa::Float> &sp, double reftime, casa::Vector<casa::Double> &v, const casa::CountedPtr<Scantable>& s, string mode = "before" ) ; 405 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" ) ; 408 406 vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ; 407 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" ) ; 409 408 vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ; 409 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 410 vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ; 411 411 vector<int> getRowIdFromTime( double reftime, const casa::Vector<casa::Double>& t ) ; 412 412 413 413 // Chopper-Wheel type calibration 414 void calibrateCW( casa::CountedPtr<Scantable> &out, 415 const casa::CountedPtr<Scantable> &on, 416 const casa::CountedPtr<Scantable> &off, 417 const casa::CountedPtr<Scantable> &sky, 418 const casa::CountedPtr<Scantable> &hot, 419 const casa::CountedPtr<Scantable> &cold, 420 const casa::Vector<casa::uInt> &rows, 421 const casa::String &antname ) ; 414 422 vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on, 415 423 casa::CountedPtr<Scantable>& off, … … 420 428 string antname ) ; 421 429 // Tsys * (ON-OFF)/OFF 422 vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,423 casa::CountedPtr<Scantable>& off,424 int index ) ;425 430 void calibrateALMA( casa::CountedPtr<Scantable>& out, 426 431 const casa::CountedPtr<Scantable>& on, 427 432 const casa::CountedPtr<Scantable>& off, 428 c asa::Vector<casa::uInt>& rows ) ;433 const casa::Vector<casa::uInt>& rows ) ; 429 434 vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig, 430 435 casa::CountedPtr<Scantable>& ref,
Note:
See TracChangeset
for help on using the changeset viewer.