Changeset 1673 for branches/alma
- Timestamp:
- 01/13/10 12:50:13 (15 years ago)
- Location:
- branches/alma
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/alma/python/asapmath.py
r1659 r1673 872 872 asaplog.push( 'Calibrating %s position-switched data.' % antname ) 873 873 print_log() 874 if ( antname.find( 'APEX' ) != -1 or antname.find( 'ALMA' ) != -1 ): 874 if ( antname.find( 'APEX' ) != -1 ): 875 scal = apexcal( scantab, scannos, calmode, verify ) 876 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ): 875 877 scal = almacal( scantab, scannos, calmode, verify ) 876 878 else: … … 879 881 asaplog.push( 'Calibrating %s frequency-switched data.' % antname ) 880 882 print_log() 881 if ( antname.find( 'APEX' ) != -1 or antname.find( 'ALMA' ) != -1 ): 883 if ( antname.find( 'APEX' ) != -1 ): 884 scal = apexcal( scantab, scannos, calmode, verify ) 885 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ): 882 886 scal = almacal( scantab, scannos, calmode, verify ) 883 887 else: … … 893 897 return scal 894 898 895 def a lmacal( scantab, scannos=[], calmode='none', verify=False ):896 """ 897 Calibrate APEX and ALMAdata899 def apexcal( scantab, scannos=[], calmode='none', verify=False ): 900 """ 901 Calibrate APEX data 898 902 899 903 Parameters: … … 909 913 ssub = scantab.get_scan( scannos ) 910 914 scal = scantable( stm.cwcal( ssub, calmode, antname ) ) 915 return scal 916 917 def almacal( scantab, scannos=[], calmode='none', verify=False ): 918 """ 919 Calibrate ALMA data 920 921 Parameters: 922 scantab: scantable 923 scannos: list of scan number 924 calmode: calibration mode 925 926 verify: verify calibration 927 """ 928 from asap._asap import stmath 929 stm = stmath() 930 ssub = scantab.get_scan( scannos ) 931 scal = scantable( stm.almacal( ssub, calmode ) ) 911 932 return scal 912 933 -
branches/alma/python/scantable.py
r1666 r1673 38 38 average = rcParams['scantable.autoaverage'] 39 39 if getpt is None: 40 getpt = False40 getpt = True 41 41 varlist = vars() 42 42 from asap._asap import stmath -
branches/alma/src/STMath.cpp
r1658 r1673 59 59 // tolerance for direction comparison (rad) 60 60 #define TOL_OTF 1.0e-15 61 #define TOL_POINT 9.6963e-5 // 20 arcsec61 #define TOL_POINT 2.9088821e-4 // 1 arcmin 62 62 63 63 STMath::STMath(bool insitu) : … … 1448 1448 // change FREQ_ID to unshifted IF setting (only for APEX?) 1449 1449 if ( choffset2 != 0.0 ) { 1450 int freqid = fidColOut( 0 ) ; // assume single-IF data1450 uInt freqid = fidColOut( 0 ) ; // assume single-IF data 1451 1451 double refpix, refval, increment ; 1452 1452 out->frequencies().getEntry( refpix, refval, increment, freqid ) ; … … 2499 2499 Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5); 2500 2500 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 ) { 2502 2502 if (flag[k] > 0) { 2503 2503 spec[k] = 0.0; … … 3404 3404 } 3405 3405 3406 CountedPtr<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 3406 3507 CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s, 3407 3508 const String antname ) … … 3767 3868 ScalarColumn<uInt> sigfidCol ; 3768 3869 ScalarColumn<uInt> reffidCol ; 3769 uInt nchan = (uInt)ssig->nchan() ;3870 Int nchan = (Int)ssig->nchan() ; 3770 3871 sigifnoCol.attach( sigtab, "IFNO" ) ; 3771 3872 refifnoCol.attach( reftab, "IFNO" ) ; … … 3852 3953 } 3853 3954 3955 CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s ) 3956 { 3957 CountedPtr<Scantable> out ; 3958 3959 return out ; 3960 } 3961 3854 3962 vector<float> STMath::getSpectrumFromTime( string reftime, 3855 3963 CountedPtr<Scantable>& s, … … 3864 3972 } 3865 3973 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 ; 3867 3975 return s->getSpectrum( 0 ) ; 3868 3976 } … … 3878 3986 id = idx[1] ; 3879 3987 } 3880 os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;3988 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 3881 3989 sp = s->getSpectrum( id ) ; 3882 3990 } … … 3890 3998 id = idx[1] ; 3891 3999 } 3892 os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4000 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 3893 4001 sp = s->getSpectrum( id ) ; 3894 4002 } … … 3915 4023 } 3916 4024 } 3917 os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;4025 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ; 3918 4026 sp = s->getSpectrum( id ) ; 3919 4027 } … … 3923 4031 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ; 3924 4032 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 ; 3926 4034 sp = s->getSpectrum( id ) ; 3927 4035 } … … 3930 4038 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ; 3931 4039 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 ; 3933 4041 sp = s->getSpectrum( id ) ; 3934 4042 } 3935 4043 else if ( idx[0] == idx[1] ) { 3936 4044 // use before 3937 os << "No need to interporate." << LogIO::POST ;4045 //os << "No need to interporate." << LogIO::POST ; 3938 4046 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 ; 3940 4048 sp = s->getSpectrum( id ) ; 3941 4049 } 3942 4050 else { 3943 4051 // 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 ; 3945 4053 double t0 = getMJD( s->getTime( idx[0] ) ) ; 3946 4054 double t1 = getMJD( s->getTime( idx[1] ) ) ; … … 4038 4146 else if ( s->nrow() == 1 ) { 4039 4147 uInt tcalid = s->getTcalId( 0 ) ; 4040 os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;4148 //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4041 4149 tcalTable.getEntry( time, tcalval, tcalid ) ; 4042 4150 tcalval.tovector( tcal ) ; … … 4055 4163 } 4056 4164 uInt tcalid = s->getTcalId( id ) ; 4057 os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4165 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4058 4166 tcalTable.getEntry( time, tcalval, tcalid ) ; 4059 4167 tcalval.tovector( tcal ) ; … … 4069 4177 } 4070 4178 uInt tcalid = s->getTcalId( id ) ; 4071 os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4179 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4072 4180 tcalTable.getEntry( time, tcalval, tcalid ) ; 4073 4181 tcalval.tovector( tcal ) ; … … 4096 4204 } 4097 4205 uInt tcalid = s->getTcalId( id ) ; 4098 os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4206 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4099 4207 tcalTable.getEntry( time, tcalval, tcalid ) ; 4100 4208 tcalval.tovector( tcal ) ; … … 4106 4214 int id = idx[1] ; 4107 4215 uInt tcalid = s->getTcalId( id ) ; 4108 os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4216 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4109 4217 tcalTable.getEntry( time, tcalval, tcalid ) ; 4110 4218 tcalval.tovector( tcal ) ; … … 4115 4223 int id = idx[0] ; 4116 4224 uInt tcalid = s->getTcalId( id ) ; 4117 os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4225 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4118 4226 tcalTable.getEntry( time, tcalval, tcalid ) ; 4119 4227 tcalval.tovector( tcal ) ; … … 4121 4229 else if ( idx[0] == idx[1] ) { 4122 4230 // use before 4123 os << "No need to interporate." << LogIO::POST ;4231 //os << "No need to interporate." << LogIO::POST ; 4124 4232 int id = idx[0] ; 4125 4233 uInt tcalid = s->getTcalId( id ) ; 4126 os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;4234 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ; 4127 4235 tcalTable.getEntry( time, tcalval, tcalid ) ; 4128 4236 tcalval.tovector( tcal ) ; … … 4130 4238 else { 4131 4239 // 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 ; 4133 4241 double t0 = getMJD( s->getTime( idx[0] ) ) ; 4134 4242 double t1 = getMJD( s->getTime( idx[1] ) ) ; … … 4170 4278 } 4171 4279 else if ( s->nrow() == 1 ) { 4172 os << "use row " << 0 << LogIO::POST ;4280 //os << "use row " << 0 << LogIO::POST ; 4173 4281 tsysval = tsysCol( 0 ) ; 4174 4282 tsysval.tovector( tsys ) ; … … 4186 4294 id = idx[1] ; 4187 4295 } 4188 os << "use row " << id << LogIO::POST ;4296 //os << "use row " << id << LogIO::POST ; 4189 4297 tsysval = tsysCol( id ) ; 4190 4298 tsysval.tovector( tsys ) ; … … 4199 4307 id = idx[1] ; 4200 4308 } 4201 os << "use row " << id << LogIO::POST ;4309 //os << "use row " << id << LogIO::POST ; 4202 4310 tsysval = tsysCol( id ) ; 4203 4311 tsysval.tovector( tsys ) ; … … 4225 4333 } 4226 4334 } 4227 os << "use row " << id << LogIO::POST ;4335 //os << "use row " << id << LogIO::POST ; 4228 4336 tsysval = tsysCol( id ) ; 4229 4337 tsysval.tovector( tsys ) ; … … 4234 4342 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ; 4235 4343 int id = idx[1] ; 4236 os << "use row " << id << LogIO::POST ;4344 //os << "use row " << id << LogIO::POST ; 4237 4345 tsysval = tsysCol( id ) ; 4238 4346 tsysval.tovector( tsys ) ; … … 4242 4350 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ; 4243 4351 int id = idx[0] ; 4244 os << "use row " << id << LogIO::POST ;4352 //os << "use row " << id << LogIO::POST ; 4245 4353 tsysval = tsysCol( id ) ; 4246 4354 tsysval.tovector( tsys ) ; … … 4248 4356 else if ( idx[0] == idx[1] ) { 4249 4357 // use before 4250 os << "No need to interporate." << LogIO::POST ;4358 //os << "No need to interporate." << LogIO::POST ; 4251 4359 int id = idx[0] ; 4252 os << "use row " << id << LogIO::POST ;4360 //os << "use row " << id << LogIO::POST ; 4253 4361 tsysval = tsysCol( id ) ; 4254 4362 tsysval.tovector( tsys ) ; … … 4256 4364 else { 4257 4365 // 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 ; 4259 4367 double t0 = getMJD( s->getTime( idx[0] ) ) ; 4260 4368 double t1 = getMJD( s->getTime( idx[1] ) ) ; … … 4346 4454 hot->unsetSelection() ; 4347 4455 cold->unsetSelection() ; 4456 off->unsetSelection() ; 4457 4458 return sp ; 4459 } 4460 4461 vector<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() ; 4348 4499 off->unsetSelection() ; 4349 4500 -
branches/alma/src/STMath.h
r1652 r1673 238 238 casa::Double choffset = 0.0 ); 239 239 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 240 247 casa::CountedPtr<Scantable> 241 248 freqSwitch( const casa::CountedPtr<Scantable>& in ); … … 362 369 vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ; 363 370 vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ; 371 372 // Chopper-Wheel type calibration 364 373 vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on, 365 374 casa::CountedPtr<Scantable>& off, … … 369 378 int index, 370 379 string antname ) ; 380 // Tsys * (ON-OFF)/OFF 381 vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on, 382 casa::CountedPtr<Scantable>& off, 383 int index ) ; 371 384 vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig, 372 385 casa::CountedPtr<Scantable>& ref, -
branches/alma/src/STMathWrapper.h
r1633 r1673 216 216 return ScantableWrapper( STMath::cwcal( tab, mode, name ) ) ; 217 217 } 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 } 218 226 }; 219 227 -
branches/alma/src/python_STMath.cpp
r1633 r1673 78 78 // cwcal 79 79 .def("cwcal", &STMathWrapper::cwcal) 80 .def("almacal", &STMathWrapper::almacal) 80 81 ; 81 82 };
Note:
See TracChangeset
for help on using the changeset viewer.