Changeset 2563 for branches/hpc33/src


Ignore:
Timestamp:
06/12/12 18:45:52 (12 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Rewrote cwcal based on update of almacal.


Location:
branches/hpc33/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/hpc33/src/STMath.cpp

    r2561 r2563  
    38433843    vector<int> types ;
    38443844
     3845    // save original table selection
     3846    Table torg  = s->table_ ;
     3847
    38453848    // 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" ) ;
    38563859    // 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" ) ;
    38673865    // cold scan
    38683866    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" ) ;
    38783872
    38793873    // 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" ) ;
    38893879   
    38903880    // 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 ) ;
    39023887   
    39033888    // 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() ;
    39203912      asky->unsetSelection() ;
    39213913      sel.reset() ;
    3922     }
     3914      iter->next() ;
     3915    }
     3916    delete iter ;
     3917    s->table_ = torg ;
     3918    s->attach() ;
    39233919
    39243920    // flux unit
    … …  
    45384534
    45394535Vector<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,
    45424539                                           string mode )
    45434540{
    45444541  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
    45454542  Vector<Float> sp ;
    4546 
    4547   if ( s->nrow() == 0 ) {
     4543  uInt ncol = spectra.ncolumn() ;
     4544
     4545  if ( ncol == 0 ) {
    45484546    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
    45494547    return sp ;
    45504548  }
    4551   else if ( s->nrow() == 1 ) {
     4549  else if ( ncol == 1 ) {
    45524550    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
    4553     s->specCol_.get( 0, sp ) ;
     4551    sp.reference( spectra.column( 0 ) ) ;
    45544552    return sp ;
    45554553  }
    45564554  else {
    4557     vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
    45584555    if ( mode == "before" ) {
    45594556      int id = -1 ;
    … …  
    45664563      }
    45674564      //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 ) ) ;
    45704566    }
    45714567    else if ( mode == "after" ) {
    … …  
    45794575      }
    45804576      //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 ) ) ;
    45834578    }
    45844579    else if ( mode == "nearest" ) {
    … …  
    46054600      }
    46064601      //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 ) ) ;
    46094603    }
    46104604    else if ( mode == "linear" ) {
    … …  
    46144608        int id = idx[1] ;
    46154609        //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 ) ) ;
    46184611      }
    46194612      else if ( idx[1] == -1 ) {
    … …  
    46224615        int id = idx[0] ;
    46234616        //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 ) ) ;
    46264618      }
    46274619      else if ( idx[0] == idx[1] ) {
    … …  
    46304622        int id = idx[0] ;
    46314623        //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 ) ) ;
    46344625      }
    46354626      else {
    … …  
    46394630        double t1 = timeVec[idx[1]] ;
    46404631        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 after
    4737         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 before
    4746         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 before
    4755         //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 interpolation
    4764         //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] ) ;
    47714632        sp = spectra.column( idx[0] ).copy() ;
    4772         //Vector<Float> sp1 = s->specCol_( idx[1] ) ;
    47734633        Vector<Float> sp1( spectra.column( idx[1] ) ) ;
    47744634        double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;
    … …  
    47824642    }
    47834643    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 after
    4863         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 before
    4871         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 before
    4879         //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 interpolation
    4887         //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 ;
    49064644  }
    49074645}
    … …  
    51534891}
    51544892
     4893Vector<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
    51555013vector<float> STMath::getTsysFromTime( string reftime,
    51565014                                       CountedPtr<Scantable>& s,
    … …  
    52675125        tsysval.tovector( tsys0 ) ;
    52685126        tsysval = tsysCol( idx[1] ) ;
    5269         tsysval.tovector( tsys1 ) ;       
     5127        tsysval.tovector( tsys1 ) ;
    52705128        for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {
    52715129          float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ;
    … …  
    52785136    }
    52795137    return tsys ;
     5138  }
     5139}
     5140
     5141Vector<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
     5252void 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++ ;
    52805313  }
    52815314}
    … …  
    53335366}
    53345367
    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 Calibration
    5361   //
    5362   // Ta* = Tsys * ( ON - OFF ) / OFF
    5363   //
    5364   // 2010/01/07 Takeshi Nakazato
    5365   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     else
    5372       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 
    53825368void STMath::calibrateALMA( CountedPtr<Scantable>& out,
    53835369                            const CountedPtr<Scantable>& on,
    53845370                            const CountedPtr<Scantable>& off,
    5385                             Vector<uInt>& rows )
     5371                            const Vector<uInt>& rows )
    53865372{
    53875373  // 2012/05/22 TN
    … …  
    54015387  // I know that the data is contiguous
    54025388  const uInt *p = rows.data() ;
     5389  vector<int> ids( 2 ) ;
    54035390  for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    54045391    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" ) ;
    54065394    //Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
    54075395    Vector<Float> spec = on->specCol_( *p ) ;

  • branches/hpc33/src/STMath.h

    r2560 r2563  
    403403
    404404  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" ) ;
    408406  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" ) ;
    409408  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" ) ;
    410410  vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;
    411411  vector<int> getRowIdFromTime( double reftime, const casa::Vector<casa::Double>& t ) ;
    412412
    413413  // 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 ) ;
    414422  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
    415423                                      casa::CountedPtr<Scantable>& off,
    … …  
    420428                                      string antname ) ;
    421429  // Tsys * (ON-OFF)/OFF
    422   vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
    423                                       casa::CountedPtr<Scantable>& off,
    424                                       int index ) ;
    425430  void calibrateALMA( casa::CountedPtr<Scantable>& out,
    426431                      const casa::CountedPtr<Scantable>& on,
    427432                      const casa::CountedPtr<Scantable>& off,
    428                       casa::Vector<casa::uInt>& rows ) ;
     433                      const casa::Vector<casa::uInt>& rows ) ;
    429434  vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
    430435                                        casa::CountedPtr<Scantable>& ref,
Note: See TracChangeset for help on using the changeset viewer.