source: trunk/src/RowAccumulator.cpp@ 1317

Last change on this file since 1317 was 1314, checked in by mar637, 18 years ago

fixed defect, where scantable averaging output failed if first row in ouput was all flagged

File size: 3.5 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 return timeSum_/max(n_);
128}
129
130casa::Double asap::RowAccumulator::getInterval( ) const
131{
132 return intervalSum_;
133}
134
135casa::Vector< casa::Bool > RowAccumulator::getMask( ) const
136{
137 return spectrum_.getMask();
138}
139
140casa::Vector< casa::Float > asap::RowAccumulator::getTsys( ) const
141{
142 // @fixme this assumes tsys.nelements() == 1
143 return tsysSum_/max(n_);
144}
145
146void asap::RowAccumulator::setUserMask( const casa::Vector< casa::Bool > & m )
147{
148 userMask_.resize();
149 userMask_ = m;
150}
Note: See TracBrowser for help on using the repository browser.