Line data Source code
1 : //# TsysGainCal.cc: Implementation of Tsys calibration
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 :
27 : #include <synthesis/MeasurementComponents/TsysGainCal.h>
28 :
29 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
30 : // not yet: #include <synthesis/MeasurementComponents/CalCorruptor.h>
31 :
32 : #include <casacore/ms/MeasurementSets/MSColumns.h>
33 : #include <casacore/casa/BasicMath/Math.h>
34 : #include <casacore/tables/Tables/Table.h>
35 : #include <casacore/tables/Tables/TableIter.h>
36 : #include <synthesis/CalTables/CTGlobals.h>
37 :
38 : #include <casacore/casa/Arrays/MaskArrMath.h>
39 : #include <casacore/casa/BasicSL/String.h>
40 : #include <casacore/casa/Utilities/Assert.h>
41 : #include <casacore/casa/Utilities/GenSort.h>
42 : #include <casacore/casa/Exceptions/Error.h>
43 : #include <casacore/casa/System/Aipsrc.h>
44 :
45 : #include <sstream>
46 :
47 : #include <casacore/casa/Logging/LogMessage.h>
48 : #include <casacore/casa/Logging/LogSink.h>
49 : // math.h ?
50 :
51 : using namespace casacore;
52 : namespace casa { //# NAMESPACE CASA - BEGIN
53 :
54 :
55 : // **********************************************************
56 : // StandardTsys Implementations
57 : //
58 :
59 26 : StandardTsys::StandardTsys(VisSet& vs) :
60 : VisCal(vs), // virtual base
61 : VisMueller(vs), // virtual base
62 : BJones(vs), // immediate parent
63 26 : sysCalTabName_(vs.sysCalTableName()),
64 26 : freqDepCalWt_(false),
65 26 : freqDepTsys_(true)
66 : {
67 26 : if (prtlev()>2) cout << "StandardTsys::StandardTsys(vs)" << endl;
68 :
69 : // TBD: get these directly from the MS
70 26 : nChanParList()=vs.numberChan();
71 26 : startChanList()=vs.startChan(); // should be 0?
72 :
73 26 : }
74 :
75 0 : StandardTsys::StandardTsys(String msname,Int MSnAnt,Int MSnSpw) :
76 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
77 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
78 : BJones(msname,MSnAnt,MSnSpw), // immediate parent
79 0 : sysCalTabName_(""),
80 0 : freqDepCalWt_(false),
81 0 : freqDepTsys_(true)
82 : {
83 0 : if (prtlev()>2) cout << "StandardTsys::StandardTsys(msname,MSnAnt,MSnSpw)" << endl;
84 :
85 : // Set the SYSCAL table name
86 0 : MeasurementSet ms(msname);
87 0 : sysCalTabName_ = ms.sysCalTableName();
88 :
89 : // OK?
90 0 : MSColumns mscol(ms);
91 0 : const MSSpWindowColumns& spwcols = mscol.spectralWindow();
92 0 : nChanParList()=spwcols.numChan().getColumn();
93 0 : startChanList().set(0);
94 :
95 0 : }
96 :
97 29 : StandardTsys::StandardTsys(const MSMetaInfoForCal& msmc) :
98 : VisCal(msmc), // virtual base
99 : VisMueller(msmc), // virtual base
100 : BJones(msmc), // immediate parent
101 29 : sysCalTabName_(""),
102 29 : freqDepCalWt_(false),
103 29 : freqDepTsys_(true)
104 : {
105 29 : if (prtlev()>2) cout << "StandardTsys::StandardTsys(msmc)" << endl;
106 :
107 : // Set the SYSCAL table name
108 29 : MeasurementSet ms(this->msmc().msname());
109 29 : sysCalTabName_ = ms.sysCalTableName();
110 :
111 : // OK?
112 29 : MSColumns mscol(ms);
113 29 : const MSSpWindowColumns& spwcols = mscol.spectralWindow();
114 29 : nChanParList()=spwcols.numChan().getColumn();
115 29 : startChanList().set(0);
116 :
117 29 : }
118 :
119 0 : StandardTsys::StandardTsys(const Int& nAnt) :
120 : VisCal(nAnt),
121 : VisMueller(nAnt),
122 0 : BJones(nAnt)
123 : {
124 0 : if (prtlev()>2) cout << "StandardTsys::StandardTsys(nAnt)" << endl;
125 0 : }
126 :
127 86 : StandardTsys::~StandardTsys() {
128 55 : if (prtlev()>2) cout << "StandardTsys::~StandardTsys()" << endl;
129 86 : }
130 :
131 18 : void StandardTsys::setSpecify(const Record& specify) {
132 :
133 36 : LogMessage message(LogOrigin("StandardTsys","setSpecify"));
134 :
135 : // Escape if SYSCAL table absent
136 18 : if (!Table::isReadable(sysCalTabName_))
137 0 : throw(AipsError("The SYSCAL subtable is not present in the specified MS."));
138 :
139 : // Not actually applying or solving
140 18 : setSolved(false);
141 18 : setApplied(false);
142 :
143 : // Collect Cal table parameters
144 18 : if (specify.isDefined("caltable")) {
145 18 : calTableName()=specify.asString("caltable");
146 :
147 18 : if (Table::isReadable(calTableName()))
148 : logSink() << "FYI: We are going to overwrite an existing CalTable: "
149 0 : << calTableName()
150 0 : << LogIO::POST;
151 : }
152 :
153 : // we are creating a table from scratch
154 18 : logSink() << "Creating " << typeName()
155 : << " table from MS SYSCAL subtable."
156 18 : << LogIO::POST;
157 :
158 :
159 18 : Table sysCalTab(sysCalTabName_,Table::Old);
160 :
161 : // Verify required columns in SYSCAL
162 : {
163 18 : MSSysCal mssc(sysCalTab);
164 18 : MSSysCalColumns sscol(mssc);
165 18 : if ( (sscol.spectralWindowId().isNull() ||
166 36 : !sscol.spectralWindowId().isDefined(0)) ||
167 18 : (sscol.time().isNull() ||
168 36 : !sscol.time().isDefined(0)) ||
169 18 : (sscol.interval().isNull() ||
170 72 : !sscol.interval().isDefined(0)) ||
171 18 : (sscol.antennaId().isNull() ||
172 18 : !sscol.antennaId().isDefined(0)) )
173 0 : throw(AipsError("SYSCAL table is incomplete. Cannot proceed."));
174 :
175 18 : if ( !sscol.tsysSpectrum().isNull() && sscol.tsysSpectrum().isDefined(0) )
176 16 : freqDepTsys_ = True;
177 2 : else if ( !sscol.tsys().isNull() && sscol.tsys().isDefined(0) )
178 2 : freqDepTsys_ = False;
179 : else
180 0 : throw(AipsError("SYSCAL table has no Tsys measurements. Cannot proceed."));
181 18 : }
182 :
183 18 : if (!freqDepTsys_)
184 2 : nChanParList() = 1;
185 :
186 : // Create a new caltable to fill up
187 18 : createMemCalTable();
188 :
189 : // Setup shape of solveAllRPar
190 18 : initSolvePar();
191 :
192 18 : }
193 :
194 :
195 :
196 18 : void StandardTsys::specify(const Record& specify) {
197 :
198 : // Escape if SYSCAL table absent
199 18 : if (!Table::isReadable(sysCalTabName_))
200 0 : throw(AipsError("The SYSCAL subtable is not present in the specified MS. Tsys unavailable."));
201 :
202 : // Keep a count of the number of Tsys found per spw/ant
203 18 : Matrix<Int> tsyscount(nSpw(),nAnt(),0);
204 18 : Matrix<Int> negTsys(nSpw(),nAnt(),0);
205 :
206 18 : Block<String> columns(2);
207 18 : columns[0] = "TIME";
208 18 : columns[1] = "SPECTRAL_WINDOW_ID";
209 18 : Table sysCalTab(sysCalTabName_,Table::Old);
210 18 : TableIterator sysCalIter(sysCalTab,columns);
211 :
212 : // Iterate
213 29619 : while (!sysCalIter.pastEnd()) {
214 :
215 : // First extract info from SYSCAL
216 29601 : MSSysCal mssc(sysCalIter.table());
217 29601 : MSSysCalColumns sccol(mssc);
218 :
219 29601 : Int ispw=sccol.spectralWindowId()(0);
220 29601 : currSpw()=ispw; // registers everything else!
221 29601 : Double timestamp=sccol.time()(0);
222 29601 : Double interval=sccol.interval()(0);
223 :
224 29601 : Vector<Int> ants;
225 29601 : sccol.antennaId().getColumn(ants);
226 :
227 29601 : Cube<Float> tsys;
228 29601 : if (freqDepTsys_) {
229 105 : sccol.tsysSpectrum().getColumn(tsys);
230 : } else {
231 29496 : Array<Float> tmp;
232 29496 : sccol.tsys().getColumn(tmp);
233 29496 : IPosition shape=tmp.shape();
234 29496 : tsys=tmp.reform(IPosition(3, shape(0), 1, shape(1)));
235 29496 : }
236 29601 : IPosition tsysshape(tsys.shape());
237 :
238 : // Insist only that channel axis matches
239 29601 : if (tsysshape(1)!=nChanPar())
240 0 : throw(AipsError("SYSCAL Tsys Spectrum channel axis shape doesn't match data! Cannot proceed."));
241 :
242 : // ...and that tsys pol axis makes sense
243 29601 : if (tsysshape(0)>2)
244 0 : throw(AipsError("Tsys pol axis is implausible"));
245 :
246 : // Now prepare to store in a caltable
247 29601 : currSpw()=ispw; // registers everything else!
248 29601 : refTime()=timestamp-interval/2.0;
249 29601 : currField()=msmc().fieldIdAtTime(refTime());
250 29601 : currScan()=msmc().scanNumberAtTime(refTime());
251 :
252 : // Initialize solveAllRPar, etc.
253 29601 : solveAllRPar()=0.0;
254 29601 : solveAllParOK()=true; // Assume all ok
255 29601 : solveAllParErr()=0.1; // what should we use here? ~1/bandwidth?
256 29601 : solveAllParSNR()=1.0;
257 :
258 29601 : Int npol=tsysshape(0);
259 29601 : IPosition blc(3,0,0,0), trc(3,npol-1,nChanPar()-1,0);
260 78482 : for (uInt iant=0;iant<ants.nelements();++iant) {
261 48881 : Int thisant=ants(iant);
262 48881 : blc(2)=trc(2)=thisant; // the MS antenna index (not loop index)
263 48881 : Array<Float> currtsys(tsys.xyPlane(iant));
264 48881 : solveAllRPar()(blc,trc).nonDegenerate(2)=currtsys;
265 48881 : solveAllParOK()(blc,trc)=true;
266 :
267 : // Increment tsys counter
268 48881 : ++tsyscount(ispw,thisant);
269 :
270 48881 : negTsys(ispw,thisant)+=nfalse(currtsys>=FLT_MIN);
271 :
272 : // Issue warnings for completely bogus Tsys spectra (per pol)
273 146643 : for (Int ipol=0;ipol<npol;++ipol) {
274 97762 : if (ntrue(Matrix<Float>(currtsys).row(ipol)>=FLT_MIN)==0)
275 : logSink() << " Tsys data for ant id="
276 : << thisant << " (pol=" << ipol<< ")"
277 : << " in spw " << ispw
278 96 : << " at t=" << MVTime(refTime()/C::day).string(MVTime::YMD,7)
279 : << " are all negative or zero will be entirely flagged."
280 96 : << LogIO::WARN << LogIO::POST;
281 : }
282 48881 : }
283 : // Flag any Tsys<=0.0
284 :
285 29601 : LogicalArray mask((solveAllRPar()>=FLT_MIN));
286 29601 : MaskedArray<Bool> negs(solveAllParOK(),!mask);
287 29601 : negs=false;
288 :
289 29601 : Bool uniform=True;
290 29601 : if (specify.isDefined("uniform"))
291 29601 : uniform=specify.asBool("uniform");
292 :
293 29601 : if (uniform)
294 105 : keepNCT();
295 : else
296 29496 : keepNCT(ants);
297 :
298 29601 : sysCalIter.next();
299 29601 : }
300 :
301 : // Assign scan and fieldid info
302 : // 2016Aug23 (gmoellen): restore ad hoc method, since
303 : // Tsys timestamps are quirky
304 18 : assignCTScanField(*ct_,msName());
305 :
306 18 : logSink() << "Tsys spectra counts per spw for antenna Ids 0-"<<nElem()-1<<" (per pol):" << LogIO::POST;
307 290 : for (Int ispw=0;ispw<nSpw();++ispw) {
308 272 : Vector<Int> tsyscountspw(tsyscount.row(ispw));
309 272 : if (sum(tsyscountspw)>0) {
310 : logSink() << "Spw " << ispw << ": " << tsyscountspw
311 : << " (=" << sum(tsyscountspw) << " spectra;"
312 69 : << " " << nChanParList()(ispw) << " chans per spectra, per pol)"
313 138 : << LogIO::POST;
314 350 : for (Int iant=0;iant<nAnt();++iant) {
315 281 : if (negTsys(ispw,iant)>0)
316 20 : logSink() << " (Found and flagged " << negTsys(ispw,iant)
317 : << " spurious negative (or zero) Tsys channels for ant id="
318 : << iant << " in spw " << ispw << ".)"
319 10 : << LogIO::POST;
320 : }
321 : }
322 : // else
323 : // logSink() << "Spw " << ispw << ": NONE." << LogIO::POST;
324 272 : }
325 :
326 18 : }
327 :
328 29496 : void StandardTsys::keepNCT(const Vector<Int>& ants) {
329 :
330 78008 : for (uInt iant=0;iant<ants.nelements();++iant) {
331 48512 : Int thisant=ants(iant);
332 :
333 48512 : ct_->addRow(1);
334 :
335 48512 : CTMainColumns ncmc(*ct_);
336 :
337 : // We are adding to the most-recently added row
338 48512 : Int row=ct_->nrow()-1;
339 :
340 : // Meta-info
341 48512 : ncmc.time().put(row,refTime());
342 48512 : ncmc.fieldId().put(row,currField());
343 48512 : ncmc.spwId().put(row,currSpw());
344 48512 : ncmc.scanNo().put(row,currScan());
345 48512 : ncmc.obsId().put(row,currObs());
346 48512 : ncmc.interval().put(row,0.0);
347 48512 : ncmc.antenna1().put(row,thisant);
348 48512 : ncmc.antenna2().put(row,-1);
349 :
350 : // Params
351 48512 : ncmc.fparam().put(row,solveAllRPar().xyPlane(thisant));
352 :
353 : // Stats
354 48512 : ncmc.paramerr().put(row,solveAllParErr().xyPlane(thisant));
355 48512 : ncmc.snr().put(row,solveAllParErr().xyPlane(thisant));
356 48512 : ncmc.flag().put(row,!solveAllParOK().xyPlane(thisant));
357 48512 : }
358 :
359 : // This spw now has some solutions in it
360 29496 : spwOK_(currSpw())=True;
361 29496 : }
362 :
363 2596 : void StandardTsys::correct2(vi::VisBuffer2& vb, Bool trial,
364 : Bool doWtSp, Bool dosync) {
365 :
366 : // Signal channelized weight calibration downstream
367 2596 : freqDepCalWt_=doWtSp;
368 :
369 : // Call parent:
370 2596 : BJones::correct2(vb,trial,doWtSp,dosync);
371 2596 : }
372 :
373 :
374 : // Specialized calcPar that does some sanity checking
375 2596 : void StandardTsys::calcPar() {
376 :
377 : // Call parent (this drives generalized interpolation)
378 2596 : SolvableVisCal::calcPar();
379 :
380 : // Since some interpolation types may unwittingly yield
381 : // negative Tsys, we'll trap that here by flagging and zeroing them
382 2596 : Cube<Bool> mask(currRPar()<Float(0.0));
383 2596 : currParOK()(mask)=false;
384 2596 : currRPar()(mask)=0.0; // avoids NaN generation in sqrt, even for flagged points
385 2596 : }
386 :
387 :
388 :
389 : // Specialized Jones Matrix calculation
390 997 : void StandardTsys::calcAllJones() {
391 :
392 : // Antenna-based factors are the sqrt(Tsys)
393 997 : convertArray(currJElem(),sqrt(currRPar()));
394 997 : currJElemOK()=currParOK();
395 :
396 997 : }
397 :
398 385 : void StandardTsys::syncWtScale() {
399 :
400 : // cout << "VJ::syncWtScale (" << typeName() << ")" << endl;
401 :
402 385 : if (freqDepCalWt_) {
403 :
404 : // We _will_ do a channelized weight calibration, so shape currWtScale() appropriately
405 :
406 : // Ensure proper size according to Jones matrix type
407 180 : switch (this->jonesType()) {
408 180 : case Jones::Scalar:
409 : case Jones::Diagonal: {
410 180 : currWtScale().resize(currRPar().shape()); // same as Tsys cube
411 180 : currWtScale()=0.0;
412 180 : break;
413 : }
414 0 : default: {
415 : // Only diag and scalar versions can adjust weights
416 : // cout<< "Turning off calWt()" << endl;
417 0 : calWt()=false;
418 0 : return;
419 : break;
420 : }
421 : }
422 :
423 : // Calculate the weight scale factors spectrally
424 180 : calcWtScale2();
425 :
426 :
427 : // cout << "VJ::syncWtScale: currWtScale() = " << currWtScale() << endl;
428 :
429 :
430 : }
431 : else
432 :
433 : // Just do the single-chan thing
434 205 : BJones::syncWtScale();
435 :
436 : }
437 :
438 :
439 : // Calculate weight update
440 205 : void StandardTsys::calcWtScale() {
441 :
442 : // Initialize to 1.0
443 : // (will only be replaced if a valid calculation is possible)
444 205 : currWtScale().set(1.0);
445 :
446 205 : IPosition ash(currRPar().shape());
447 205 : ash(1)=1; // only on channel weight (so far)
448 205 : Cube<Float> cWS(currWtScale());
449 :
450 : // For each pol and antenna, form 1/mean(Tsys(f))
451 205 : IPosition it3(2,0,2);
452 205 : ArrayIterator<Float> Tsys(currRPar(),it3,false);
453 205 : ArrayIterator<Bool> Tok(currParOK(),it3,false);
454 205 : ArrayIterator<Float> cWSi(cWS,it3,false);
455 :
456 1263 : while (!Tsys.pastEnd()) {
457 :
458 : // mask out flagged channels
459 1058 : MaskedArray<Float> Tsysm(Tsys.array()(Tok.array()));
460 :
461 : // If any good Tsys this ant/pol, calc the wt scale (else it remains 1.0)
462 1058 : if (Tsysm.nelementsValid()>0) {
463 1058 : Float meanTsys=mean(Tsysm);
464 1058 : if (meanTsys>0.0)
465 1058 : cWSi.array().set(1./meanTsys);
466 : }
467 :
468 1058 : Tsys.next();
469 1058 : Tok.next();
470 1058 : cWSi.next();
471 1058 : }
472 :
473 205 : }
474 :
475 : // Calculate weight update
476 180 : void StandardTsys::calcWtScale2() {
477 :
478 : // 1/Tsys (only where ok)
479 180 : currWtScale()(currParOK())=1.f/currRPar()(currParOK());
480 :
481 180 : }
482 :
483 : } //# NAMESPACE CASA - END
|
