Changeset 1633 for branches


Ignore:
Timestamp:
09/16/09 17:11:16 (15 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-1422

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: Defined calibrate() and almacal() in python/asapmath.py

Defined cwcal() in STMath class

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): task_sdaverage

Description: Describe your changes here...

asapmath.py is changed to be able to calibrate data for APEX.
This modification is tentatively considered a calibration of ALMA
single-dish data. However, it must be updated in the future since
I don't know how raw ALMA data are provided and I have to change
code along with read ALMA data.

The calibrate() function takes calibration mode from its argument and
looks antenna name of the input scantable, and calls appropriate calibration
function depending on the calibration mode and antenna name.
If antenna name include 'APEX' or 'ALMA', newly defined calibration function
apexcal() is called. For other antenna name, one of the existing calibration
functions (calps, calnod, calfs, auto_quotient) is called.

Location:
branches/alma
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/python/asapmath.py

    r1631 r1633  
    850850    return s
    851851
    852 ## def apexcal( scantab, calmode ):
    853 ##     """
    854 ##     Calibrate APEX data.
     852def calibrate( scantab, scannos=[], calmode='none', verify=None ):
     853    """
     854    Calibrate data.
    855855   
    856 ##     Parameters:
    857 ##         scantab:       scantable
    858 ##         calmode:       calibration mode
    859 ##     """
    860 ##     from asap._asap import stmath
    861 ##     stm = stmath()
    862 ##     if ( calmode == 'ps' ):
    863 ##         asaplog.push( 'APEX position-switch calibration' )
    864 ##         print_log()
    865 ##     elif ( calmode == 'fs' ):
    866 ##         asaplog.push( 'APEX frequency-switch calibration' )
    867 ##         print_log()
    868 ##     elif ( calmode == 'wob' ):
    869 ##         asaplog.push( 'APEX wobbler-switch calibration' )
    870 ##         print_log()
    871 ##     elif ( calmode == 'otf' ):
    872 ##         asaplog.push( 'APEX On-The-Fly calibration' )
    873 ##         print_log()
    874 ##     else:
    875 ##         asaplog.push( '%s: unknown calibration mode' % calmode )
    876 ##         print_log('ERROR')
     856    Parameters:
     857        scantab:       scantable
     858        scannos:       list of scan number
     859        calmode:       calibration mode
     860        verify:        verify calibration     
     861    """
     862    print 'param: %s, %s, %s' %(scannos, calmode, verify)
     863    antname = scantab.get_antennaname()
     864    if ( calmode == 'nod' ):
     865        asaplog.push( 'Calibrating nod data.' )
     866        print_log()
     867        scal = calnod( scantab, scannos=scannos, verify=verify )
     868    elif ( calmode == 'quotient' ):
     869        asaplog.push( 'Calibrating using quotient.' )
     870        print_log()
     871        scal = scantab.auto_quotient( verify=verify )
     872    elif ( calmode == 'ps' ):
     873        asaplog.push( 'Calibrating %s position-switched data.' % antname )
     874        print_log()
     875        if ( antname.find( 'APEX' ) != -1 or antname.find( 'ALMA' ) != -1 ):
     876            scal = almacal( scantab, scannos, calmode, verify )
     877        else:
     878            scal = calps( scantab, scannos=scannos, verify=verify )
     879    elif ( calmode == 'fs' or calmode == 'fsotf' ):
     880        asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
     881        print_log()
     882        if ( antname.find( 'APEX' ) != -1 or antname.find( 'ALMA' ) != -1 ):
     883            scal = almacal( scantab, scannos, calmode, verify )
     884        else:
     885            scal = calfs( scantab, scannos=scannos, verify=verify )
     886    elif ( calmode == 'otf' ):
     887        asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
     888        print_log()
     889        scal = almacal( scantab, scannos, calmode, verify )
     890    else:
     891        asaplog.push( 'No calibration.' )
     892        scal = scantab.copy()
     893
     894    return scal
     895
     896def almacal( scantab, scannos=[], calmode='none', verify=False ):
     897    """
     898    Calibrate APEX and ALMA data
     899
     900    Parameters:
     901        scantab:       scantable
     902        scannos:       list of scan number
     903        calmode:       calibration mode
     904        verify:        verify calibration     
     905    """
     906    from asap._asap import stmath
     907    stm = stmath()
     908    antname = scantab.get_antennaname()
     909    ssub = scantab.get_scan( scannos )
     910    scal = scantable( stm.cwcal( ssub, calmode, antname ) )
     911    return scal

  • branches/alma/src/STMath.cpp

    r1616 r1633  
    11991199    out = merge(tabs);
    12001200  }
    1201   else { //folding is not implemented yet
     1201  else {
    12021202    //out = out1;
    1203     Int choffset = static_cast<Int>((rv1-rv2)/inc2) ;
     1203    Double choffset = ( rv1 - rv2 ) / inc2 ;
    12041204    out = dofold( out1, out2, choffset ) ;
    12051205  }
    … …  
    12101210CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
    12111211                                      const CountedPtr<Scantable> &ref,
    1212                                       Int choffset )
    1213 {
     1212                                      Double choffset,
     1213                                      Double choffset2 )
     1214{
     1215  LogIO os( LogOrigin( "STMath", "dofold", WHERE ) ) ;
     1216  os << "choffset=" << choffset << " choffset2=" << choffset2 << LogIO::POST ;
     1217
    12141218  // output scantable
    12151219  CountedPtr<Scantable> out = getScantable( sig, false ) ;
     1220
     1221  // separate choffset to integer part and decimal part
     1222  Int ioffset = (Int)choffset ;
     1223  Double doffset = choffset - ioffset ;
     1224  Int ioffset2 = (Int)choffset2 ;
     1225  Double doffset2 = choffset2 - ioffset2 ;
     1226  os << "ioffset=" << ioffset << " doffset=" << doffset << LogIO::POST ;
     1227  os << "ioffset2=" << ioffset2 << " doffset2=" << doffset2 << LogIO::POST ; 
    12161228
    12171229  // get column
    … …  
    12281240
    12291241  // check
    1230   if ( choffset == 0 ) {
     1242  if ( ioffset == 0 ) {
    12311243    LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
    12321244    os << "channel offset is zero, no folding" << LogIO::POST ;
    … …  
    12341246  }
    12351247  int nchan = ref->nchan() ;
    1236   if ( abs(choffset) >= nchan ) {
     1248  if ( abs(ioffset) >= nchan ) {
    12371249    LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
    12381250    os << "out-band frequency switching, no folding" << LogIO::POST ;
    … …  
    12461258  ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
    12471259  ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
     1260  ScalarColumn<uInt> fidColOut( out->table(), "FREQ_ID" ) ;
    12481261
    12491262  // for each row
    … …  
    12751288    // shift reference spectra
    12761289    int refchan = spref.nelements() ;
    1277     if ( choffset > 0 ) {
    1278       for ( int j = 0 ; j < refchan-choffset ; j++ ) {
    1279         spref[j] = spref[j+choffset] ;
    1280         tsysref[j] = tsysref[j+choffset] ;
    1281         flagref[j] = flagref[j+choffset] ;
    1282       }
    1283       for ( int j = refchan-choffset ; j < refchan ; j++ ) {
    1284         spref[j] = spref[j-refchan+choffset] ;
    1285         tsysref[j] = tsysref[j-refchan+choffset] ;
    1286         flagref[j] = flagref[j-refchan+choffset] ;
     1290    Vector<Float> sspref( spref.nelements() ) ;
     1291    Vector<Float> stsysref( tsysref.nelements() ) ;
     1292    Vector<uChar> sflagref( flagref.nelements() ) ;
     1293    if ( ioffset > 0 ) {
     1294      // SPECTRA and FLAGTRA
     1295      for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
     1296        sspref[j] = spref[j+ioffset] ;
     1297        sflagref[j] = flagref[j+ioffset] ;
     1298      }
     1299      for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
     1300        sspref[j] = spref[j-refchan+ioffset] ;
     1301        sflagref[j] = flagref[j-refchan+ioffset] ;
     1302      }
     1303      spref = sspref.copy() ;
     1304      flagref = sflagref.copy() ;
     1305      for ( int j = 0 ; j < refchan - 1 ; j++ ) {
     1306        sspref[j] = doffset * spref[j+1] + ( 1.0 - doffset ) * spref[j] ;
     1307        sflagref[j] = flagref[j+1] + flagref[j] ;
     1308      }
     1309      sspref[refchan-1] = doffset * spref[0] + ( 1.0 - doffset ) * spref[refchan-1] ;
     1310      sflagref[refchan-1] = flagref[0] + flagref[refchan-1] ;
     1311
     1312      // TSYS
     1313      if ( spref.nelements() == tsysref.nelements() ) {
     1314        for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
     1315          stsysref[j] = tsysref[j+ioffset] ;
     1316        }
     1317        for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
     1318          stsysref[j] = tsysref[j-refchan+ioffset] ;
     1319        }
     1320        tsysref = stsysref.copy() ;
     1321        for ( int j = 0 ; j < refchan - 1 ; j++ ) {
     1322          stsysref[j] = doffset * tsysref[j+1] + ( 1.0 - doffset ) * tsysref[j] ;
     1323        }
     1324        stsysref[refchan-1] = doffset * tsysref[0] + ( 1.0 - doffset ) * tsysref[refchan-1] ;
    12871325      }
    12881326    }
    12891327    else {
    1290       for ( int j = 0 ; j < abs(choffset) ; j++ ) {
    1291         spref[j] = spref[refchan+choffset+j] ;
    1292         tsysref[j] = tsysref[refchan+choffset+j] ;
    1293         flagref[j] = flagref[refchan+choffset+j] ;
    1294       }
    1295       for ( int j = abs(choffset) ; j < refchan ; j++ ) {
    1296         spref[j] = spref[j+choffset] ;
    1297         tsysref[j] = tsysref[j+choffset] ;
    1298         flagref[j] = flagref[j+choffset] ;
     1328      // SPECTRA and FLAGTRA
     1329      for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
     1330        sspref[j] = spref[refchan+ioffset+j] ;
     1331        sflagref[j] = flagref[refchan+ioffset+j] ;
     1332      }
     1333      for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
     1334        sspref[j] = spref[j+ioffset] ;
     1335        sflagref[j] = flagref[j+ioffset] ;
     1336      }
     1337      spref = sspref.copy() ;
     1338      flagref = sflagref.copy() ;
     1339      sspref[0] = doffset * spref[refchan-1] + ( 1.0 - doffset ) * spref[0] ;
     1340      sflagref[0] = flagref[0] + flagref[refchan-1] ;
     1341      for ( int j = 1 ; j < refchan ; j++ ) {
     1342        sspref[j] = doffset * spref[j-1] + ( 1.0 - doffset ) * spref[j] ;
     1343        sflagref[j] = flagref[j-1] + flagref[j] ;
     1344      }
     1345      // TSYS
     1346      if ( spref.nelements() == tsysref.nelements() ) {
     1347        for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
     1348          stsysref[j] = tsysref[refchan+ioffset+j] ;
     1349        }
     1350        for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
     1351          stsysref[j] = tsysref[j+ioffset] ;
     1352        }
     1353        tsysref = stsysref.copy() ;
     1354        stsysref[0] = doffset * tsysref[refchan-1] + ( 1.0 - doffset ) * tsysref[0] ;
     1355        for ( int j = 1 ; j < refchan ; j++ ) {
     1356          stsysref[j] = doffset * tsysref[j-1] + ( 1.0 - doffset ) * tsysref[j] ;
     1357        }
     1358      }
     1359    }
     1360
     1361    // shift signal spectra if necessary (only for APEX?)
     1362    if ( choffset2 != 0.0 ) {
     1363      int sigchan = spsig.nelements() ;
     1364      Vector<Float> sspsig( spsig.nelements() ) ;
     1365      Vector<Float> stsyssig( tsyssig.nelements() ) ;
     1366      Vector<uChar> sflagsig( flagsig.nelements() ) ;
     1367      if ( ioffset2 > 0 ) {
     1368        // SPECTRA and FLAGTRA
     1369        for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
     1370          sspsig[j] = spsig[j+ioffset2] ;
     1371          sflagsig[j] = flagsig[j+ioffset2] ;
     1372        }
     1373        for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
     1374          sspsig[j] = spsig[j-sigchan+ioffset2] ;
     1375          sflagsig[j] = flagsig[j-sigchan+ioffset2] ;
     1376        }
     1377        spsig = sspsig.copy() ;
     1378        flagsig = sflagsig.copy() ;
     1379        for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
     1380          sspsig[j] = doffset2 * spsig[j+1] + ( 1.0 - doffset2 ) * spsig[j] ;
     1381          sflagsig[j] = flagsig[j+1] || flagsig[j] ;
     1382        }
     1383        sspsig[sigchan-1] = doffset2 * spsig[0] + ( 1.0 - doffset2 ) * spsig[sigchan-1] ;
     1384        sflagsig[sigchan-1] = flagsig[0] || flagsig[sigchan-1] ;
     1385        // TSTS
     1386        if ( spsig.nelements() == tsyssig.nelements() ) {
     1387          for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
     1388            stsyssig[j] = tsyssig[j+ioffset2] ;
     1389          }
     1390          for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
     1391            stsyssig[j] = tsyssig[j-sigchan+ioffset2] ;
     1392          }
     1393          tsyssig = stsyssig.copy() ;
     1394          for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
     1395            stsyssig[j] = doffset2 * tsyssig[j+1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
     1396          }
     1397          stsyssig[sigchan-1] = doffset2 * tsyssig[0] + ( 1.0 - doffset2 ) * tsyssig[sigchan-1] ;
     1398        }
     1399      }
     1400      else {
     1401        // SPECTRA and FLAGTRA
     1402        for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
     1403          sspsig[j] = spsig[sigchan+ioffset2+j] ;
     1404          sflagsig[j] = flagsig[sigchan+ioffset2+j] ;
     1405        }
     1406        for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
     1407          sspsig[j] = spsig[j+ioffset2] ;
     1408          sflagsig[j] = flagsig[j+ioffset2] ;
     1409        }
     1410        spsig = sspsig.copy() ;
     1411        flagsig = sflagsig.copy() ;
     1412        sspsig[0] = doffset2 * spsig[sigchan-1] + ( 1.0 - doffset2 ) * spsig[0] ;
     1413        sflagsig[0] = flagsig[0] + flagsig[sigchan-1] ;
     1414        for ( int j = 1 ; j < sigchan ; j++ ) {
     1415          sspsig[j] = doffset2 * spsig[j-1] + ( 1.0 - doffset2 ) * spsig[j] ;
     1416          sflagsig[j] = flagsig[j-1] + flagsig[j] ;
     1417        }
     1418        // TSYS
     1419        if ( spsig.nelements() == tsyssig.nelements() ) {
     1420          for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
     1421            stsyssig[j] = tsyssig[sigchan+ioffset2+j] ;
     1422          }
     1423          for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
     1424            stsyssig[j] = tsyssig[j+ioffset2] ;
     1425          }
     1426          tsyssig = stsyssig.copy() ;
     1427          stsyssig[0] = doffset2 * tsyssig[sigchan-1] + ( 1.0 - doffset2 ) * tsyssig[0] ;
     1428          for ( int j = 1 ; j < sigchan ; j++ ) {
     1429            stsyssig[j] = doffset2 * tsyssig[j-1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
     1430          }
     1431        }
    12991432      }
    13001433    }
    … …  
    13021435    // folding
    13031436    acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
    1304     acc.add( spref, !flagref, tsysref, intref, timeref ) ;
     1437    acc.add( sspref, !sflagref, stsysref, intref, timeref ) ;
    13051438   
    13061439    // put result
    … …  
    13131446    intervalColOut.put( i, acc.getInterval() ) ;
    13141447    mjdColOut.put( i, acc.getTime() ) ;
     1448    // change FREQ_ID to unshifted IF setting (only for APEX?)
     1449    if ( choffset2 != 0.0 ) {
     1450      int freqid = fidColOut( 0 ) ; // assume single-IF data
     1451      double refpix, refval, increment ;
     1452      out->frequencies().getEntry( refpix, refval, increment, freqid ) ;
     1453      refval -= choffset * increment ;
     1454      uInt newfreqid = out->frequencies().addEntry( refpix, refval, increment ) ;
     1455      Vector<uInt> freqids = fidColOut.getColumn() ;
     1456      for ( uInt j = 0 ; j < freqids.nelements() ; j++ ) {
     1457        if ( freqids[j] == freqid )
     1458          freqids[j] = newfreqid ;
     1459      }
     1460      fidColOut.putColumn( freqids ) ;
     1461    }
    13151462
    13161463    acc.reset() ;
    … …  
    31493296  return out ;
    31503297}
     3298
     3299CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
     3300                                     const String calmode,
     3301                                     const String antname )
     3302{
     3303  // frequency switch
     3304  if ( calmode == "fs" ) {
     3305    return cwcalfs( s, antname ) ;
     3306  }
     3307  else {
     3308    string skystr = "*_sky" ;
     3309    string hotstr = "*_hot" ;
     3310    string coldstr = "*_cold" ;
     3311    string onstr = "*_" ;
     3312    string offstr = "*_" ;
     3313   
     3314    if ( calmode == "ps" || calmode == "otf" ) {
     3315      onstr += "pson" ;
     3316      offstr += "psoff" ;
     3317    }
     3318    else if ( calmode == "wob" ) {
     3319      onstr += "wobon" ;
     3320      offstr += "woboff" ;
     3321    }
     3322   
     3323    vector<bool> masks = s->getMask( 0 ) ;
     3324   
     3325    // sky scan
     3326    STSelector sel = STSelector() ;
     3327    sel.setName( skystr ) ;
     3328    s->setSelection( sel ) ;
     3329    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
     3330    CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
     3331    s->unsetSelection() ;
     3332    sel.reset() ;
     3333
     3334    // hot scan
     3335    sel.setName( hotstr ) ;
     3336    s->setSelection( sel ) ;
     3337    tmp[0] = getScantable( s, false ) ;
     3338    CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
     3339    s->unsetSelection() ;
     3340    sel.reset() ;
     3341   
     3342    // cold scan
     3343    sel.setName( coldstr ) ;
     3344    s->setSelection( sel ) ;
     3345    tmp[0] = getScantable( s, false ) ;
     3346    CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCNAN" ) ;
     3347    s->unsetSelection() ;
     3348    sel.reset() ;
     3349
     3350    // off scan
     3351    sel.setName( offstr ) ;
     3352    s->setSelection( sel ) ;
     3353    tmp[0] = getScantable( s, false ) ;
     3354    CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
     3355    s->unsetSelection() ;
     3356    sel.reset() ;
     3357   
     3358    // on scan
     3359    bool insitu = insitu_ ;
     3360    insitu_ = false ;
     3361    CountedPtr<Scantable> out = getScantable( s, true ) ;
     3362    insitu_ = insitu ;
     3363    sel.setName( onstr ) ;
     3364    s->setSelection( sel ) ;
     3365    TableCopy::copyRows( out->table(), s->table() ) ;
     3366    s->unsetSelection() ;
     3367    sel.reset() ;
     3368   
     3369    // process each on scan
     3370    for ( int i = 0 ; i < out->nrow() ; i++ ) {
     3371      vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ;
     3372      out->setSpectrum( sp, i ) ;
     3373    }
     3374
     3375    // remove additional string from SRCNAME
     3376    ScalarColumn<String> srcnameCol ;
     3377    srcnameCol.attach( out->table(), "SRCNAME" ) ;
     3378    Vector<String> srcnames( srcnameCol.getColumn() ) ;
     3379    for ( uInt i = 0 ; i < srcnames.nelements() ; i++ ) {
     3380      srcnames[i] = srcnames[i].substr( 0, srcnames[i].find( onstr.substr(1,onstr.size()-1) ) ) ;
     3381    }
     3382    srcnameCol.putColumn( srcnames ) ;
     3383
     3384    return out ;
     3385  }
     3386}
     3387 
     3388CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
     3389                                       const String antname )
     3390{
     3391  string skystr = "*_sky" ;
     3392  string skystr1 = "*_fslo_sky" ;
     3393  string skystr2 = "*_fshi_sky" ;
     3394  string hotstr = "*_hot" ;
     3395  string hotstr1 = "*_fslo_hot" ;
     3396  string hotstr2 = "*_fshi_hot" ;
     3397  string coldstr = "*_cold" ;
     3398  string coldstr1 = "*_fslo_cold" ;
     3399  string coldstr2 = "*_fshi_cold" ;
     3400  string offstr1 = "*_fslo_off" ;
     3401  string offstr2 = "*_fshi_off" ;
     3402  string sigstr = "*_" ;
     3403  string refstr = "*_" ;
     3404
     3405  // APEX calibration mode
     3406  int apexcalmode = 1 ;
     3407 
     3408  if ( antname.find( "APEX" ) != string::npos ) {
     3409    sigstr += "fslo" ;
     3410    refstr += "fshi" ;   
     3411
     3412    // check if off scan exists or not
     3413    STSelector sel = STSelector() ;
     3414    sel.setName( offstr1 ) ;
     3415    try {
     3416      s->setSelection( sel ) ;
     3417    }
     3418    catch ( AipsError &e ) {
     3419      apexcalmode = 0 ;
     3420    }
     3421  }
     3422  else {
     3423    sigstr += "fssig" ;
     3424    refstr += "fsref" ;
     3425  }
     3426
     3427  vector<bool> masks = s->getMask( 0 ) ;
     3428  CountedPtr<Scantable> ssig, sref ;
     3429  CountedPtr<Scantable> out ;
     3430
     3431  if ( antname.find( "APEX" ) != string::npos && apexcalmode == 0 ) {
     3432    // APEX fs data without off scan
     3433    // sky scan
     3434    STSelector sel = STSelector() ;
     3435    sel.setName( skystr1 ) ;
     3436    s->setSelection( sel ) ;
     3437    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
     3438    CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ;
     3439    s->unsetSelection() ;
     3440    sel.reset() ;
     3441    sel.setName( skystr2 ) ;
     3442    s->setSelection( sel ) ;
     3443    tmp[0] = getScantable( s, false ) ;
     3444    CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ;
     3445    s->unsetSelection() ;
     3446    sel.reset() ;
     3447   
     3448    // hot scan
     3449    sel.setName( hotstr1 ) ;
     3450    s->setSelection( sel ) ;
     3451    tmp[0] = getScantable( s, false ) ;
     3452    CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ;
     3453    s->unsetSelection() ;
     3454    sel.reset() ;
     3455    sel.setName( hotstr2 ) ;
     3456    s->setSelection( sel ) ;
     3457    tmp[0] = getScantable( s, false ) ;
     3458    CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ;
     3459    s->unsetSelection() ;
     3460    sel.reset() ;
     3461   
     3462    // cold scan
     3463    sel.setName( coldstr1 ) ;
     3464    s->setSelection( sel ) ;
     3465    tmp[0] = getScantable( s, false ) ;
     3466    CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ;
     3467    s->unsetSelection() ;
     3468    sel.reset() ;
     3469    sel.setName( coldstr2 ) ;
     3470    s->setSelection( sel ) ;
     3471    tmp[0] = getScantable( s, false ) ;
     3472    CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ;
     3473    s->unsetSelection() ;
     3474    sel.reset() ;
     3475   
     3476    // ref scan
     3477    bool insitu = insitu_ ;
     3478    insitu_ = false ;
     3479    sref = getScantable( s, true ) ;
     3480    insitu_ = insitu ;
     3481    sel.setName( refstr ) ;
     3482    s->setSelection( sel ) ;
     3483    TableCopy::copyRows( sref->table(), s->table() ) ;
     3484    s->unsetSelection() ;
     3485    sel.reset() ;
     3486   
     3487    // sig scan
     3488    insitu_ = false ;
     3489    ssig = getScantable( s, true ) ;
     3490    insitu_ = insitu ;
     3491    sel.setName( sigstr ) ;
     3492    s->setSelection( sel ) ;
     3493    TableCopy::copyRows( ssig->table(), s->table() ) ;
     3494    s->unsetSelection() ;
     3495    sel.reset() ; 
     3496
     3497    // process each sig and ref scan
     3498    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
     3499      vector< CountedPtr<Scantable> > sky( 2 ) ;
     3500      sky[0] = askylo ;
     3501      sky[1] = askyhi ;
     3502      vector< CountedPtr<Scantable> > hot( 2 ) ;
     3503      hot[0] = ahotlo ;
     3504      hot[1] = ahothi ;
     3505      vector< CountedPtr<Scantable> > cold( 2 ) ;
     3506      cold[0] = acoldlo ;
     3507      cold[1] = acoldhi ;
     3508      vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ;
     3509      ssig->setSpectrum( sp, i ) ;
     3510      sky[0] = askyhi ;
     3511      sky[1] = askylo ;
     3512      hot[0] = ahothi ;
     3513      hot[1] = ahotlo ;
     3514      cold[0] = acoldhi ;
     3515      cold[1] = acoldlo ;
     3516      sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ;
     3517      sref->setSpectrum( sp, i ) ;
     3518    }
     3519
     3520  }
     3521  else if ( antname.find( "APEX" ) != string::npos && apexcalmode == 1 ) {
     3522    // APEX fs data with off scan
     3523    // sky scan
     3524    STSelector sel = STSelector() ;
     3525    sel.setName( skystr1 ) ;
     3526    s->setSelection( sel ) ;
     3527    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
     3528    CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ;
     3529    s->unsetSelection() ;
     3530    sel.reset() ;
     3531    sel.setName( skystr2 ) ;
     3532    s->setSelection( sel ) ;
     3533    tmp[0] = getScantable( s, false ) ;
     3534    CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ;
     3535    s->unsetSelection() ;
     3536    sel.reset() ;
     3537   
     3538    // hot scan
     3539    sel.setName( hotstr1 ) ;
     3540    s->setSelection( sel ) ;
     3541    tmp[0] = getScantable( s, false ) ;
     3542    CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ;
     3543    s->unsetSelection() ;
     3544    sel.reset() ;
     3545    sel.setName( hotstr2 ) ;
     3546    s->setSelection( sel ) ;
     3547    tmp[0] = getScantable( s, false ) ;
     3548    CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ;
     3549    s->unsetSelection() ;
     3550    sel.reset() ;
     3551   
     3552    // cold scan
     3553    sel.setName( coldstr1 ) ;
     3554    s->setSelection( sel ) ;
     3555    tmp[0] = getScantable( s, false ) ;
     3556    CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ;
     3557    s->unsetSelection() ;
     3558    sel.reset() ;
     3559    sel.setName( coldstr2 ) ;
     3560    s->setSelection( sel ) ;
     3561    tmp[0] = getScantable( s, false ) ;
     3562    CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ;
     3563    s->unsetSelection() ;
     3564    sel.reset() ;
     3565
     3566    // off scan
     3567    sel.setName( offstr1 ) ;
     3568    s->setSelection( sel ) ;
     3569    tmp[0] = getScantable( s, false ) ;
     3570    CountedPtr<Scantable> aofflo = average( tmp, masks, "TINT", "SCAN" ) ;
     3571    s->unsetSelection() ;
     3572    sel.reset() ;
     3573    sel.setName( offstr2 ) ;
     3574    s->setSelection( sel ) ;
     3575    tmp[0] = getScantable( s, false ) ;
     3576    CountedPtr<Scantable> aoffhi = average( tmp, masks, "TINT", "SCAN" ) ;
     3577    s->unsetSelection() ;
     3578    sel.reset() ;
     3579   
     3580    // ref scan
     3581    bool insitu = insitu_ ;
     3582    insitu_ = false ;
     3583    sref = getScantable( s, true ) ;
     3584    insitu_ = insitu ;
     3585    sel.setName( refstr ) ;
     3586    s->setSelection( sel ) ;
     3587    TableCopy::copyRows( sref->table(), s->table() ) ;
     3588    s->unsetSelection() ;
     3589    sel.reset() ;
     3590   
     3591    // sig scan
     3592    insitu_ = false ;
     3593    ssig = getScantable( s, true ) ;
     3594    insitu_ = insitu ;
     3595    sel.setName( sigstr ) ;
     3596    s->setSelection( sel ) ;
     3597    TableCopy::copyRows( ssig->table(), s->table() ) ;
     3598    s->unsetSelection() ;
     3599    sel.reset() ;
     3600
     3601    // process each sig and ref scan
     3602    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
     3603      vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ;
     3604      ssig->setSpectrum( sp, i ) ;
     3605      sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ;
     3606      sref->setSpectrum( sp, i ) ;
     3607    }
     3608  }
     3609  else if ( antname.find( "APEX" ) == string::npos ) {
     3610    // non-APEX fs data
     3611    // sky scan
     3612    STSelector sel = STSelector() ;
     3613    sel.setName( skystr ) ;
     3614    s->setSelection( sel ) ;
     3615    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
     3616    CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
     3617    s->unsetSelection() ;
     3618    sel.reset() ;
     3619   
     3620    // hot scan
     3621    sel.setName( hotstr ) ;
     3622    s->setSelection( sel ) ;
     3623    tmp[0] = getScantable( s, false ) ;
     3624    CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
     3625    s->unsetSelection() ;
     3626    sel.reset() ;
     3627
     3628     // cold scan
     3629    sel.setName( coldstr ) ;
     3630    s->setSelection( sel ) ;
     3631    tmp[0] = getScantable( s, false ) ;
     3632    CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCAN" ) ;
     3633    s->unsetSelection() ;
     3634    sel.reset() ;
     3635   
     3636    // ref scan
     3637    bool insitu = insitu_ ;
     3638    insitu_ = false ;
     3639    sref = getScantable( s, true ) ;
     3640    insitu_ = insitu ;
     3641    sel.setName( refstr ) ;
     3642    s->setSelection( sel ) ;
     3643    TableCopy::copyRows( sref->table(), s->table() ) ;
     3644    s->unsetSelection() ;
     3645    sel.reset() ;
     3646   
     3647    // sig scan
     3648    insitu_ = false ;
     3649    ssig = getScantable( s, true ) ;
     3650    insitu_ = insitu ;
     3651    sel.setName( sigstr ) ;
     3652    s->setSelection( sel ) ;
     3653    TableCopy::copyRows( ssig->table(), s->table() ) ;
     3654    s->unsetSelection() ;
     3655    sel.reset() ;
     3656
     3657    // process each sig and ref scan
     3658    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
     3659      vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ;
     3660      ssig->setSpectrum( sp, i ) ;
     3661      sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ;
     3662      sref->setSpectrum( sp, i ) ;
     3663    }
     3664  }
     3665
     3666  // do folding if necessary
     3667  Table sigtab = ssig->table() ;
     3668  Table reftab = sref->table() ;
     3669  ScalarColumn<uInt> sigifnoCol ;
     3670  ScalarColumn<uInt> refifnoCol ;
     3671  ScalarColumn<uInt> sigfidCol ;
     3672  ScalarColumn<uInt> reffidCol ;
     3673  uInt nchan = (uInt)ssig->nchan() ;
     3674  sigifnoCol.attach( sigtab, "IFNO" ) ;
     3675  refifnoCol.attach( reftab, "IFNO" ) ;
     3676  sigfidCol.attach( sigtab, "FREQ_ID" ) ;
     3677  reffidCol.attach( reftab, "FREQ_ID" ) ;
     3678  Vector<uInt> sfids( sigfidCol.getColumn() ) ;
     3679  Vector<uInt> rfids( reffidCol.getColumn() ) ;
     3680  vector<uInt> sfids_unique ;
     3681  vector<uInt> rfids_unique ;
     3682  vector<uInt> sifno_unique ;
     3683  vector<uInt> rifno_unique ;
     3684  for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
     3685    if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
     3686      sfids_unique.push_back( sfids[i] ) ;
     3687      sifno_unique.push_back( ssig->getIF( i ) ) ;
     3688    }
     3689    if ( count( rfids_unique.begin(), rfids_unique.end(),  rfids[i] ) == 0 ) {
     3690      rfids_unique.push_back( rfids[i] ) ;
     3691      rifno_unique.push_back( sref->getIF( i ) ) ;
     3692    }
     3693  }
     3694  double refpix_sig, refval_sig, increment_sig ;
     3695  double refpix_ref, refval_ref, increment_ref ;
     3696  vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
     3697  for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
     3698    ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
     3699    sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
     3700    if ( refpix_sig == refpix_ref ) {
     3701      double foffset = refval_ref - refval_sig ;
     3702      int choffset = static_cast<int>(foffset/increment_sig) ;
     3703      double doffset = foffset / increment_sig ;
     3704      if ( abs(choffset) >= nchan ) {
     3705        LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
     3706        os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
     3707        os << "Just return signal data" << LogIO::POST ;
     3708        //std::vector< CountedPtr<Scantable> > tabs ;
     3709        //tabs.push_back( ssig ) ;
     3710        //tabs.push_back( sref ) ;
     3711        //out = merge( tabs ) ;
     3712        tmp[i] = ssig ;
     3713      }
     3714      else {
     3715        STSelector sel = STSelector() ;
     3716        vector<int> v( 1, sifno_unique[i] ) ;
     3717        sel.setIFs( v ) ;
     3718        ssig->setSelection( sel ) ;
     3719        sel.reset() ;
     3720        v[0] = rifno_unique[i] ;
     3721        sel.setIFs( v ) ;
     3722        sref->setSelection( sel ) ;
     3723        sel.reset() ;
     3724        if ( antname.find( "APEX" ) != string::npos ) {
     3725          tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ;
     3726          //tmp[i] = dofold( ssig, sref, doffset ) ;
     3727        }
     3728        else {
     3729          tmp[i] = dofold( ssig, sref, doffset ) ;
     3730        }
     3731        // remove additional string from SRCNAME
     3732        ScalarColumn<String> srcnameCol ;
     3733        srcnameCol.attach( tmp[i]->table(), "SRCNAME" ) ;
     3734        Vector<String> srcnames( srcnameCol.getColumn() ) ;
     3735        for ( uInt i = 0 ; i < srcnames.nelements() ; i++ ) {
     3736          srcnames[i] = srcnames[i].substr( 0, srcnames[i].find( sigstr.substr(1,sigstr.size()-1) ) ) ;
     3737        }
     3738        srcnameCol.putColumn( srcnames ) ;
     3739        ssig->unsetSelection() ;
     3740        sref->unsetSelection() ;
     3741      }
     3742    }
     3743  }
     3744
     3745  out = merge( tmp ) ;
     3746
     3747  return out ;
     3748}
     3749
     3750vector<float> STMath::getSpectrumFromTime( string reftime,
     3751                                           CountedPtr<Scantable>& s,
     3752                                           string mode )
     3753{
     3754  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
     3755  vector<float> sp ;
     3756
     3757  if ( s->nrow() == 0 ) {
     3758    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
     3759    return sp ;
     3760  }
     3761  else if ( s->nrow() == 1 ) {
     3762    os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
     3763    return s->getSpectrum( 0 ) ;
     3764  }
     3765  else {
     3766    vector<int> idx = getRowIdFromTime( reftime, s ) ;
     3767    if ( mode == "before" ) {
     3768      int id = -1 ;
     3769      if ( idx[0] != -1 ) {
     3770        id = idx[0] ;
     3771      }
     3772      else if ( idx[1] != -1 ) {
     3773        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
     3774        id = idx[1] ;
     3775      }
     3776      os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     3777      sp = s->getSpectrum( id ) ;
     3778    }
     3779    else if ( mode == "after" ) {
     3780      int id = -1 ;
     3781      if ( idx[1] != -1 ) {
     3782        id = idx[1] ;
     3783      }
     3784      else if ( idx[0] != -1 ) {
     3785        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
     3786        id = idx[1] ;
     3787      }
     3788      os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     3789      sp = s->getSpectrum( id ) ;
     3790    }
     3791    else if ( mode == "nearest" ) {
     3792      int id = -1 ;
     3793      if ( idx[0] == -1 ) {
     3794        id = idx[1] ;
     3795      }
     3796      else if ( idx[1] == -1 ) {
     3797        id = idx[0] ;
     3798      }
     3799      else if ( idx[0] == idx[1] ) {
     3800        id = idx[0] ;
     3801      }
     3802      else {
     3803        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     3804        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     3805        double tref = getMJD( reftime ) ;
     3806        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     3807          id = idx[1] ;
     3808        }
     3809        else {
     3810          id = idx[0] ;
     3811        }
     3812      }
     3813      os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     3814      sp = s->getSpectrum( id ) ;     
     3815    }
     3816    else if ( mode == "linear" ) {
     3817      if ( idx[0] == -1 ) {
     3818        // use after
     3819        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
     3820        int id = idx[1] ;
     3821        os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     3822        sp = s->getSpectrum( id ) ;
     3823      }
     3824      else if ( idx[1] == -1 ) {
     3825        // use before
     3826        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
     3827        int id = idx[0] ;
     3828        os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     3829        sp = s->getSpectrum( id ) ;
     3830      }
     3831      else if ( idx[0] == idx[1] ) {
     3832        // use before
     3833        os << "No need to interporate." << LogIO::POST ;
     3834        int id = idx[0] ;
     3835        os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     3836        sp = s->getSpectrum( id ) ;
     3837      }
     3838      else {
     3839        // do interpolation
     3840        os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     3841        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     3842        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     3843        double tref = getMJD( reftime ) ;
     3844        vector<float> sp0 = s->getSpectrum( idx[0] ) ;
     3845        vector<float> sp1 = s->getSpectrum( idx[1] ) ;
     3846        for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
     3847          float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
     3848          sp.push_back( v ) ;
     3849        }
     3850      }
     3851    }
     3852    else {
     3853      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     3854    }
     3855    return sp ;
     3856  }
     3857}
     3858
     3859double STMath::getMJD( string strtime )
     3860{
     3861  if ( strtime.find("/") == string::npos ) {
     3862    // MJD time string
     3863    return atof( strtime.c_str() ) ;
     3864  }
     3865  else {
     3866    // string in YYYY/MM/DD/HH:MM:SS format
     3867    uInt year = atoi( strtime.substr( 0, 4 ).c_str() ) ;
     3868    uInt month = atoi( strtime.substr( 5, 2 ).c_str() ) ;
     3869    uInt day = atoi( strtime.substr( 8, 2 ).c_str() ) ;
     3870    uInt hour = atoi( strtime.substr( 11, 2 ).c_str() ) ;
     3871    uInt minute = atoi( strtime.substr( 14, 2 ).c_str() ) ;
     3872    uInt sec = atoi( strtime.substr( 17, 2 ).c_str() ) ;
     3873    Time t( year, month, day, hour, minute, sec ) ;
     3874    return t.modifiedJulianDay() ;
     3875  }
     3876}
     3877
     3878vector<int> STMath::getRowIdFromTime( string reftime, CountedPtr<Scantable> &s )
     3879{
     3880  double reft = getMJD( reftime ) ;
     3881  double dtmin = 1.0e100 ;
     3882  double dtmax = -1.0e100 ;
     3883  vector<double> dt ;
     3884  int just_before = -1 ;
     3885  int just_after = -1 ;
     3886  for ( int i = 0 ; i < s->nrow() ; i++ ) {
     3887    dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;
     3888  }
     3889  for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
     3890    if ( dt[i] > 0.0 ) {
     3891      // after reftime
     3892      if ( dt[i] < dtmin ) {
     3893        just_after = i ;
     3894        dtmin = dt[i] ;
     3895      }
     3896    }
     3897    else if ( dt[i] < 0.0 ) {
     3898      // before reftime
     3899      if ( dt[i] > dtmax ) {
     3900        just_before = i ;
     3901        dtmax = dt[i] ;
     3902      }
     3903    }
     3904    else {
     3905      // just a reftime
     3906      just_before = i ;
     3907      just_after = i ;
     3908      dtmax = 0 ;
     3909      dtmin = 0 ;
     3910      break ;
     3911    }
     3912  }
     3913
     3914  vector<int> v ;
     3915  v.push_back( just_before ) ;
     3916  v.push_back( just_after ) ;
     3917
     3918  return v ;
     3919}
     3920
     3921vector<float> STMath::getTcalFromTime( string reftime,
     3922                                       CountedPtr<Scantable>& s,
     3923                                       string mode )
     3924{
     3925  LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ;
     3926  vector<float> tcal ;
     3927  STTcal tcalTable = s->tcal() ;
     3928  String time ;
     3929  Vector<Float> tcalval ;
     3930  if ( s->nrow() == 0 ) {
     3931    os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;
     3932    return tcal ;
     3933  }
     3934  else if ( s->nrow() == 1 ) {
     3935    uInt tcalid = s->getTcalId( 0 ) ;
     3936    os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     3937    tcalTable.getEntry( time, tcalval, tcalid ) ;
     3938    tcalval.tovector( tcal ) ;
     3939    return tcal ;
     3940  }
     3941  else {
     3942    vector<int> idx = getRowIdFromTime( reftime, s ) ;
     3943    if ( mode == "before" ) {
     3944      int id = -1 ;
     3945      if ( idx[0] != -1 ) {
     3946        id = idx[0] ;
     3947      }
     3948      else if ( idx[1] != -1 ) {
     3949        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
     3950        id = idx[1] ;
     3951      }
     3952      uInt tcalid = s->getTcalId( id ) ;
     3953      os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     3954      tcalTable.getEntry( time, tcalval, tcalid ) ;
     3955      tcalval.tovector( tcal ) ;
     3956    }
     3957    else if ( mode == "after" ) {
     3958      int id = -1 ;
     3959      if ( idx[1] != -1 ) {
     3960        id = idx[1] ;
     3961      }
     3962      else if ( idx[0] != -1 ) {
     3963        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
     3964        id = idx[1] ;
     3965      }
     3966      uInt tcalid = s->getTcalId( id ) ;
     3967      os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     3968      tcalTable.getEntry( time, tcalval, tcalid ) ;
     3969      tcalval.tovector( tcal ) ;
     3970    }
     3971    else if ( mode == "nearest" ) {
     3972      int id = -1 ;
     3973      if ( idx[0] == -1 ) {
     3974        id = idx[1] ;
     3975      }
     3976      else if ( idx[1] == -1 ) {
     3977        id = idx[0] ;
     3978      }
     3979      else if ( idx[0] == idx[1] ) {
     3980        id = idx[0] ;
     3981      }
     3982      else {
     3983        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     3984        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     3985        double tref = getMJD( reftime ) ;
     3986        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     3987          id = idx[1] ;
     3988        }
     3989        else {
     3990          id = idx[0] ;
     3991        }
     3992      }
     3993      uInt tcalid = s->getTcalId( id ) ;
     3994      os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     3995      tcalTable.getEntry( time, tcalval, tcalid ) ;
     3996      tcalval.tovector( tcal ) ;
     3997    }
     3998    else if ( mode == "linear" ) {
     3999      if ( idx[0] == -1 ) {
     4000        // use after
     4001        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
     4002        int id = idx[1] ;
     4003        uInt tcalid = s->getTcalId( id ) ;
     4004        os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4005        tcalTable.getEntry( time, tcalval, tcalid ) ;
     4006        tcalval.tovector( tcal ) ;
     4007      }
     4008      else if ( idx[1] == -1 ) {
     4009        // use before
     4010        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
     4011        int id = idx[0] ;
     4012        uInt tcalid = s->getTcalId( id ) ;
     4013        os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4014        tcalTable.getEntry( time, tcalval, tcalid ) ;
     4015        tcalval.tovector( tcal ) ;
     4016      }
     4017      else if ( idx[0] == idx[1] ) {
     4018        // use before
     4019        os << "No need to interporate." << LogIO::POST ;
     4020        int id = idx[0] ;
     4021        uInt tcalid = s->getTcalId( id ) ;
     4022        os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4023        tcalTable.getEntry( time, tcalval, tcalid ) ;
     4024        tcalval.tovector( tcal ) ;
     4025      }
     4026      else {
     4027        // do interpolation
     4028        os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4029        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4030        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4031        double tref = getMJD( reftime ) ;
     4032        vector<float> tcal0 ;
     4033        vector<float> tcal1 ;
     4034        uInt tcalid0 = s->getTcalId( idx[0] ) ;
     4035        uInt tcalid1 = s->getTcalId( idx[1] ) ;
     4036        tcalTable.getEntry( time, tcalval, tcalid0 ) ;
     4037        tcalval.tovector( tcal0 ) ;
     4038        tcalTable.getEntry( time, tcalval, tcalid1 ) ;
     4039        tcalval.tovector( tcal1 ) ;       
     4040        for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) {
     4041          float v = ( tcal1[i] - tcal0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tcal0[i] ;
     4042          tcal.push_back( v ) ;
     4043        }
     4044      }
     4045    }
     4046    else {
     4047      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     4048    }
     4049    return tcal ;
     4050  }
     4051}
     4052
     4053vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
     4054                                            CountedPtr<Scantable>& off,
     4055                                            CountedPtr<Scantable>& sky,
     4056                                            CountedPtr<Scantable>& hot,
     4057                                            CountedPtr<Scantable>& cold,
     4058                                            int index,
     4059                                            string antname )
     4060{
     4061  string reftime = on->getTime( index ) ;
     4062  vector<int> ii( 1, on->getIF( index ) ) ;
     4063  vector<int> ib( 1, on->getBeam( index ) ) ;
     4064  vector<int> ip( 1, on->getPol( index ) ) ;
     4065  vector<int> ic( 1, on->getScan( index ) ) ;
     4066  STSelector sel = STSelector() ;
     4067  sel.setIFs( ii ) ;
     4068  sel.setBeams( ib ) ;
     4069  sel.setPolarizations( ip ) ;
     4070  sky->setSelection( sel ) ;
     4071  hot->setSelection( sel ) ;
     4072  cold->setSelection( sel ) ;
     4073  off->setSelection( sel ) ;
     4074  vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
     4075  vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
     4076  vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
     4077  vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
     4078  vector<float> spec = on->getSpectrum( index ) ;
     4079  vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
     4080  vector<float> sp( tcal.size() ) ;
     4081  if ( antname.find( "APEX" ) != string::npos ) {
     4082    // using gain array
     4083    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4084      float v = ( ( spec[j] - spoff[j] ) / spoff[j] )
     4085        * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
     4086      sp[j] = v ;
     4087    }
     4088  }
     4089  else if ( antname.find( "ALMA" ) != string::npos ) {
     4090    // two-load calibration
     4091    // from equation 5 and 12 of ALMA memo 318 (Mangum 2000)
     4092    //
     4093    // 2009/09/09 Takeshi Nakazato
     4094    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4095      //
     4096      // in case that Tcal is assumed as signal gain
     4097      //
     4098      //float K = ( sphot[j] - spcold[j] ) / ( thot - tcold ) ;
     4099      //float v = ( spec[j] - spoff[j] ) * exp( tsig ) / ( K * tcal[j] * eta ) ;
     4100      //
     4101      // in case that Tcal is defined as
     4102      //
     4103      //    Tcal = ( K * Gs * eta_l * exp( - tau_s ) )^-1
     4104      //
     4105      float v = tcal[j] * ( spec[j] - spsky[j] ) ;
     4106      sp[j] = v ;
     4107    }
     4108  }
     4109  else {
     4110    // Chopper-Wheel calibration (Ulich & Haas 1976)
     4111    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4112      float v = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
     4113      sp[j] = v ;
     4114    }
     4115  }
     4116  sel.reset() ;
     4117  sky->unsetSelection() ;
     4118  hot->unsetSelection() ;
     4119  cold->unsetSelection() ;
     4120  off->unsetSelection() ;
     4121
     4122  return sp ;
     4123}
     4124
     4125vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
     4126                                              CountedPtr<Scantable>& ref,
     4127                                              CountedPtr<Scantable>& sky,
     4128                                              CountedPtr<Scantable>& hot,
     4129                                              CountedPtr<Scantable>& cold,
     4130                                              int index )
     4131{
     4132  string reftime = sig->getTime( index ) ;
     4133  vector<int> ii( 1, sig->getIF( index ) ) ;
     4134  vector<int> ib( 1, sig->getBeam( index ) ) ;
     4135  vector<int> ip( 1, sig->getPol( index ) ) ;
     4136  vector<int> ic( 1, sig->getScan( index ) ) ;
     4137  STSelector sel = STSelector() ;
     4138  sel.setIFs( ii ) ;
     4139  sel.setBeams( ib ) ;
     4140  sel.setPolarizations( ip ) ;
     4141  sky->setSelection( sel ) ;
     4142  hot->setSelection( sel ) ;
     4143  cold->setSelection( sel ) ;
     4144  vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
     4145  vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
     4146  vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
     4147  vector<float> spref = ref->getSpectrum( index ) ;
     4148  vector<float> spsig = sig->getSpectrum( index ) ;
     4149  vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
     4150  vector<float> sp( tcal.size() ) ;
     4151  for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4152    float v = tcal[j] * spsky[j] / ( sphot[j] - spsky[j] ) * ( spsig[j] - spref[j] ) / spref[j] ;
     4153    sp[j] = v ;
     4154  }
     4155  sel.reset() ;
     4156  sky->unsetSelection() ;
     4157  hot->unsetSelection() ;
     4158  cold->unsetSelection() ;
     4159
     4160  return sp ;
     4161}
     4162
     4163vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
     4164                                              CountedPtr<Scantable>& ref,
     4165                                              vector< CountedPtr<Scantable> >& sky,
     4166                                              vector< CountedPtr<Scantable> >& hot,
     4167                                              vector< CountedPtr<Scantable> >& cold,
     4168                                              int index )
     4169{
     4170  string reftime = sig->getTime( index ) ;
     4171  vector<int> ii( 1, sig->getIF( index ) ) ;
     4172  vector<int> ib( 1, sig->getBeam( index ) ) ;
     4173  vector<int> ip( 1, sig->getPol( index ) ) ;
     4174  vector<int> ic( 1, sig->getScan( index ) ) ;
     4175  STSelector sel = STSelector() ;
     4176  sel.setIFs( ii ) ;
     4177  sel.setBeams( ib ) ;
     4178  sel.setPolarizations( ip ) ;
     4179  sky[0]->setSelection( sel ) ;
     4180  hot[0]->setSelection( sel ) ;
     4181  cold[0]->setSelection( sel ) ;
     4182  vector<float> spskys = getSpectrumFromTime( reftime, sky[0], "linear" ) ;
     4183  vector<float> sphots = getSpectrumFromTime( reftime, hot[0], "linear" ) ;
     4184  vector<float> spcolds = getSpectrumFromTime( reftime, cold[0], "linear" ) ;
     4185  vector<float> tcals = getTcalFromTime( reftime, sky[0], "linear" ) ;
     4186  sel.reset() ;
     4187  ii[0] = ref->getIF( index ) ;
     4188  sel.setIFs( ii ) ;
     4189  sel.setBeams( ib ) ;
     4190  sel.setPolarizations( ip ) ;
     4191  sky[1]->setSelection( sel ) ;
     4192  hot[1]->setSelection( sel ) ;
     4193  cold[1]->setSelection( sel ) ;
     4194  vector<float> spskyr = getSpectrumFromTime( reftime, sky[1], "linear" ) ;
     4195  vector<float> sphotr = getSpectrumFromTime( reftime, hot[1], "linear" ) ;
     4196  vector<float> spcoldr = getSpectrumFromTime( reftime, cold[1], "linear" ) ;
     4197  vector<float> tcalr = getTcalFromTime( reftime, sky[1], "linear" ) ; 
     4198  vector<float> spref = ref->getSpectrum( index ) ;
     4199  vector<float> spsig = sig->getSpectrum( index ) ;
     4200  vector<float> sp( tcals.size() ) ;
     4201  for ( unsigned int j = 0 ; j < tcals.size() ; j++ ) {
     4202    float v = tcals[j] * spsig[j] / ( sphots[j] - spskys[j] ) - tcalr[j] * spref[j] / ( sphotr[j] - spskyr[j] ) ;
     4203    sp[j] = v ;
     4204  }
     4205  sel.reset() ;
     4206  sky[0]->unsetSelection() ;
     4207  hot[0]->unsetSelection() ;
     4208  cold[0]->unsetSelection() ;
     4209  sky[1]->unsetSelection() ;
     4210  hot[1]->unsetSelection() ;
     4211  cold[1]->unsetSelection() ;
     4212
     4213  return sp ;
     4214}

  • branches/alma/src/STMath.h

    r1609 r1633  
    200200
    201201  /**
     202   * Calibrate data with Chopper-Wheel like calibration method
     203   * which adopts position switching by antenna motion,
     204   * wobbler (nutator) switching and On-The-Fly observation.
     205   *
     206   * The method is applicable to APEX, and other telescopes other than GBT.
     207   *
     208   * @param a Scantable which contains ON and OFF scans
     209   * @param a string that indicates calibration mode
     210   * @param a string that indicates antenna name
     211   **/
     212  casa::CountedPtr<Scantable> cwcal( const casa::CountedPtr<Scantable>& s,
     213                                       const casa::String calmode,
     214                                       const casa::String antname );
     215
     216  /**
     217   * Calibrate frequency switched scans with Chopper-Wheel like
     218   * calibration method.
     219   *
     220   * The method is applicable to APEX, and other telescopes other than GBT.
     221   *
     222   * @param a Scantable which contains ON and OFF scans
     223   * @param a string that indicates antenna name
     224   **/
     225  casa::CountedPtr<Scantable> cwcalfs( const casa::CountedPtr<Scantable>& s,
     226                                       const casa::String antname );
     227
     228
     229  /**
    202230   * Folding frequency-switch data
    203231   * @param sig
    … …  
    207235  casa::CountedPtr<Scantable> dofold( const casa::CountedPtr<Scantable> &sig,
    208236                                      const casa::CountedPtr<Scantable> &ref,
    209                                       casa::Int choffset );
     237                                      casa::Double choffset,
     238                                      casa::Double choffset = 0.0 );
    210239
    211240  casa::CountedPtr<Scantable>
    … …  
    330359    flagsFromMA(const casa::MaskedArray<casa::Float>& ma);
    331360
     361  vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
     362  vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
     363  vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;
     364  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
     365                                      casa::CountedPtr<Scantable>& off,
     366                                      casa::CountedPtr<Scantable>& sky,
     367                                      casa::CountedPtr<Scantable>& hot,
     368                                      casa::CountedPtr<Scantable>& cold,
     369                                      int index,
     370                                      string antname ) ;
     371  vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
     372                                        casa::CountedPtr<Scantable>& ref,
     373                                        casa::CountedPtr<Scantable>& sky,
     374                                        casa::CountedPtr<Scantable>& hot,
     375                                        casa::CountedPtr<Scantable>& cold,
     376                                        int index ) ;
     377  vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
     378                                        casa::CountedPtr<Scantable>& ref,
     379                                        vector< casa::CountedPtr<Scantable> >& sky,
     380                                        vector< casa::CountedPtr<Scantable> >& hot,
     381                                        vector< casa::CountedPtr<Scantable> >& cold,
     382                                        int index ) ;
     383  double getMJD( string strtime ) ;
     384
    332385  bool insitu_;
    333386};

  • branches/alma/src/STMathWrapper.h

    r1516 r1633  
    206206  }
    207207
     208  // cwcal
     209  ScantableWrapper cwcal( const ScantableWrapper &in,
     210                          const std::string calmode,
     211                          const std::string antname )
     212  {
     213    casa::CountedPtr<Scantable> tab = in.getCP() ;
     214    casa::String mode( calmode ) ;
     215    casa::String name( antname ) ;
     216    return ScantableWrapper( STMath::cwcal( tab, mode, name ) ) ;
     217  }
    208218};
    209219

  • branches/alma/src/python_STMath.cpp

    r1516 r1633  
    7676        // testing average spectra with different channel/resolution
    7777        .def("_new_average", &STMathWrapper::new_average)
     78        // cwcal
     79        .def("cwcal", &STMathWrapper::cwcal)
    7880          ;
    7981    };
Note: See TracChangeset for help on using the changeset viewer.