Changeset 2502 for branches/hpc33

Timestamp:

05/11/12 18:38:39 (13 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: No

Ready for Test: 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...

Re-write almacal to improve performance.
Changed data selection.

Location:

branches/hpc33/src

Files:

• STMath.cpp (modified) (5 diffs)
• STMath.h (modified) (1 diff)

branches/hpc33/src/STMath.cpp

@@ -3923 +3923 @@
     
     // process each on scan
-    ArrayColumn<Float> tsysCol ;
-    tsysCol.attach( out->table(), "TSYS" ) ;
-    for ( int i = 0 ; i < out->nrow() ; i++ ) {
-      vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
-      out->setSpectrum( sp, i ) ;
+//     for ( int i = 0 ; i < out->nrow() ; i++ ) {
+//       vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
+//       out->setSpectrum( sp, i ) ;
+//     }
+
+    // selection first, then process
+    vector<unsigned int> ifs = out->getIFNos() ;
+    vector<unsigned int> beams = out->getBeamNos() ;
+    vector<unsigned int> pols = out->getPolNos() ;
+    vector<int> v( 1 ) ;
+    for ( vector<unsigned int>::iterator i = ifs.begin() ; \
+          i != ifs.end() ; i++ ) {
+      v[0] = *i ;
+      sel.setIFs( v ) ;
+      for ( vector<unsigned int>::iterator j = pols.begin() ; \
+            j != pols.end() ; j++ ) {
+        v[0] = *j ;
+        sel.setPolarizations( v ) ;
+        for ( vector<unsigned int>::iterator k = beams.begin() ; \
+              k != beams.end() ; k++ ) {
+          v[0] = *k ;
+          sel.setBeams( v ) ;
+          out->setSelection( sel ) ;
+          aoff->setSelection( sel ) ;
+
+          // calibrate data
+          calibrateALMA( out, aoff ) ;
+
+          out->unsetSelection() ;
+          aoff->unsetSelection() ;
+          sel.reset() ;
+        }
+      }
     }
 
@@ -4575 +4603 @@
         vector<float> sp0 = s->getSpectrum( idx[0] ) ;
         vector<float> sp1 = s->getSpectrum( idx[1] ) ;
+        sp.resize( sp0.size() ) ;
         for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
-          float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
-          sp.push_back( v ) ;
+          sp[i] = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
         }
       }
@@ -4975 +5003 @@
   vector<int> ib( 1, on->getBeam( index ) ) ;
   vector<int> ip( 1, on->getPol( index ) ) ;
-  vector<int> ic( 1, on->getScan( index ) ) ;
   STSelector sel = STSelector() ;
   sel.setIFs( ii ) ;
@@ -5025 +5052 @@
   vector<int> ib( 1, on->getBeam( index ) ) ;
   vector<int> ip( 1, on->getPol( index ) ) ;
-  vector<int> ic( 1, on->getScan( index ) ) ;
   STSelector sel = STSelector() ;
   sel.setIFs( ii ) ;
@@ -5058 +5084 @@
 
   return sp ;
+}
+
+void STMath::calibrateALMA( CountedPtr<Scantable>& on,
+                                        CountedPtr<Scantable>& off )
+{
+//   string reftime = on->getTime( index ) ;
+  ROTableColumn timeCol( on->table(), "TIME" ) ;
+  ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
+  vector<float> sp( on->nchan( on->getIF(0) ) ) ;
+  unsigned int spsize = sp.size() ;
+  for ( int index = 0 ; index < on->nrow() ; index++ ) {
+    double reftime = timeCol.asdouble(index) ;
+    vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
+    vector<float> spec = on->getSpectrum( index ) ;
+    Vector<Float> tsys = tsysCol( index ) ;
+    // ALMA Calibration
+    // 
+    // Ta* = Tsys * ( ON - OFF ) / OFF
+    //
+    // 2010/01/07 Takeshi Nakazato
+    unsigned int tsyssize = tsys.nelements() ;
+    for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
+      float tscale = 0.0 ;
+      if ( tsyssize == spsize ) 
+        tscale = tsys[j] ;
+      else 
+        tscale = tsys[0] ;
+      float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
+      sp[j] = v ;
+    }
+    on->setSpectrum( sp, index ) ;
+  }
 }
 

branches/hpc33/src/STMath.h

@@ -421 +421 @@
                                       casa::CountedPtr<Scantable>& off,
                                       int index ) ;
+  void calibrateALMA( casa::CountedPtr<Scantable>& on,
+                      casa::CountedPtr<Scantable>& off ) ;
   vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
                                         casa::CountedPtr<Scantable>& ref,