Changeset 1634 for branches/alma/src


Ignore:
Timestamp:
09/18/09 18:42:35 (15 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-1422

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix

  1. Changed code to update Tsys after calibration
  2. Changed code to update FluxUnit after calibration
  3. Removed additional log messages


Location:
branches/alma/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/src/STMath.cpp

    r1633 r1634  
    33683368   
    33693369    // process each on scan
     3370    ArrayColumn<Float> tsysCol ;
     3371    tsysCol.attach( out->table(), "TSYS" ) ;
    33703372    for ( int i = 0 ; i < out->nrow() ; i++ ) {
    33713373      vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ;
    33723374      out->setSpectrum( sp, i ) ;
     3375      string reftime = out->getTime( i ) ;
     3376      vector<int> ii( 1, out->getIF( i ) ) ;
     3377      vector<int> ib( 1, out->getBeam( i ) ) ;
     3378      vector<int> ip( 1, out->getPol( i ) ) ;
     3379      sel.setIFs( ii ) ;
     3380      sel.setBeams( ib ) ;
     3381      sel.setPolarizations( ip ) ;
     3382      asky->setSelection( sel ) ;   
     3383      vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
     3384      const Vector<Float> Vtsys( sptsys ) ;
     3385      tsysCol.put( i, Vtsys ) ;
     3386      asky->unsetSelection() ;
     3387      sel.reset() ;
    33733388    }
    33743389
    … …  
    33813396    }
    33823397    srcnameCol.putColumn( srcnames ) ;
     3398
     3399    // flux unit
     3400    out->setFluxUnit( "K" ) ;
    33833401
    33843402    return out ;
    … …  
    34963514
    34973515    // process each sig and ref scan
     3516    ArrayColumn<Float> tsysCollo ;
     3517    tsysCollo.attach( ssig->table(), "TSYS" ) ;
     3518    ArrayColumn<Float> tsysColhi ;
     3519    tsysColhi.attach( sref->table(), "TSYS" ) ;
    34983520    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
    34993521      vector< CountedPtr<Scantable> > sky( 2 ) ;
    … …  
    35083530      vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ;
    35093531      ssig->setSpectrum( sp, i ) ;
     3532      string reftime = ssig->getTime( i ) ;
     3533      vector<int> ii( 1, ssig->getIF( i ) ) ;
     3534      vector<int> ib( 1, ssig->getBeam( i ) ) ;
     3535      vector<int> ip( 1, ssig->getPol( i ) ) ;
     3536      sel.setIFs( ii ) ;
     3537      sel.setBeams( ib ) ;
     3538      sel.setPolarizations( ip ) ;
     3539      askylo->setSelection( sel ) ;
     3540      vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
     3541      const Vector<Float> Vtsyslo( sptsys ) ;
     3542      tsysCollo.put( i, Vtsyslo ) ;
     3543      askylo->unsetSelection() ;
     3544      sel.reset() ;
    35103545      sky[0] = askyhi ;
    35113546      sky[1] = askylo ;
    … …  
    35163551      sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ;
    35173552      sref->setSpectrum( sp, i ) ;
     3553      reftime = sref->getTime( i ) ;
     3554      ii[0] = sref->getIF( i )  ;
     3555      ib[0] = sref->getBeam( i ) ;
     3556      ip[0] = sref->getPol( i ) ;
     3557      sel.setIFs( ii ) ;
     3558      sel.setBeams( ib ) ;
     3559      sel.setPolarizations( ip ) ;
     3560      askyhi->setSelection( sel ) ;   
     3561      sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
     3562      const Vector<Float> Vtsyshi( sptsys ) ;
     3563      tsysColhi.put( i, Vtsyshi ) ;
     3564      askyhi->unsetSelection() ;
     3565      sel.reset() ;
    35183566    }
    35193567
    … …  
    36003648
    36013649    // process each sig and ref scan
     3650    ArrayColumn<Float> tsysCollo ;
     3651    tsysCollo.attach( ssig->table(), "TSYS" ) ;
     3652    ArrayColumn<Float> tsysColhi ;
     3653    tsysColhi.attach( sref->table(), "TSYS" ) ;
    36023654    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
    36033655      vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ;
    36043656      ssig->setSpectrum( sp, i ) ;
    36053657      sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ;
     3658      string reftime = ssig->getTime( i ) ;
     3659      vector<int> ii( 1, ssig->getIF( i ) ) ;
     3660      vector<int> ib( 1, ssig->getBeam( i ) ) ;
     3661      vector<int> ip( 1, ssig->getPol( i ) ) ;
     3662      sel.setIFs( ii ) ;
     3663      sel.setBeams( ib ) ;
     3664      sel.setPolarizations( ip ) ;
     3665      askylo->setSelection( sel ) ;
     3666      vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
     3667      const Vector<Float> Vtsyslo( sptsys ) ;
     3668      tsysCollo.put( i, Vtsyslo ) ;
     3669      askylo->unsetSelection() ;
     3670      sel.reset() ;
    36063671      sref->setSpectrum( sp, i ) ;
     3672      reftime = sref->getTime( i ) ;
     3673      ii[0] = sref->getIF( i )  ;
     3674      ib[0] = sref->getBeam( i ) ;
     3675      ip[0] = sref->getPol( i ) ;
     3676      sel.setIFs( ii ) ;
     3677      sel.setBeams( ib ) ;
     3678      sel.setPolarizations( ip ) ;
     3679      askyhi->setSelection( sel ) ;   
     3680      sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
     3681      const Vector<Float> Vtsyshi( sptsys ) ;
     3682      tsysColhi.put( i, Vtsyshi ) ;
     3683      askyhi->unsetSelection() ;
     3684      sel.reset() ;
    36073685    }
    36083686  }
    … …  
    36563734
    36573735    // process each sig and ref scan
     3736    ArrayColumn<Float> tsysColsig ;
     3737    tsysColsig.attach( ssig->table(), "TSYS" ) ;
     3738    ArrayColumn<Float> tsysColref ;
     3739    tsysColref.attach( ssig->table(), "TSYS" ) ;
    36583740    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
    36593741      vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ;
    36603742      ssig->setSpectrum( sp, i ) ;
     3743      string reftime = ssig->getTime( i ) ;
     3744      vector<int> ii( 1, ssig->getIF( i ) ) ;
     3745      vector<int> ib( 1, ssig->getBeam( i ) ) ;
     3746      vector<int> ip( 1, ssig->getPol( i ) ) ;
     3747      sel.setIFs( ii ) ;
     3748      sel.setBeams( ib ) ;
     3749      sel.setPolarizations( ip ) ;
     3750      asky->setSelection( sel ) ;
     3751      vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
     3752      const Vector<Float> Vtsys( sptsys ) ;
     3753      tsysColsig.put( i, Vtsys ) ;
     3754      asky->unsetSelection() ;
     3755      sel.reset() ;
    36613756      sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ;
    36623757      sref->setSpectrum( sp, i ) ;
     3758      tsysColref.put( i, Vtsys ) ;
    36633759    }
    36643760  }
    … …  
    37433839  }
    37443840
    3745   out = merge( tmp ) ;
     3841  if ( tmp.size() > 1 ) {
     3842    out = merge( tmp ) ;
     3843  }
     3844  else {
     3845    out = tmp[0] ;
     3846  }
     3847
     3848  // flux unit
     3849  out->setFluxUnit( "K" ) ;
    37463850
    37473851  return out ;
    … …  
    40484152    }
    40494153    return tcal ;
     4154  }
     4155}
     4156
     4157vector<float> STMath::getTsysFromTime( string reftime,
     4158                                       CountedPtr<Scantable>& s,
     4159                                       string mode )
     4160{
     4161  LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ;
     4162  ArrayColumn<Float> tsysCol ;
     4163  tsysCol.attach( s->table(), "TSYS" ) ;
     4164  vector<float> tsys ;
     4165  String time ;
     4166  Vector<Float> tsysval ;
     4167  if ( s->nrow() == 0 ) {
     4168    os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ;
     4169    return tsys ;
     4170  }
     4171  else if ( s->nrow() == 1 ) {
     4172    os << "use row " << 0 << LogIO::POST ;
     4173    tsysval = tsysCol( 0 ) ;
     4174    tsysval.tovector( tsys ) ;
     4175    return tsys ;
     4176  }
     4177  else {
     4178    vector<int> idx = getRowIdFromTime( reftime, s ) ;
     4179    if ( mode == "before" ) {
     4180      int id = -1 ;
     4181      if ( idx[0] != -1 ) {
     4182        id = idx[0] ;
     4183      }
     4184      else if ( idx[1] != -1 ) {
     4185        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
     4186        id = idx[1] ;
     4187      }
     4188      os << "use row " << id << LogIO::POST ;
     4189      tsysval = tsysCol( id ) ;
     4190      tsysval.tovector( tsys ) ;
     4191    }
     4192    else if ( mode == "after" ) {
     4193      int id = -1 ;
     4194      if ( idx[1] != -1 ) {
     4195        id = idx[1] ;
     4196      }
     4197      else if ( idx[0] != -1 ) {
     4198        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
     4199        id = idx[1] ;
     4200      }
     4201      os << "use row " << id << LogIO::POST ;
     4202      tsysval = tsysCol( id ) ;
     4203      tsysval.tovector( tsys ) ;
     4204    }
     4205    else if ( mode == "nearest" ) {
     4206      int id = -1 ;
     4207      if ( idx[0] == -1 ) {
     4208        id = idx[1] ;
     4209      }
     4210      else if ( idx[1] == -1 ) {
     4211        id = idx[0] ;
     4212      }
     4213      else if ( idx[0] == idx[1] ) {
     4214        id = idx[0] ;
     4215      }
     4216      else {
     4217        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4218        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4219        double tref = getMJD( reftime ) ;
     4220        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     4221          id = idx[1] ;
     4222        }
     4223        else {
     4224          id = idx[0] ;
     4225        }
     4226      }
     4227      os << "use row " << id << LogIO::POST ;
     4228      tsysval = tsysCol( id ) ;
     4229      tsysval.tovector( tsys ) ;
     4230    }
     4231    else if ( mode == "linear" ) {
     4232      if ( idx[0] == -1 ) {
     4233        // use after
     4234        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
     4235        int id = idx[1] ;
     4236        os << "use row " << id << LogIO::POST ;
     4237        tsysval = tsysCol( id ) ;
     4238        tsysval.tovector( tsys ) ;
     4239      }
     4240      else if ( idx[1] == -1 ) {
     4241        // use before
     4242        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
     4243        int id = idx[0] ;
     4244        os << "use row " << id << LogIO::POST ;
     4245        tsysval = tsysCol( id ) ;
     4246        tsysval.tovector( tsys ) ;
     4247      }
     4248      else if ( idx[0] == idx[1] ) {
     4249        // use before
     4250        os << "No need to interporate." << LogIO::POST ;
     4251        int id = idx[0] ;
     4252        os << "use row " << id << LogIO::POST ;
     4253        tsysval = tsysCol( id ) ;
     4254        tsysval.tovector( tsys ) ;
     4255      }
     4256      else {
     4257        // do interpolation
     4258        os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4259        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4260        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4261        double tref = getMJD( reftime ) ;
     4262        vector<float> tsys0 ;
     4263        vector<float> tsys1 ;
     4264        tsysval = tsysCol( idx[0] ) ;
     4265        tsysval.tovector( tsys0 ) ;
     4266        tsysval = tsysCol( idx[1] ) ;
     4267        tsysval.tovector( tsys1 ) ;       
     4268        for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {
     4269          float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ;
     4270          tsys.push_back( v ) ;
     4271        }
     4272      }
     4273    }
     4274    else {
     4275      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     4276    }
     4277    return tsys ;
    40504278  }
    40514279}
    … …  
    42134441  return sp ;
    42144442}
     4443
     4444CountedPtr<Scantable> STMath::pressedOut( const casa::CountedPtr<Scantable> &in,
     4445                                          const float numpoly,
     4446                                          const float radius,
     4447                                          const float threshold )
     4448{
     4449  // scan pattern analysis
     4450
     4451
     4452  // for each channel
     4453  int nchan = in->nchan() ;
     4454  for ( int i = 0 ; i < nchan ; i++ ) {
     4455    // smoothing
     4456   
     4457    // get dTij
     4458   
     4459    // polynomial fitting
     4460   
     4461    // remove scanning effect from the map
     4462
     4463  }
     4464
     4465  return in ;
     4466}

  • branches/alma/src/STMath.h

    r1633 r1634  
    317317    throw (casa::AipsError) ;
    318318
     319  /***
     320   * "Pressed-Out" method (Sofe & Reich 1979)
     321   * @param input scantable
     322   * @param order of polynomial fitting function
     323   * @param smoothing beam radius [arcsec]
     324   * @param threshold valus in the unit of sigma (standard deviation)
     325   ***/
     326  casa::CountedPtr<Scantable>
     327    pressedOut( const casa::CountedPtr<Scantable> &in,
     328                const float numpoly,
     329                const float radius,
     330                const float threshold ) ;
    319331
    320332private:
    … …  
    361373  vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
    362374  vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
     375  vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
    363376  vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;
    364377  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,

  • branches/alma/src/Scantable.h

    r1603 r1634  
    309309  float getParAngle(int whichrow) const
    310310    { return paraCol_(whichrow); }
     311  int getTcalId(int whichrow) const
     312    { return mtcalidCol_(whichrow); }
    311313
    312314  std::string getSourceName(int whichrow) const
Note: See TracChangeset for help on using the changeset viewer.