Changeset 2978


Ignore:
Timestamp:
07/25/14 23:06:41 (10 years ago)
Author:
WataruKawasaki
Message:

New Development: No

JIRA Issue: Yes CAS-6584, CAS-6787

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: test_sdcal

Put in Release Notes:

Module(s): sd

Description: (1) modified STMath::almacal(), STMath::cwcal(), and the relevant functions so that flag information are properly handled in sdcal's 'ps'(for ALMA and APEX), 'otf', and 'otfraster' modes.

(2) found a bug in SimpleInterpolationHelper::GetFromTime() and fixed it (the bug was causing a serious issue reported in CAS-6787).


Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/CalibrationHelper.h

    r2919 r2978  
     1#include <iostream>
    12#include <vector>
    23
    … …  
    106107}
    107108 
     109vector<int> getRowIdFromTime2(double reftime,
     110                              const Vector<Double> &t,
     111                              const Vector<uInt> &flagrow,
     112                              const Matrix<uChar> &flagtra)
     113{
     114  unsigned int nchan = flagtra[0].nelements();
     115  vector<int> v(2*nchan);
     116
     117  for (unsigned int j = 0; j < nchan; ++j) {
     118    //   double reft = reftime ;
     119    double dtmin = 1.0e100 ;
     120    double dtmax = -1.0e100 ;
     121    //   vector<double> dt ;
     122    int just_before = -1 ;
     123    int just_after = -1 ;
     124    Vector<Double> dt = t - reftime ;
     125    for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
     126      if ( flagrow[i] > 0 ) continue;
     127      if ( flagtra.column(i)[j] == 1 << 7) continue;
     128
     129      if ( dt[i] > 0.0 ) {
     130        // after reftime
     131        if ( dt[i] < dtmin ) {
     132          just_after = i ;
     133          dtmin = dt[i] ;
     134        }
     135      }
     136      else if ( dt[i] < 0.0 ) {
     137        // before reftime
     138        if ( dt[i] > dtmax ) {
     139          just_before = i ;
     140          dtmax = dt[i] ;
     141        }
     142      }
     143      else {
     144        // just a reftime
     145        just_before = i ;
     146        just_after = i ;
     147        dtmax = 0 ;
     148        dtmin = 0 ;
     149        break ;
     150      }
     151    }
     152
     153    v[j*2]   = just_before ;
     154    v[j*2+1] = just_after ;
     155  }
     156 
     157  return v ;
     158}
     159 
    108160template<class T>
    109161class SimpleInterpolationHelper
    … …  
    116168                                   const string mode)
    117169  {
    118     Vector<Float> return_value;
     170    Vector<Float> return_value(idx.size()/2);
    119171    LogIO os_;
    120172    LogIO os( LogOrigin( "STMath", data.method_name(), WHERE ) ) ;
    … …  
    126178    }
    127179    else {
    128       if ( mode == "before" ) {
    129         int id = -1 ;
    130         if ( idx[0] != -1 ) {
    131           id = idx[0] ;
    132         }
    133         else if ( idx[1] != -1 ) {
    134           os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
    135           id = idx[1] ;
    136         }
     180      for (unsigned int i = 0; i < idx.size()/2; ++i) {
     181        unsigned int idx0 = 2*i;
     182        unsigned int idx1 = 2*i + 1;
     183        //no off data available. calibration impossible.
     184        if ( ( idx[idx0] == -1 ) && ( idx[idx1] == -1 ) ) continue;
     185
     186        if ( mode == "before" ) {
     187          int id = -1 ;
     188          if ( idx[idx0] != -1 ) {
     189            id = idx[idx0] ;
     190          }
     191          else if ( idx[idx1] != -1 ) {
     192            os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
     193            id = idx[idx1] ;
     194          }
    137195       
    138         return_value = data.GetEntry(id);
    139       }
    140       else if ( mode == "after" ) {
    141         int id = -1 ;
    142         if ( idx[1] != -1 ) {
    143           id = idx[1] ;
    144         }
    145         else if ( idx[0] != -1 ) {
    146           os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
    147           id = idx[1] ;
    148         }
     196          return_value[i] = data.GetEntry(id)[i];
     197        }
     198        else if ( mode == "after" ) {
     199          int id = -1 ;
     200          if ( idx[idx1] != -1 ) {
     201            id = idx[idx1] ;
     202          }
     203          else if ( idx[idx0] != -1 ) {
     204            os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
     205            id = idx[idx1] ;
     206          }
    149207       
    150         return_value = data.GetEntry(id);
    151       }
    152       else if ( mode == "nearest" ) {
    153         int id = -1 ;
    154         if ( idx[0] == -1 ) {
    155           id = idx[1] ;
    156         }
    157         else if ( idx[1] == -1 ) {
    158           id = idx[0] ;
    159         }
    160         else if ( idx[0] == idx[1] ) {
    161           id = idx[0] ;
    162         }
    163         else {
    164           double t0 = timeVec[idx[0]] ;
    165           double t1 = timeVec[idx[1]] ;
    166           if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) {
    167             id = idx[1] ;
     208          return_value[i] = data.GetEntry(id)[i];
     209        }
     210        else if ( mode == "nearest" ) {
     211          int id = -1 ;
     212          if ( idx[idx0] == -1 ) {
     213            id = idx[idx1] ;
     214          }
     215          else if ( idx[idx1] == -1 ) {
     216            id = idx[idx0] ;
     217          }
     218          else if ( idx[idx0] == idx[idx1] ) {
     219            id = idx[idx0] ;
    168220          }
    169221          else {
    170             id = idx[0] ;
    171           }
    172         }
    173         return_value = data.GetEntry(id);
    174       }
    175       else if ( mode == "linear" ) {
    176         if ( idx[0] == -1 ) {
    177           // use after
    178           os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
    179           int id = idx[1] ;
    180           return_value = data.GetEntry(id);
    181         }
    182         else if ( idx[1] == -1 ) {
    183           // use before
    184           os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
    185           int id = idx[0] ;
    186           return_value = data.GetEntry(id);
    187         }
    188         else if ( idx[0] == idx[1] ) {
    189           // use before
    190           //os << "No need to interporate." << LogIO::POST ;
    191           int id = idx[0] ;
    192           return_value = data.GetEntry(id);
    193         }
    194         else {
    195           // do interpolation
    196 
    197           double t0 = timeVec[idx[0]] ;
    198           double t1 = timeVec[idx[1]] ;
    199           Vector<Float> value0 = data.GetEntry(idx[0]);
    200           Vector<Float> value1 = data.GetEntry(idx[1]);
    201           double tfactor = (reftime - t0) / (t1 - t0) ;
    202           for ( unsigned int i = 0 ; i < value0.size() ; i++ ) {
    203             value1[i] = ( value1[i] - value0[i] ) * tfactor + value0[i] ;
    204           }
    205           return_value = value1;
    206         }
    207       }
    208       else {
    209         os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     222            double t0 = timeVec[idx[idx0]] ;
     223            double t1 = timeVec[idx[idx1]] ;
     224            if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) {
     225              id = idx[idx1] ;
     226            }
     227            else {
     228              id = idx[idx0] ;
     229            }
     230          }
     231          return_value[i] = data.GetEntry(id)[i];
     232        }
     233        else if ( mode == "linear" ) {
     234          if ( idx[idx0] == -1 ) {
     235            // use after
     236            os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
     237            int id = idx[idx1] ;
     238            return_value[i] = data.GetEntry(id)[i];
     239          }
     240          else if ( idx[idx1] == -1 ) {
     241            // use before
     242            os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
     243            int id = idx[idx0] ;
     244            return_value[i] = data.GetEntry(id)[i];
     245          }
     246          else if ( idx[idx0] == idx[idx1] ) {
     247            // use before
     248            //os << "No need to interporate." << LogIO::POST ;
     249            int id = idx[idx0] ;
     250            return_value[i] = data.GetEntry(id)[i];
     251          }
     252          else {
     253            // do interpolation
     254            double t0 = timeVec[idx[idx0]] ;
     255            double t1 = timeVec[idx[idx1]] ;
     256            Vector<Float> value0 = data.GetEntry(idx[idx0]);
     257            Vector<Float> value1 = data.GetEntry(idx[idx1]);
     258            double tfactor = (reftime - t0) / (t1 - t0) ;
     259            return_value[i] = ( value1[i] - value0[i] ) * tfactor + value0[i] ;
     260          }
     261        }
     262        else {
     263          os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     264        }
    210265      }
    211266    }
    … …  
    238293    Vector<Double> timeVec = GetScalarColumn<Double>(off->table(), "TIME");
    239294    Vector<Double> refTimeVec = GetScalarColumn<Double>(on->table(), "TIME");
     295    Vector<uInt> flagrowVec = GetScalarColumn<uInt>(off->table(), "FLAGROW");
     296    Vector<uInt> refFlagrowVec = GetScalarColumn<uInt>(on->table(), "FLAGROW");
     297    Matrix<uChar> flagtraMtx = GetArrayColumn<uChar>(off->table(), "FLAGTRA");
    240298    SpectralData offspectra(Matrix<Float>(GetArrayColumn<Float>(off->table(), "SPECTRA")));
    241299    unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
    242     vector<int> ids( 2 ) ;
    243     for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
     300    vector<int> ids( 2 * spsize ) ;
     301
     302    for ( unsigned int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    244303      uInt row = rows[irow];
    245304      double reftime = refTimeVec[row];
    246       ids = getRowIdFromTime( reftime, timeVec ) ;
     305      ids = getRowIdFromTime2( reftime, timeVec, flagrowVec, flagtraMtx ) ;
    247306      Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeVec, ids, offspectra, "linear");
    248307      Vector<Float> spec = in_spectra_column(row);
    249308      Vector<Float> tsys = in_tsys_column(row);
    250309      Vector<uChar> flag = in_flagtra_column(row);
     310
    251311      // ALMA Calibration
    252312      //
    … …  
    256316      unsigned int tsyssize = tsys.nelements() ;
    257317      for ( unsigned int j = 0 ; j < spsize ; j++ ) {
    258         if ( spoff[j] == 0.0 ) {
    259           spec[j] = 0.0 ;
    260           flag[j] = (uChar)True;
     318        //if there is no off data available for a channel, just flag the channel.(2014/7/18 WK)
     319        if ((ids[2*j] == -1)&&(ids[2*j+1] == -1)) {
     320          flag[j] = 1 << 7;
     321          continue;
    261322        }
    262         else {
    263           spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
     323
     324        if (refFlagrowVec[row] == 0) {
     325          if ( spoff[j] == 0.0 ) {
     326            spec[j] = 0.0 ;
     327            flag[j] = (uChar)True;
     328          }
     329          else {
     330            spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
     331          }
     332          spec[j] *= (tsyssize == spsize) ? tsys[j] : tsys[0];
    264333        }
    265         spec[j] *= (tsyssize == spsize) ? tsys[j] : tsys[0];
    266334      }
    267335      out_spectra_column.put(row, spec);
    … …  
    295363    Vector<Double> timeHot = GetScalarColumn<Double>(hot->table(), "TIME");
    296364    Vector<Double> timeOn = GetScalarColumn<Double>(on->table(), "TIME");
     365    Vector<uInt> flagrowOff = GetScalarColumn<uInt>(off->table(), "FLAGROW");
     366    Vector<uInt> flagrowSky = GetScalarColumn<uInt>(sky->table(), "FLAGROW");
     367    Vector<uInt> flagrowHot = GetScalarColumn<uInt>(hot->table(), "FLAGROW");
     368    Vector<uInt> flagrowOn = GetScalarColumn<uInt>(on->table(), "FLAGROW");
     369    Matrix<uChar> flagtraOff = GetArrayColumn<uChar>(off->table(), "FLAGTRA");
     370    Matrix<uChar> flagtraSky = GetArrayColumn<uChar>(sky->table(), "FLAGTRA");
     371    Matrix<uChar> flagtraHot = GetArrayColumn<uChar>(hot->table(), "FLAGTRA");
    297372    SpectralData offspectra(Matrix<Float>(GetArrayColumn<Float>(off->table(), "SPECTRA")));
    298373    SpectralData skyspectra(Matrix<Float>(GetArrayColumn<Float>(sky->table(), "SPECTRA")));
    … …  
    301376    TsysData tsysdata(sky);
    302377    unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
    303     vector<int> ids( 2 ) ;
    304     for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
     378    vector<int> idsOff( 2 * spsize ) ;
     379    vector<int> idsSky( 2 * spsize ) ;
     380    vector<int> idsHot( 2 * spsize ) ;
     381    for ( unsigned int irow = 0 ; irow < rows.nelements() ; irow++ ) {
    305382      uInt row = rows[irow];
    306383      double reftime = timeOn[row];
    307       ids = getRowIdFromTime( reftime, timeOff ) ;
    308       Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeOff, ids, offspectra, "linear");
    309       ids = getRowIdFromTime( reftime, timeSky ) ;
    310       Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
    311       Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
    312       Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
    313       ids = getRowIdFromTime( reftime, timeHot ) ;
    314       Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
     384      idsOff = getRowIdFromTime2( reftime, timeOff, flagrowOff, flagtraOff ) ;
     385      Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeOff, idsOff, offspectra, "linear");
     386      idsSky = getRowIdFromTime2( reftime, timeSky, flagrowSky, flagtraSky ) ;
     387      Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, idsSky, skyspectra, "linear");
     388      Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, idsSky, tcaldata, "linear");
     389      Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, idsSky, tsysdata, "linear");
     390      idsHot = getRowIdFromTime2( reftime, timeHot, flagrowHot, flagtraHot ) ;
     391      Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, idsHot, hotspectra, "linear");
    315392      Vector<Float> spec = in_spectra_column(row);
    316393      Vector<uChar> flag = in_flagtra_column(row);
    … …  
    318395        // using gain array
    319396        for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    320           if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
    321             spec[j] = 0.0 ;
     397          //if at least one of off/sky/hot data unavailable, just flag the channel.
     398          if (((idsOff[2*j] == -1)&&(idsOff[2*j+1] == -1))||
     399              ((idsSky[2*j] == -1)&&(idsSky[2*j+1] == -1))||
     400              ((idsHot[2*j] == -1)&&(idsHot[2*j+1] == -1))) {
    322401            flag[j] = (uChar)True;
    323           }
    324           else {
    325             spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
    326               * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
     402            continue;
     403          }
     404          if (flagrowOn[row] == 0) {
     405            if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
     406              spec[j] = 0.0 ;
     407              flag[j] = (uChar)True;
     408            }
     409            else {
     410              spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
     411                * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
     412            }
    327413          }
    328414        }
    … …  
    331417        // Chopper-Wheel calibration (Ulich & Haas 1976)
    332418        for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
    333           if ( (sphot[j]-spsky[j]) == 0.0 ) {
    334             spec[j] = 0.0 ;
     419          //if at least one of off/sky/hot data unavailable, just flag the channel.
     420          if (((idsOff[2*j] == -1)&&(idsOff[2*j+1] == -1))||
     421              ((idsSky[2*j] == -1)&&(idsSky[2*j+1] == -1))||
     422              ((idsHot[2*j] == -1)&&(idsHot[2*j+1] == -1))) {
    335423            flag[j] = (uChar)True;
    336           }
    337           else {
    338             spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
     424            continue;
     425          }
     426          if (flagrowOn[row] == 0) {
     427            if ( (sphot[j]-spsky[j]) == 0.0 ) {
     428              spec[j] = 0.0 ;
     429              flag[j] = (uChar)True;
     430            }
     431            else {
     432              spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
     433            }
    339434          }
    340435        }

  • trunk/src/STMath.cpp

    r2974 r2978  
    1212
    1313#include <sstream>
     14#include <iostream>
    1415
    1516#include <casa/iomanip.h>
    … …  
    42664267      // accumulate data
    42674268      s->flagsCol_.get( irow, flag ) ;
     4269      //if row-flagged, all channels set flagged
     4270      if (s->getFlagRow(irow)) {
     4271        for (uInt k = 0; k < nchan; ++k) {
     4272          flag(k) = 1 << 7;
     4273        }
     4274      }
    42684275      convertArray( bflag, flag ) ;
    42694276      s->specCol_.get( irow, spec ) ;
    42704277      tsys.assign( s->tsysCol_( irow ) ) ;
    4271       if ( !allEQ(bflag,True) )
    4272         acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
     4278      //if ( !allEQ(bflag,True) )
     4279      acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
    42734280      double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ;
    42744281      //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
    … …  
    43034310    irow = rows[len-1] ;
    43044311    s->flagsCol_.get( irow, flag ) ;
     4312    //if row-flagged, all channels set flagged
     4313    if (s->getFlagRow(irow)) {
     4314      for (uInt k = 0; k < nchan; ++k) {
     4315        flag(k) = 1 << 7;
     4316      }
     4317    }
    43054318    convertArray( bflag, flag ) ;
    43064319    s->specCol_.get( irow, spec ) ;
    43074320    tsys.assign( s->tsysCol_( irow ) ) ;
    4308     if (!allEQ(bflag,True) )
    4309       acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
     4321    //if (!allEQ(bflag,True) )
     4322    acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
    43104323    if ( acc.state() ) {
    43114324      atab.addRow() ;
Note: See TracChangeset for help on using the changeset viewer.