source: trunk/src/RowAccumulator.cpp@ 1892

Last change on this file since 1892 was 1819, checked in by Kana Sugimoto, 14 years ago

New Development: No

JIRA Issue: No (merge alma branch to trunk)

Ready for Test: Yes

Interface Changes: No

Test Programs: regressions may work

Module(s): all single dish modules

Description:

Merged all changes in alma (r1386:1818) and newfiller (r1774:1818) branch.


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<uInt>(v.nelements(), 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
53 // only add these if not everything masked
54 if ( !allEQ(m, False) ) {
55 totalweight *= addTsys(tsys);
56 totalweight *= addInterval(interval);
57 addTime(time);
58 }
59 addSpectrum(v, m, totalweight);
60 initialized_ = True;
61}
62
63void RowAccumulator::addSpectrum( const Vector< Float >& v,
64 const Vector< Bool >& m,
65 Float weight)
66{
67 Float totalweight = weight;
68 MaskedArray<Float> data(v,m);
69 if ( weightType_ == asap::W_VAR ) {
70 if (m.nelements() == userMask_.nelements()) {
71 Float fac = 1.0/variance(data(userMask_));
72 totalweight *= fac;
73 } else {
74 Float fac = 1.0/variance(data);
75 totalweight *= fac;
76 }
77 }
78 data *= totalweight;
79 MaskedArray<Float> wadd(Vector<Float>(m.nelements(),totalweight), m);
80 weightSum_ += wadd;
81 spectrum_ += data;
82 const MaskedArray<uInt> inc(Vector<uInt>(m.nelements(),1), m);
83 n_ += inc;
84}
85
86Float RowAccumulator::addTsys( const casa::Vector< casa::Float > & v )
87{
88 // @fixme this assume tsys is the same for all channels
89
90 Float w = 1.0;
91 tsysSum_ += v[0];
92 if ( weightType_ == asap::W_TSYS || weightType_ == asap::W_TINTSYS ) {
93 w /= (v[0]*v[0]);
94 }
95 return w;
96}
97
98void asap::RowAccumulator::addTime( casa::Double t )
99{
100 timeSum_ += t;
101}
102
103Float asap::RowAccumulator::addInterval( casa::Double inter )
104{
105 Float w = 1.0;
106 intervalSum_ += inter;
107 if ( weightType_ == asap::W_TINT || weightType_ == asap::W_TINTSYS ) {
108 w /= Float(inter);
109 }
110 return w;
111}
112
113void asap::RowAccumulator::reset( )
114{
115 initialized_ = False;
116 intervalSum_ = 0.0;
117 tsysSum_.resize();
118 timeSum_ = 0.0;
119}
120
121casa::Vector< casa::Float > RowAccumulator::getSpectrum( ) const
122{
123 return (spectrum_/weightSum_).getArray();
124}
125
126casa::Double asap::RowAccumulator::getTime( ) const
127{
128 uInt n = max(n_);
129 return timeSum_/Float(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() > uInt(0));
141}
142
143casa::Vector< casa::Float > asap::RowAccumulator::getTsys( ) const
144{
145 // @fixme this assumes tsys.nelements() == 1
146 return tsysSum_/Float(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}
160
Note: See TracBrowser for help on using the repository browser.