Changeset 532 for trunk


Ignore:
Timestamp:
03/02/05 11:03:24 (20 years ago)
Author:
kil064
Message:

add Tsys weighting to average_pol

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/python/scantable.py

    r530 r532  
    988988                         averaging will be applied. The output will have all
    989989                         specified points masked.
    990             weight:      Weighting scheme. 'none' (default), or
    991                          'var' (variance weighted)
     990            weight:      Weighting scheme. 'none' (default), 'var' (1/var(spec)
     991                         weighted), or 'tsys' (1/Tsys**2 weighted)
    992992            insitu:      if False a new scantable is returned.
    993993                         Otherwise, the scaling is done in-situ
  • trunk/src/SDMath.cc

    r519 r532  
    813813{
    814814   WeightType wtType = NONE;
    815    convertWeightString(wtType, weightStr, False);
    816 
    817    const uInt nRows = in.nRow();
     815   convertWeightString(wtType, weightStr, True);
    818816
    819817// Create output Table and reshape number of polarizations
     
    824822  header.npol = 1;
    825823  pTabOut->putSDHeader(header);
     824//
     825  const Table& tabIn = in.table();
    826826
    827827// Shape of input and output data
     
    834834  }
    835835//
     836  const uInt nRows = in.nRow();
    836837  const uInt nChan = shapeIn(asap::ChanAxis);
     838  AlwaysAssert(asap::nAxes==4,AipsError);
    837839  const IPosition vecShapeOut(4,1,1,1,nChan);     // A multi-dim form of a Vector shape
    838840  IPosition start(4), end(4);
     
    843845  Array<Bool> outMask(shapeOut, True);
    844846  const IPosition axes(2, asap::PolAxis, asap::ChanAxis);              // pol-channel plane
     847
     848// Attach Tsys column if needed
     849
     850  ROArrayColumn<Float> tSysCol;
     851  Array<Float> tSys;
     852  if (wtType==TSYS) {
     853     tSysCol.attach(tabIn,"TSYS");
     854  }
    845855//
    846856  const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis));
     
    855865      Array<Float>& arr = marr.getRWArray();
    856866      const Array<Bool>& barr = marr.getMask();
     867     
     868// Get Tsys
     869
     870      if (wtType==TSYS) {
     871         tSysCol.get(iRow,tSys);
     872      }
    857873
    858874// Make iterators to iterate by pol-channel planes
     875// The tSys array is empty unless wtType=TSYS so only
     876// access the iterator is that is the case
    859877
    860878      ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
    861879      ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
     880      ReadOnlyArrayIterator<Float>* pItTsysPlane = 0;
     881      if (wtType==TSYS) pItTsysPlane = new ReadOnlyArrayIterator<Float>(tSys, axes);
    862882
    863883// Accumulations
     
    874894        Vector<Float> t1(nChan); t1 = 0.0;
    875895        Vector<Bool> t2(nChan); t2 = True;
     896        Float tSys = 0.0;
    876897        MaskedArray<Float> vecSum(t1,t2);
    877898        Float norm = 0.0;
     
    879900           ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
    880901           ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
     902//
     903           ReadOnlyVectorIterator<Float>* pItTsysVec = 0;
     904           if (wtType==TSYS) {
     905              pItTsysVec = new ReadOnlyVectorIterator<Float>(pItTsysPlane->array(), 1);
     906           }             
     907//
    881908           while (!itDataVec.pastEnd()) {     
    882909
     
    885912              if (useMask) {
    886913                 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
    887                  if (wtType==VAR) fac = 1.0 / variance(spec);
     914                 if (wtType==VAR) {
     915                    fac = 1.0 / variance(spec);
     916                 } else if (wtType==TSYS) {
     917                    tSys = pItTsysVec->vector()[0];      // Drop pseudo channel dependency
     918                    fac = 1.0 / tSys / tSys;
     919                 }                   
    888920              } else {
    889921                 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
    890                  if (wtType==VAR) fac = 1.0 / variance(spec);
     922                 if (wtType==VAR) {
     923                    fac = 1.0 / variance(spec);
     924                 } else if (wtType==TSYS) {
     925                    tSys = pItTsysVec->vector()[0];      // Drop pseudo channel dependency
     926                    fac = 1.0 / tSys / tSys;
     927                 }
    891928              }
    892929
     
    901938              itDataVec.next();
    902939              itMaskVec.next();
     940              if (wtType==TSYS) pItTsysVec->next();
    903941           }
     942           
     943// Clean up
     944
     945           if (pItTsysVec) {
     946              delete pItTsysVec;
     947              pItTsysVec = 0;
     948           }           
    904949        }
    905950
     
    926971        itDataPlane.next();
    927972        itMaskPlane.next();
     973        if (wtType==TSYS) pItTsysPlane->next();
    928974      }
    929975
     
    935981      putDataInSDC(sc, outData, outMask);
    936982      pTabOut->putSDContainer(sc);
     983//
     984      if (wtType==TSYS) {
     985         delete pItTsysPlane;
     986         pItTsysPlane = 0;
     987      }
    937988   }
    938989
  • trunk/src/SDMath.h

    r518 r532  
    121121                            casa::Bool doAll, casa::uInt what, casa::Bool tSys) const;
    122122
    123 // Average polarizations
     123// Average polarizations.
    124124   SDMemTable* averagePol(const SDMemTable& in, const casa::Vector<casa::Bool>& mask,
    125125                          const casa::String& wtStr) const;
Note: See TracChangeset for help on using the changeset viewer.