Changeset 2919 for trunk/src/STMath.cpp


Ignore:
Timestamp:
04/04/14 20:08:17 (10 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...

Refactoring.


File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STMath.cpp

    r2918 r2919  
    5858
    5959#include "CalibrationHelper.h"
     60#include "IterationHelper.h"
    6061
    6162using namespace casa;
    6263using namespace asap;
    63 
    64 // namespace {
    65 // class CalibrationIterator
    66 // {
    67 // public:
    68 //   CalibrationIterator()
    69 //     : table_(table)
    70 //   {}
    71 //   ~CalibrationIterator() {}
    72 //   void Iterate() {
    73 //     STSelector sel;
    74 //     vector<string> cols( 3 ) ;
    75 //     cols[0] = "BEAMNO" ;
    76 //     cols[1] = "POLNO" ;
    77 //     cols[2] = "IFNO" ;
    78 //     STIdxIter2 iter(out, cols);
    79 //     while ( !iter.pastEnd() ) {
    80 //       Record ids = iter.currentValue();
    81 //       stringstream ss ;
    82 //       ss << "SELECT FROM $1 WHERE "
    83 //          << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
    84 //          << "POLNO==" << ids.asuInt("POLNO") << "&&"
    85 //          << "IFNO==" << ids.asuInt("IFNO") ;
    86 //       //cout << "TaQL string: " << ss.str() << endl ;
    87 //       sel.setTaQL( ss.str() ) ;
    88 //       ////
    89 //       // begin individual processing
    90 //       aoff->setSelection( sel ) ;
    91 //       ahot->setSelection( sel ) ;
    92 //       asky->setSelection( sel ) ;
    93 //       Vector<uInt> rows = iter.getRows( SHARE ) ;
    94 //       // out should be an exact copy of s except that SPECTRA column is empty
    95 //       calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ;
    96 //       aoff->unsetSelection() ;
    97 //       ahot->unsetSelection() ;
    98 //       asky->unsetSelection() ;
    99 //       // end individual processing
    100 //       ////
    101 //       sel.reset() ;
    102 //       iter.next() ;
    103 //     }   
    104 //   };
    105 // protected:
    106 //   CountedPtr<Scantable> table_;
    107 // };
    108 
    109 
    110 // } // anonymous namespace
    11164
    11265// 2012/02/17 TN
    … …  
    35813534    copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ;
    35823535   
    3583     // process each on scan
    3584     STSelector sel ;
    3585     vector<string> cols( 3 ) ;
    3586     cols[0] = "BEAMNO" ;
    3587     cols[1] = "POLNO" ;
    3588     cols[2] = "IFNO" ;
    3589     STIdxIter2 iter(out, cols);
    3590     while ( !iter.pastEnd() ) {
    3591       Record ids = iter.currentValue();
    3592       stringstream ss ;
    3593       ss << "SELECT FROM $1 WHERE "
    3594          << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
    3595          << "POLNO==" << ids.asuInt("POLNO") << "&&"
    3596          << "IFNO==" << ids.asuInt("IFNO") ;
    3597       //cout << "TaQL string: " << ss.str() << endl ;
    3598       sel.setTaQL( ss.str() ) ;
    3599       aoff->setSelection( sel ) ;
    3600       ahot->setSelection( sel ) ;
    3601       asky->setSelection( sel ) ;
    3602       Vector<uInt> rows = iter.getRows( SHARE ) ;
    3603       // out should be an exact copy of s except that SPECTRA column is empty
    3604       calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ;
    3605       aoff->unsetSelection() ;
    3606       ahot->unsetSelection() ;
    3607       asky->unsetSelection() ;
    3608       sel.reset() ;
    3609       iter.next() ;
    3610     }
     3536    // iterate throgh IterationHelper
     3537    ChopperWheelCalibrator cal(out, s, aoff, asky, ahot, acold);
     3538    IterationHelper<ChopperWheelCalibrator>::Iterate(cal, "BEAMNO,POLNO,IFNO");
     3539
    36113540    s->table_ = torg ;
    36123541    s->attach() ;
    … …  
    36643593//     t0 = mathutil::gettimeofday_sec() ;
    36653594
    3666     // using STIdxIter2
    3667     vector<string> cols( 3 ) ;
    3668     cols[0] = "BEAMNO" ;
    3669     cols[1] = "POLNO" ;
    3670     cols[2] = "IFNO" ;
    3671     STIdxIter2 iter( out, cols ) ;
    3672     STSelector sel ;
    3673     while ( !iter.pastEnd() ) {
    3674       Record ids = iter.currentValue() ;
    3675       stringstream ss ;
    3676       ss << "SELECT FROM $1 WHERE "
    3677          << "BEAMNO==" << ids.asuInt(cols[0]) << "&&"
    3678          << "POLNO==" << ids.asuInt(cols[1]) << "&&"
    3679          << "IFNO==" << ids.asuInt(cols[2]) ;
    3680       //cout << "TaQL string: " << ss.str() << endl ;
    3681       sel.setTaQL( ss.str() ) ;
    3682       aoff->setSelection( sel ) ;
    3683       Vector<uInt> rows = iter.getRows( SHARE ) ;
    3684       // out should be an exact copy of s except that SPECTRA column is empty
    3685       calibrateALMA( out, s, aoff, rows ) ;
    3686       aoff->unsetSelection() ;
    3687       sel.reset() ;
    3688       iter.next() ;
    3689     }
     3595    // iterate throgh IterationHelper
     3596    AlmaCalibrator cal(out, s, aoff);
     3597    IterationHelper<AlmaCalibrator>::Iterate(cal, "BEAMNO,POLNO,IFNO");
     3598
     3599    // finalize
    36903600    s->table_ = torg ;
    36913601    s->attach() ;
    … …  
    38483758     
    38493759      // process each sig and ref scan
    3850 //       STSelector sel ;
    3851       vector<string> cols( 3 ) ;
    3852       cols[0] = "BEAMNO" ;
    3853       cols[1] = "POLNO" ;
    3854       cols[2] = "IFNO" ;
    3855       STIdxIter2 iter(ssig, cols);
    3856       STSelector sel ;
    3857       while ( !iter.pastEnd() ) {
    3858         Record ids = iter.currentValue() ;
    3859         stringstream ss ;
    3860         ss << "SELECT FROM $1 WHERE "
    3861            << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
    3862            << "POLNO==" << ids.asuInt("POLNO") << "&&"
    3863            << "IFNO==" << ids.asuInt("IFNO") ;
    3864         //cout << "TaQL string: " << ss.str() << endl ;
    3865         sel.setTaQL( ss.str() ) ;
    3866         aofflo->setSelection( sel ) ;
    3867         ahotlo->setSelection( sel ) ;
    3868         askylo->setSelection( sel ) ;
    3869         Vector<uInt> rows = iter.getRows( SHARE ) ;
    3870         calibrateCW( ssig, rsig, aofflo, askylo, ahotlo, acoldlo, rows, antname ) ;
    3871         aofflo->unsetSelection() ;
    3872         ahotlo->unsetSelection() ;
    3873         askylo->unsetSelection() ;
    3874         sel.reset() ;
    3875         iter.next() ;
    3876       }
    3877       iter = STIdxIter2(sref, cols);
    3878       while ( !iter.pastEnd() ) {
    3879         Record ids = iter.currentValue() ;
    3880         stringstream ss ;
    3881         ss << "SELECT FROM $1 WHERE "
    3882            << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
    3883            << "POLNO==" << ids.asuInt("POLNO") << "&&"
    3884            << "IFNO==" << ids.asuInt("IFNO") ;
    3885         //cout << "TaQL string: " << ss.str() << endl ;
    3886         sel.setTaQL( ss.str() ) ;
    3887         aoffhi->setSelection( sel ) ;
    3888         ahothi->setSelection( sel ) ;
    3889         askyhi->setSelection( sel ) ;
    3890         Vector<uInt> rows = iter.getRows( SHARE ) ;
    3891         calibrateCW( sref, rref, aoffhi, askyhi, ahothi, acoldhi, rows, antname ) ;
    3892         aoffhi->unsetSelection() ;
    3893         ahothi->unsetSelection() ;
    3894         askyhi->unsetSelection() ;
    3895         sel.reset() ;
    3896         iter.next() ;
    3897       }
     3760      // iterate throgh IterationHelper
     3761      ChopperWheelCalibrator cal_sig(ssig, rsig, aofflo, askylo, ahotlo, acoldlo);
     3762      IterationHelper<ChopperWheelCalibrator>::Iterate(cal_sig, "BEAMNO,POLNO,IFNO");
     3763      ChopperWheelCalibrator cal_ref(sref, rref, aoffhi, askyhi, ahothi, acoldhi);
     3764      IterationHelper<ChopperWheelCalibrator>::Iterate(cal_ref, "BEAMNO,POLNO,IFNO");
    38983765    }
    38993766  }
    … …  
    40533920
    40543921  return out ;
    4055 }
    4056 
    4057 void STMath::calibrateCW( CountedPtr<Scantable> &out,
    4058                           const CountedPtr<Scantable>& on,
    4059                           const CountedPtr<Scantable>& off,
    4060                           const CountedPtr<Scantable>& sky,
    4061                           const CountedPtr<Scantable>& hot,
    4062                           const CountedPtr<Scantable>& cold,
    4063                           const Vector<uInt> &rows,
    4064                           const String &antname )
    4065 {
    4066   // 2012/05/22 TN
    4067   // Assume that out has empty SPECTRA column
    4068 
    4069   // if rows is empty, just return
    4070   if ( rows.nelements() == 0 )
    4071     return ;
    4072   ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
    4073   Vector<Double> timeOff = timeCol.getColumn() ;
    4074   timeCol.attach( sky->table(), "TIME" ) ;
    4075   Vector<Double> timeSky = timeCol.getColumn() ;
    4076   timeCol.attach( hot->table(), "TIME" ) ;
    4077   Vector<Double> timeHot = timeCol.getColumn() ;
    4078   timeCol.attach( on->table(), "TIME" ) ;
    4079   ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
    4080   SpectralData offspectra(arrayFloatCol.getColumn());
    4081   arrayFloatCol.attach( sky->table(), "SPECTRA" ) ;
    4082   SpectralData skyspectra(arrayFloatCol.getColumn());
    4083   arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
    4084   SpectralData hotspectra(arrayFloatCol.getColumn());
    4085   TcalData tcaldata(sky);
    4086   TsysData tsysdata(sky);
    4087   unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
    4088   // I know that the data is contiguous
    4089   const uInt *p = rows.data() ;
    4090   vector<int> ids( 2 ) ;
    4091   Block<uInt> flagchan( spsize ) ;
    4092   uInt nflag = 0 ;
    4093   for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    4094     double reftime = timeCol.asdouble(*p) ;
    4095     ids = getRowIdFromTime( reftime, timeOff ) ;
    4096     Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeOff, ids, offspectra, "linear");
    4097     ids = getRowIdFromTime( reftime, timeSky ) ;
    4098     Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
    4099     Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
    4100     Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
    4101     ids = getRowIdFromTime( reftime, timeHot ) ;
    4102     Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
    4103     Vector<Float> spec = on->specCol_( *p ) ;
    4104     if ( antname.find( "APEX" ) != String::npos ) {
    4105       // using gain array
    4106       for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    4107         if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
    4108           spec[j] = 0.0 ;
    4109           flagchan[nflag++] = j ;
    4110         }
    4111         else {
    4112           spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
    4113             * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
    4114         }
    4115       }
    4116     }
    4117     else {
    4118       // Chopper-Wheel calibration (Ulich & Haas 1976)
    4119       for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    4120         if ( (sphot[j]-spsky[j]) == 0.0 ) {
    4121           spec[j] = 0.0 ;
    4122           flagchan[nflag++] = j ;
    4123         }
    4124         else {
    4125           spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
    4126         }
    4127       }
    4128     }
    4129     out->specCol_.put( *p, spec ) ;
    4130     out->tsysCol_.put( *p, tsys ) ;
    4131     if ( nflag > 0 ) {
    4132       Vector<uChar> fl = out->flagsCol_( *p ) ;
    4133       for ( unsigned int j = 0 ; j < nflag ; j++ ) {
    4134         fl[flagchan[j]] = (uChar)True ;
    4135       }
    4136       out->flagsCol_.put( *p, fl ) ;
    4137     }
    4138     nflag = 0 ;
    4139     p++ ;
    4140   }
    4141 }
    4142 
    4143 void STMath::calibrateALMA( CountedPtr<Scantable>& out,
    4144                             const CountedPtr<Scantable>& on,
    4145                             const CountedPtr<Scantable>& off,
    4146                             const Vector<uInt>& rows )
    4147 {
    4148   // 2012/05/22 TN
    4149   // Assume that out has empty SPECTRA column
    4150 
    4151   // if rows is empty, just return
    4152   if ( rows.nelements() == 0 )
    4153     return ;
    4154   ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
    4155   Vector<Double> timeVec = timeCol.getColumn() ;
    4156   timeCol.attach( on->table(), "TIME" ) ;
    4157   ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
    4158   SpectralData offspectra(arrayFloatCol.getColumn());
    4159   unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
    4160   // I know that the data is contiguous
    4161   const uInt *p = rows.data() ;
    4162   vector<int> ids( 2 ) ;
    4163   Block<uInt> flagchan( spsize ) ;
    4164   uInt nflag = 0 ;
    4165   for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    4166     double reftime = timeCol.asdouble(*p) ;
    4167     ids = getRowIdFromTime( reftime, timeVec ) ;
    4168     Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeVec, ids, offspectra, "linear");
    4169     Vector<Float> spec = on->specCol_( *p ) ;
    4170     Vector<Float> tsys = on->tsysCol_( *p ) ;
    4171     // ALMA Calibration
    4172     //
    4173     // Ta* = Tsys * ( ON - OFF ) / OFF
    4174     //
    4175     // 2010/01/07 Takeshi Nakazato
    4176     unsigned int tsyssize = tsys.nelements() ;
    4177     for ( unsigned int j = 0 ; j < spsize ; j++ ) {
    4178       if ( spoff[j] == 0.0 ) {
    4179         spec[j] = 0.0 ;
    4180         flagchan[nflag++] = j ;
    4181       }
    4182       else {
    4183         spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
    4184       }
    4185       if ( tsyssize == spsize )
    4186         spec[j] *= tsys[j] ;
    4187       else
    4188         spec[j] *= tsys[0] ;
    4189     }
    4190     out->specCol_.put( *p, spec ) ;
    4191     if ( nflag > 0 ) {
    4192       Vector<uChar> fl = out->flagsCol_( *p ) ;
    4193       for ( unsigned int j = 0 ; j < nflag ; j++ ) {
    4194         fl[flagchan[j]] = (uChar)True ;
    4195       }
    4196       out->flagsCol_.put( *p, fl ) ;
    4197     }
    4198     nflag = 0 ;
    4199     p++ ;
    4200   }
    42013922}
    42023923
Note: See TracChangeset for help on using the changeset viewer.