Changeset 1673 for branches/alma/src


Ignore:
Timestamp:
01/13/10 12:50:13 (15 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-1799

Ready to Release: No

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...

Calibration algorithm for ALMA data (data from OSF) is updated.


Location:
branches/alma/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/src/STMath.cpp

    r1658 r1673  
    5959// tolerance for direction comparison (rad)
    6060#define TOL_OTF    1.0e-15
    61 #define TOL_POINT  9.6963e-5  // 20 arcsec
     61#define TOL_POINT  2.9088821e-4  // 1 arcmin
    6262
    6363STMath::STMath(bool insitu) :
    … …  
    14481448    // change FREQ_ID to unshifted IF setting (only for APEX?)
    14491449    if ( choffset2 != 0.0 ) {
    1450       int freqid = fidColOut( 0 ) ; // assume single-IF data
     1450      uInt freqid = fidColOut( 0 ) ; // assume single-IF data
    14511451      double refpix, refval, increment ;
    14521452      out->frequencies().getEntry( refpix, refval, increment, freqid ) ;
    … …  
    24992499      Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
    25002500      Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
    2501       for (int k=0; k < flag.nelements(); ++k ) {
     2501      for (unsigned int k=0; k < flag.nelements(); ++k ) {
    25022502        if (flag[k] > 0) {
    25032503          spec[k] = 0.0;
    … …  
    34043404}
    34053405 
     3406CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
     3407                                       const String calmode )
     3408{
     3409  // frequency switch
     3410  if ( calmode == "fs" ) {
     3411    return almacalfs( s ) ;
     3412  }
     3413  else {
     3414    string onstr = "*_" ;
     3415    string offstr = "*_" ;
     3416   
     3417    if ( calmode == "ps" || calmode == "otf" ) {
     3418      onstr += "pson" ;
     3419      offstr += "psoff" ;
     3420    }
     3421    else if ( calmode == "wob" ) {
     3422      onstr += "wobon" ;
     3423      offstr += "woboff" ;
     3424    }
     3425   
     3426    vector<bool> masks = s->getMask( 0 ) ;
     3427   
     3428    // off scan
     3429    STSelector sel = STSelector() ;
     3430    sel.setName( offstr ) ;
     3431    s->setSelection( sel ) ;
     3432    // TODO 2010/01/08 TN
     3433    // Grouping by time should be needed before averaging.
     3434    // Each group must have own unique SCANNO (should be renumbered).
     3435    // See PIPELINE/SDCalibration.py
     3436    CountedPtr<Scantable> soff = getScantable( s, false ) ;
     3437    Table ttab = soff->table() ;
     3438    ROScalarColumn<Double> timeCol( ttab, "TIME" ) ;
     3439    uInt nrow = timeCol.nrow() ;
     3440    Vector<Double> timeSep( nrow - 1 ) ;
     3441    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
     3442      timeSep[i] = timeCol(i+1) - timeCol(i) ;
     3443    }
     3444    ScalarColumn<Double> intervalCol( ttab, "INTERVAL" ) ;
     3445    Vector<Double> interval = intervalCol.getColumn() ;
     3446    interval /= 86400.0 ;
     3447    ScalarColumn<uInt> scanCol( ttab, "SCANNO" ) ;
     3448    vector<uInt> glist ;
     3449    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
     3450      double gap = 2.0 * timeSep[i] / ( interval[i] + interval[i+1] ) ;
     3451      //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
     3452      if ( gap > 1.1 ) {
     3453        glist.push_back( i ) ;
     3454      }
     3455    }
     3456    Vector<uInt> gaplist( glist ) ;
     3457    //cout << "gaplist = " << gaplist << endl ;
     3458    uInt newid = 0 ;
     3459    for ( uInt i = 0 ; i < nrow ; i++ ) {
     3460      scanCol.put( i, newid ) ;
     3461      if ( i == gaplist[newid] ) {
     3462        newid++ ;
     3463      }
     3464    }
     3465    //cout << "new scancol = " << scanCol.getColumn() << endl ;
     3466    vector< CountedPtr<Scantable> > tmp( 1, soff ) ;
     3467    CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
     3468    //cout << "aoff.nrow = " << aoff->nrow() << endl ;
     3469    s->unsetSelection() ;
     3470    sel.reset() ;
     3471   
     3472    // on scan
     3473    bool insitu = insitu_ ;
     3474    insitu_ = false ;
     3475    CountedPtr<Scantable> out = getScantable( s, true ) ;
     3476    insitu_ = insitu ;
     3477    sel.setName( onstr ) ;
     3478    s->setSelection( sel ) ;
     3479    TableCopy::copyRows( out->table(), s->table() ) ;
     3480    s->unsetSelection() ;
     3481    sel.reset() ;
     3482   
     3483    // process each on scan
     3484    ArrayColumn<Float> tsysCol ;
     3485    tsysCol.attach( out->table(), "TSYS" ) ;
     3486    for ( int i = 0 ; i < out->nrow() ; i++ ) {
     3487      vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
     3488      out->setSpectrum( sp, i ) ;
     3489    }
     3490
     3491    // remove additional string from SRCNAME
     3492    ScalarColumn<String> srcnameCol ;
     3493    srcnameCol.attach( out->table(), "SRCNAME" ) ;
     3494    Vector<String> srcnames( srcnameCol.getColumn() ) ;
     3495    for ( uInt i = 0 ; i < srcnames.nelements() ; i++ ) {
     3496      srcnames[i] = srcnames[i].substr( 0, srcnames[i].find( onstr.substr(1,onstr.size()-1) ) ) ;
     3497    }
     3498    srcnameCol.putColumn( srcnames ) ;
     3499
     3500    // flux unit
     3501    out->setFluxUnit( "K" ) ;
     3502
     3503    return out ;
     3504  }
     3505}
     3506
    34063507CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
    34073508                                       const String antname )
    … …  
    37673868  ScalarColumn<uInt> sigfidCol ;
    37683869  ScalarColumn<uInt> reffidCol ;
    3769   uInt nchan = (uInt)ssig->nchan() ;
     3870  Int nchan = (Int)ssig->nchan() ;
    37703871  sigifnoCol.attach( sigtab, "IFNO" ) ;
    37713872  refifnoCol.attach( reftab, "IFNO" ) ;
    … …  
    38523953}
    38533954
     3955CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
     3956{
     3957  CountedPtr<Scantable> out ;
     3958
     3959  return out ;
     3960}
     3961
    38543962vector<float> STMath::getSpectrumFromTime( string reftime,
    38553963                                           CountedPtr<Scantable>& s,
    … …  
    38643972  }
    38653973  else if ( s->nrow() == 1 ) {
    3866     os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
     3974    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
    38673975    return s->getSpectrum( 0 ) ;
    38683976  }
    … …  
    38783986        id = idx[1] ;
    38793987      }
    3880       os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     3988      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    38813989      sp = s->getSpectrum( id ) ;
    38823990    }
    … …  
    38903998        id = idx[1] ;
    38913999      }
    3892       os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4000      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    38934001      sp = s->getSpectrum( id ) ;
    38944002    }
    … …  
    39154023        }
    39164024      }
    3917       os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4025      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    39184026      sp = s->getSpectrum( id ) ;     
    39194027    }
    … …  
    39234031        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
    39244032        int id = idx[1] ;
    3925         os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4033        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    39264034        sp = s->getSpectrum( id ) ;
    39274035      }
    … …  
    39304038        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
    39314039        int id = idx[0] ;
    3932         os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4040        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    39334041        sp = s->getSpectrum( id ) ;
    39344042      }
    39354043      else if ( idx[0] == idx[1] ) {
    39364044        // use before
    3937         os << "No need to interporate." << LogIO::POST ;
     4045        //os << "No need to interporate." << LogIO::POST ;
    39384046        int id = idx[0] ;
    3939         os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4047        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
    39404048        sp = s->getSpectrum( id ) ;
    39414049      }
    39424050      else {
    39434051        // do interpolation
    3944         os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4052        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
    39454053        double t0 = getMJD( s->getTime( idx[0] ) ) ;
    39464054        double t1 = getMJD( s->getTime( idx[1] ) ) ;
    … …  
    40384146  else if ( s->nrow() == 1 ) {
    40394147    uInt tcalid = s->getTcalId( 0 ) ;
    4040     os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4148    //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    40414149    tcalTable.getEntry( time, tcalval, tcalid ) ;
    40424150    tcalval.tovector( tcal ) ;
    … …  
    40554163      }
    40564164      uInt tcalid = s->getTcalId( id ) ;
    4057       os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4165      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    40584166      tcalTable.getEntry( time, tcalval, tcalid ) ;
    40594167      tcalval.tovector( tcal ) ;
    … …  
    40694177      }
    40704178      uInt tcalid = s->getTcalId( id ) ;
    4071       os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4179      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    40724180      tcalTable.getEntry( time, tcalval, tcalid ) ;
    40734181      tcalval.tovector( tcal ) ;
    … …  
    40964204      }
    40974205      uInt tcalid = s->getTcalId( id ) ;
    4098       os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4206      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    40994207      tcalTable.getEntry( time, tcalval, tcalid ) ;
    41004208      tcalval.tovector( tcal ) ;
    … …  
    41064214        int id = idx[1] ;
    41074215        uInt tcalid = s->getTcalId( id ) ;
    4108         os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4216        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    41094217        tcalTable.getEntry( time, tcalval, tcalid ) ;
    41104218        tcalval.tovector( tcal ) ;
    … …  
    41154223        int id = idx[0] ;
    41164224        uInt tcalid = s->getTcalId( id ) ;
    4117         os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4225        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    41184226        tcalTable.getEntry( time, tcalval, tcalid ) ;
    41194227        tcalval.tovector( tcal ) ;
    … …  
    41214229      else if ( idx[0] == idx[1] ) {
    41224230        // use before
    4123         os << "No need to interporate." << LogIO::POST ;
     4231        //os << "No need to interporate." << LogIO::POST ;
    41244232        int id = idx[0] ;
    41254233        uInt tcalid = s->getTcalId( id ) ;
    4126         os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4234        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
    41274235        tcalTable.getEntry( time, tcalval, tcalid ) ;
    41284236        tcalval.tovector( tcal ) ;
    … …  
    41304238      else {
    41314239        // do interpolation
    4132         os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4240        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
    41334241        double t0 = getMJD( s->getTime( idx[0] ) ) ;
    41344242        double t1 = getMJD( s->getTime( idx[1] ) ) ;
    … …  
    41704278  }
    41714279  else if ( s->nrow() == 1 ) {
    4172     os << "use row " << 0 << LogIO::POST ;
     4280    //os << "use row " << 0 << LogIO::POST ;
    41734281    tsysval = tsysCol( 0 ) ;
    41744282    tsysval.tovector( tsys ) ;
    … …  
    41864294        id = idx[1] ;
    41874295      }
    4188       os << "use row " << id << LogIO::POST ;
     4296      //os << "use row " << id << LogIO::POST ;
    41894297      tsysval = tsysCol( id ) ;
    41904298      tsysval.tovector( tsys ) ;
    … …  
    41994307        id = idx[1] ;
    42004308      }
    4201       os << "use row " << id << LogIO::POST ;
     4309      //os << "use row " << id << LogIO::POST ;
    42024310      tsysval = tsysCol( id ) ;
    42034311      tsysval.tovector( tsys ) ;
    … …  
    42254333        }
    42264334      }
    4227       os << "use row " << id << LogIO::POST ;
     4335      //os << "use row " << id << LogIO::POST ;
    42284336      tsysval = tsysCol( id ) ;
    42294337      tsysval.tovector( tsys ) ;
    … …  
    42344342        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
    42354343        int id = idx[1] ;
    4236         os << "use row " << id << LogIO::POST ;
     4344        //os << "use row " << id << LogIO::POST ;
    42374345        tsysval = tsysCol( id ) ;
    42384346        tsysval.tovector( tsys ) ;
    … …  
    42424350        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
    42434351        int id = idx[0] ;
    4244         os << "use row " << id << LogIO::POST ;
     4352        //os << "use row " << id << LogIO::POST ;
    42454353        tsysval = tsysCol( id ) ;
    42464354        tsysval.tovector( tsys ) ;
    … …  
    42484356      else if ( idx[0] == idx[1] ) {
    42494357        // use before
    4250         os << "No need to interporate." << LogIO::POST ;
     4358        //os << "No need to interporate." << LogIO::POST ;
    42514359        int id = idx[0] ;
    4252         os << "use row " << id << LogIO::POST ;
     4360        //os << "use row " << id << LogIO::POST ;
    42534361        tsysval = tsysCol( id ) ;
    42544362        tsysval.tovector( tsys ) ;
    … …  
    42564364      else {
    42574365        // do interpolation
    4258         os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4366        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
    42594367        double t0 = getMJD( s->getTime( idx[0] ) ) ;
    42604368        double t1 = getMJD( s->getTime( idx[1] ) ) ;
    … …  
    43464454  hot->unsetSelection() ;
    43474455  cold->unsetSelection() ;
     4456  off->unsetSelection() ;
     4457
     4458  return sp ;
     4459}
     4460
     4461vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
     4462                                            CountedPtr<Scantable>& off,
     4463                                            int index )
     4464{
     4465  string reftime = on->getTime( index ) ;
     4466  vector<int> ii( 1, on->getIF( index ) ) ;
     4467  vector<int> ib( 1, on->getBeam( index ) ) ;
     4468  vector<int> ip( 1, on->getPol( index ) ) ;
     4469  vector<int> ic( 1, on->getScan( index ) ) ;
     4470  STSelector sel = STSelector() ;
     4471  sel.setIFs( ii ) ;
     4472  sel.setBeams( ib ) ;
     4473  sel.setPolarizations( ip ) ;
     4474  off->setSelection( sel ) ;
     4475  vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
     4476  vector<float> spec = on->getSpectrum( index ) ;
     4477  //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
     4478  //vector<float> tsys = on->getTsysVec( index ) ;
     4479  ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
     4480  Vector<Float> tsys = tsysCol( index ) ;
     4481  vector<float> sp( spec.size() ) ;
     4482  // ALMA Calibration
     4483  //
     4484  // Ta* = Tsys * ( ON - OFF ) / OFF
     4485  //
     4486  // 2010/01/07 Takeshi Nakazato
     4487  unsigned int tsyssize = tsys.nelements() ;
     4488  unsigned int spsize = sp.size() ;
     4489  for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
     4490    float tscale = 0.0 ;
     4491    if ( tsyssize == spsize )
     4492      tscale = tsys[j] ;
     4493    else
     4494      tscale = tsys[0] ;
     4495    float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
     4496    sp[j] = v ;
     4497  }
     4498  sel.reset() ;
    43484499  off->unsetSelection() ;
    43494500

  • branches/alma/src/STMath.h

    r1652 r1673  
    238238                                      casa::Double choffset = 0.0 );
    239239
     240  /**
     241   * ALMA calibration
     242   **/
     243  casa::CountedPtr<Scantable> almacal( const casa::CountedPtr<Scantable>& s,
     244                                       const casa::String calmode ) ;
     245  casa::CountedPtr<Scantable> almacalfs( const casa::CountedPtr<Scantable>& s ) ;
     246
    240247  casa::CountedPtr<Scantable>
    241248    freqSwitch( const casa::CountedPtr<Scantable>& in );
    … …  
    362369  vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
    363370  vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;
     371
     372  // Chopper-Wheel type calibration
    364373  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
    365374                                      casa::CountedPtr<Scantable>& off,
    … …  
    369378                                      int index,
    370379                                      string antname ) ;
     380  // Tsys * (ON-OFF)/OFF
     381  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
     382                                      casa::CountedPtr<Scantable>& off,
     383                                      int index ) ;
    371384  vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
    372385                                        casa::CountedPtr<Scantable>& ref,

  • branches/alma/src/STMathWrapper.h

    r1633 r1673  
    216216    return ScantableWrapper( STMath::cwcal( tab, mode, name ) ) ;
    217217  }
     218  // almacal
     219  ScantableWrapper almacal( const ScantableWrapper &in,
     220                          const std::string calmode )
     221  {
     222    casa::CountedPtr<Scantable> tab = in.getCP() ;
     223    casa::String mode( calmode ) ;
     224    return ScantableWrapper( STMath::almacal( tab, mode ) ) ;
     225  }
    218226};
    219227

  • branches/alma/src/python_STMath.cpp

    r1633 r1673  
    7878        // cwcal
    7979        .def("cwcal", &STMathWrapper::cwcal)
     80        .def("almacal", &STMathWrapper::almacal)
    8081          ;
    8182    };
Note: See TracChangeset for help on using the changeset viewer.