source: branches/alma/src/RowAccumulator.cpp@ 1539

Last change on this file since 1539 was 1446, checked in by TakTsutsumi, 16 years ago

Merged recent updates (since 2007) from nrao-asap

File size: 3.7 KB
Line 
1//
2// C++ Implementation: RowAccumulator
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <Malte.Marquarding@csiro.au>, (C) 2005
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <casa/iomanip.h>
13
14#include <casa/Arrays/MaskArrMath.h>
15#include <casa/Arrays/MaskArrLogi.h>
16#include <casa/Arrays/ArrayMath.h>
17#include <casa/Arrays/ArrayLogical.h>
18#include "RowAccumulator.h"
19
20
21using namespace casa;
22using namespace asap;
23
24RowAccumulator::RowAccumulator(WeightType wt) :
25 weightType_(wt),
26 initialized_(False)
27{
28 reset();
29}
30
31RowAccumulator::~RowAccumulator()
32{
33}
34
35
36void RowAccumulator::add( const Vector< Float >& v,
37 const Vector< Bool >& m,
38 const Vector< Float >& tsys,
39 Double interval,
40 Double time )
41{
42 if (!initialized_) {
43 Vector<Float> dummy(v.nelements(), 0.0);
44 Vector<Bool> dummymsk(m.nelements(), True);
45 spectrum_.setData(dummy, dummymsk);
46 n_.setData(Vector<Float>(v.nelements(), 0.0), dummymsk);
47 weightSum_.setData(Vector<Float>(v.nelements(), 0.0), dummymsk);
48 tsysSum_.resize(tsys.nelements()); tsysSum_=0.0;
49 }
50 // add spectrum related weights, so far it is variance only.
51 Float totalweight = 1.0;
52 totalweight *= addTsys(tsys);
53 // only add these if not everything masked
54 if ( !allEQ(m, False) ) {
55 totalweight *= addInterval(interval);
56 addTime(time);
57 }
58 addSpectrum(v, m, totalweight);
59 initialized_ = True;
60}
61
62void RowAccumulator::addSpectrum( const Vector< Float >& v,
63 const Vector< Bool >& m,
64 Float weight)
65{
66 Float totalweight = weight;
67 MaskedArray<Float> data(v,m);
68 if ( weightType_ == asap::VAR ) {
69 if (m.nelements() == userMask_.nelements()) {
70 Float fac = 1.0/variance(data(userMask_));
71 totalweight *= fac;
72 } else {
73 Float fac = 1.0/variance(data);
74 totalweight *= fac;
75 }
76 }
77 data *= totalweight;
78 MaskedArray<Float> wadd(Vector<Float>(m.nelements(),totalweight), m);
79 weightSum_ += wadd;
80 spectrum_ += data;
81 const MaskedArray<Float> inc(Vector<Float>(m.nelements(),1.0), m);
82 n_ += inc;
83}
84
85Float RowAccumulator::addTsys( const casa::Vector< casa::Float > & v )
86{
87 // @fixme this assume tsys is the same for all channels
88
89 Float w = 1.0;
90 tsysSum_ += v[0];
91 if ( weightType_ == asap::TSYS || weightType_ == asap::TINTSYS ) {
92 w /= (v[0]*v[0]);
93 }
94 return w;
95}
96
97void asap::RowAccumulator::addTime( casa::Double t )
98{
99 timeSum_ += t;
100}
101
102Float asap::RowAccumulator::addInterval( casa::Double inter )
103{
104 Float w = 1.0;
105 intervalSum_ += inter;
106 if ( weightType_ == asap::TINT || weightType_ == asap::TINTSYS ) {
107 w /= Float(inter);
108 }
109 return w;
110}
111
112void asap::RowAccumulator::reset( )
113{
114 initialized_ = False;
115 intervalSum_ = 0.0;
116 tsysSum_.resize();
117 timeSum_ = 0.0;
118}
119
120casa::Vector< casa::Float > RowAccumulator::getSpectrum( ) const
121{
122 return (spectrum_/weightSum_).getArray();
123}
124
125casa::Double asap::RowAccumulator::getTime( ) const
126{
127 Float n = max(n_);
128 if (n < 1.0) n = 1.0;
129 return timeSum_/n;
130}
131
132casa::Double asap::RowAccumulator::getInterval( ) const
133{
134 return intervalSum_;
135}
136
137casa::Vector< casa::Bool > RowAccumulator::getMask( ) const
138{
139 // Return the "total" mask - False where no points have been accumulated.
140 return (n_.getArray() > Float(0.0));
141}
142
143casa::Vector< casa::Float > asap::RowAccumulator::getTsys( ) const
144{
145 // @fixme this assumes tsys.nelements() == 1
146 return tsysSum_/max(n_);
147}
148
149void asap::RowAccumulator::setUserMask( const casa::Vector< casa::Bool > & m )
150{
151 userMask_.resize();
152 userMask_ = m;
153}
154
155// Added by TT check the state of RowAccumulator
156casa::Bool RowAccumulator::state() const
157{
158 return initialized_;
159}
Note: See TracBrowser for help on using the repository browser.