Line data Source code
1 : //# VBGContinuumSubtractor.cc: Subtract the continuum from VisBuffGroups and
2 : //# write them to a different MS.
3 : //# Copyright (C) 2011
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 :
28 : #include <msvis/MSVis/VBGContinuumSubtractor.h>
29 : #include <msvis/MSVis/VBContinuumSubtractor.h>
30 : #include <msvis/MSVis/VBRemapper.h>
31 : #include <msvis/MSVis/SubMS.h>
32 : #include <msvis/MSVis/VisBufferComponents.h>
33 : #include <msvis/MSVis/VisBuffGroup.h>
34 : #include <msvis/MSVis/VisBuffGroupAcc.h>
35 : #include <casacore/casa/Exceptions/Error.h>
36 : #include <casacore/casa/Logging/LogIO.h>
37 : #include <casacore/ms/MSSel/MSSelection.h>
38 :
39 : using namespace casacore;
40 : namespace casa {
41 :
42 0 : VBGContinuumSubtractor::VBGContinuumSubtractor(MeasurementSet& outms,
43 : MSColumns *msc,
44 : const VBRemapper& remapper,
45 : const ROVisibilityIterator& invi,
46 : const uInt fitorder,
47 : const MS::PredefinedColumns datacol,
48 : const String& fitspw,
49 0 : const String& outspw) :
50 : GroupWriteToNewMS(outms, msc, remapper),
51 : fitorder_p(fitorder),
52 : datacol_p(datacol),
53 : outspw_p(outspw),
54 : rowsdone_p(0),
55 : tvi_debug(False),
56 0 : want_cont_p(False)
57 : {
58 0 : doWS_p = invi.existsWeightSpectrum();
59 0 : doFC_p = invi.existsFlagCategory();
60 :
61 : // Almost everything except the derived columns.
62 0 : prefetchColumns_p = asyncio::PrefetchColumns::prefetchColumns(
63 : VisBufferComponents::Ant1,
64 : VisBufferComponents::Ant2,
65 : VisBufferComponents::ArrayId,
66 : VisBufferComponents::DataDescriptionId,
67 : VisBufferComponents::Exposure,
68 : VisBufferComponents::Feed1,
69 : VisBufferComponents::Feed2,
70 : VisBufferComponents::FieldId,
71 : VisBufferComponents::FlagCube,
72 : VisBufferComponents::Flag,
73 : VisBufferComponents::FlagRow,
74 : VisBufferComponents::Freq,
75 : VisBufferComponents::ObservationId,
76 :
77 : // The cube always gets used, even if its
78 : // contents aren't.
79 : VisBufferComponents::ObservedCube,
80 :
81 : VisBufferComponents::NChannel,
82 : VisBufferComponents::NCorr,
83 : VisBufferComponents::NRow,
84 : VisBufferComponents::ProcessorId,
85 : VisBufferComponents::Scan,
86 : VisBufferComponents::SpW,
87 : VisBufferComponents::SigmaMat,
88 : VisBufferComponents::StateId,
89 : VisBufferComponents::Time,
90 : VisBufferComponents::TimeCentroid,
91 : VisBufferComponents::TimeInterval,
92 : VisBufferComponents::WeightMat,
93 : VisBufferComponents::UvwMat,
94 0 : -1);
95 0 : if(datacol == MS::MODEL_DATA)
96 0 : prefetchColumns_p.insert(VisBufferComponents::ModelCube);
97 0 : else if(datacol == MS::CORRECTED_DATA)
98 0 : prefetchColumns_p.insert(VisBufferComponents::CorrectedCube);
99 :
100 0 : if(doWS_p)
101 0 : prefetchColumns_p.insert(VisBufferComponents::WeightSpectrum);
102 0 : if(doFC_p)
103 0 : prefetchColumns_p.insert(VisBufferComponents::FlagCategory);
104 :
105 0 : VisBuffGroupAcc::fillChanMask(fitmask_p, fitspw, invi.ms());
106 :
107 0 : MSSelection mssel;
108 0 : mssel.setSpwExpr(outspw);
109 0 : Matrix<Int> chansel = mssel.getChanList(&invi.ms(), 1);
110 0 : Vector<Int> spws(chansel.column(0));
111 0 : uInt nselspws = spws.nelements();
112 :
113 0 : for(uInt i = 0; i < nselspws; ++i)
114 0 : outspws_p.insert(spws[i]);
115 0 : }
116 :
117 0 : VBGContinuumSubtractor::~VBGContinuumSubtractor()
118 : {
119 0 : VisBuffGroupAcc::clearChanMask(fitmask_p);
120 0 : }
121 :
122 : // VBGContinuumSubtractor& VBGContinuumSubtractor::operator=(const VBGContinuumSubtractor &other)
123 : // {
124 : // // trivial so far.
125 : // vi_p = other.vi_p;
126 : // return *this;
127 : // }
128 :
129 0 : Bool VBGContinuumSubtractor::process(VisBuffGroup& vbg)
130 : {
131 0 : Bool worked = true;
132 0 : uInt nvbs = vbg.nBuf();
133 0 : Int maxAnt = 0;
134 0 : Int maxSpw = 0; // VisBuffGroupAcc is 1 of those things that uses SpW when
135 : // it should use DDID.
136 0 : Int maxFld = 0;
137 :
138 : // Dagnabbit, VisBuffGroupAcc accumulates DATA and MODEL_DATA (even if it
139 : // isn't there, apparently), but not CORRECTED_DATA or FLOAT_DATA.
140 : // Compensate by moving the wanted column into DATA if necessary, before
141 : // accumulating.
142 0 : Bool otherToData = datacol_p != MS::DATA;
143 :
144 0 : for(uInt bufnum = 0; bufnum < nvbs; ++bufnum){
145 0 : if(otherToData)
146 0 : vbg(bufnum).visCube() = vbg(bufnum).dataCube(datacol_p);
147 :
148 0 : if(vbg(bufnum).numberAnt() > maxAnt) // Record maxAnt and maxFld
149 0 : maxAnt = vbg(bufnum).numberAnt(); // even for buffers that won't
150 0 : if(vbg(bufnum).fieldId() > maxFld) // be used in the fit.
151 0 : maxFld = vbg(bufnum).fieldId();
152 :
153 0 : Int spw = vbg(bufnum).spectralWindow();
154 0 : if(fitmask_p.count(spw) > 0){ // This requires fitspw to
155 : // follow the '' = nothing,
156 : // '*' = everything convention.
157 0 : if(spw > maxSpw)
158 0 : maxSpw = vbg(bufnum).spectralWindow();
159 : }
160 : }
161 :
162 : // Find the continuum
163 0 : VisBuffGroupAcc vbga(maxAnt + 1,
164 0 : maxSpw + 1, // VBContinuumSubtractor doesn't care
165 : // whether all the vbs have distinct spws or not.
166 0 : maxFld + 1, // There should only be 1 selected
167 : 0.0, // field, but its number is arbitrary.
168 0 : false); // VBGA is very Calibrater-centric when
169 : // it comes to MODEL_DATA.
170 :
171 0 : vbga.setTVIDebug(tvi_debug);
172 :
173 0 : for(uInt bufnum = 0; bufnum < nvbs; ++bufnum)
174 0 : if(fitmask_p.count(vbg(bufnum).spectralWindow()) > 0)
175 0 : vbga.accumulate(vbg(bufnum));
176 :
177 : // Select the fit channels after loading the VBs into vbga, so the VBs in vbg
178 : // are unaffected.
179 0 : vbga.applyChanMask(fitmask_p);
180 :
181 0 : vbga.finalizeAverage(); // Is this necessary when each vb is
182 : // being stored seperately?
183 :
184 0 : VBContinuumSubtractor vbcs;
185 :
186 : // It might be better later to cache the known lo and hi freqs, and use
187 : // vbcs.init(). See AMueller::selfSolveOne().
188 0 : vbcs.initFromVBGA(vbga);
189 :
190 0 : vbcs.setTVIDebug(tvi_debug);
191 :
192 : // datacol_p is in DATA now.
193 0 : vbcs.fit(vbga, fitorder_p, MS::DATA, coeffs_p, coeffsOK_p, false, true, false);
194 :
195 : //uInt oldrowsdone = rowsdone_p;
196 0 : for(uInt bufnum = 0; bufnum < nvbs; ++bufnum){
197 0 : uInt spw = vbg(bufnum).spectralWindow();
198 :
199 0 : if(outspws_p.find(spw) != outspws_p.end()){
200 : // datacol_p is in DATA now. Is this repetitious? Yes it is.
201 0 : if(!vbcs.apply(vbg(bufnum), MS::DATA, coeffs_p, coeffsOK_p, !want_cont_p,
202 0 : appliedSpWs_p.count(spw) < 1)){
203 0 : worked = false;
204 0 : break;
205 : }
206 0 : appliedSpWs_p.insert(spw);
207 :
208 : // Use SIGMA like a storage place for corrected weights.
209 0 : if(otherToData){
210 0 : vbg(bufnum).sigmaMat() = vbg(bufnum).weightMat();
211 0 : arrayTransformInPlace(vbg(bufnum).sigmaMat(), subms::wtToSigma);
212 : }
213 :
214 0 : rowsdone_p = GroupWriteToNewMS::write(outms_p, msc_p, vbg(bufnum),
215 0 : rowsdone_p, remapper_p,
216 0 : doFC_p,
217 : false, // for now
218 0 : doWS_p);
219 : //cerr << "Wrote out row IDs " << oldrowsdone << " - " << rowsdone_p - 1 << ",";
220 : }
221 : //else
222 : // cerr << "No output for";
223 : //cerr << " spw " << spw << endl;
224 : //oldrowsdone = rowsdone_p;
225 : }
226 :
227 0 : return worked;
228 : }
229 :
230 : using namespace casacore;
231 : } // end namespace casa
|
