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: aips2-request@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 0 : StandardTsys::StandardTsys(VisSet& vs) :
60 : VisCal(vs), // virtual base
61 : VisMueller(vs), // virtual base
62 : BJones(vs), // immediate parent
63 : sysCalTabName_(vs.sysCalTableName()),
64 : freqDepCalWt_(false),
65 0 : freqDepTsys_(true)
66 : {
67 0 : if (prtlev()>2) cout << "StandardTsys::StandardTsys(vs)" << endl;
68 :
69 : // TBD: get these directly from the MS
70 0 : nChanParList()=vs.numberChan();
71 0 : startChanList()=vs.startChan(); // should be 0?
72 :
73 0 : }
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 : sysCalTabName_(""),
80 : 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 0 : StandardTsys::StandardTsys(const MSMetaInfoForCal& msmc) :
98 : VisCal(msmc), // virtual base
99 : VisMueller(msmc), // virtual base
100 : BJones(msmc), // immediate parent
101 : sysCalTabName_(""),
102 : freqDepCalWt_(false),
103 0 : freqDepTsys_(true)
104 : {
105 0 : if (prtlev()>2) cout << "StandardTsys::StandardTsys(msmc)" << endl;
106 :
107 : // Set the SYSCAL table name
108 0 : MeasurementSet ms(this->msmc().msname());
109 0 : sysCalTabName_ = ms.sysCalTableName();
110 :
111 : // OK?
112 0 : MSColumns mscol(ms);
113 0 : const MSSpWindowColumns& spwcols = mscol.spectralWindow();
114 0 : nChanParList()=spwcols.numChan().getColumn();
115 0 : startChanList().set(0);
116 :
117 0 : }
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 0 : StandardTsys::~StandardTsys() {
128 0 : if (prtlev()>2) cout << "StandardTsys::~StandardTsys()" << endl;
129 0 : }
130 :
131 0 : void StandardTsys::setSpecify(const Record& specify) {
132 :
133 0 : LogMessage message(LogOrigin("StandardTsys","setSpecify"));
134 :
135 : // Escape if SYSCAL table absent
136 0 : 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 0 : setSolved(false);
141 0 : setApplied(false);
142 :
143 : // Collect Cal table parameters
144 0 : if (specify.isDefined("caltable")) {
145 0 : calTableName()=specify.asString("caltable");
146 :
147 0 : 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 0 : logSink() << "Creating " << typeName()
155 : << " table from MS SYSCAL subtable."
156 0 : << LogIO::POST;
157 :
158 :
159 0 : Table sysCalTab(sysCalTabName_,Table::Old);
160 :
161 : // Verify required columns in SYSCAL
162 : {
163 0 : MSSysCal mssc(sysCalTab);
164 0 : MSSysCalColumns sscol(mssc);
165 0 : if ( (sscol.spectralWindowId().isNull() ||
166 0 : !sscol.spectralWindowId().isDefined(0)) ||
167 0 : (sscol.time().isNull() ||
168 0 : !sscol.time().isDefined(0)) ||
169 0 : (sscol.interval().isNull() ||
170 0 : !sscol.interval().isDefined(0)) ||
171 0 : (sscol.antennaId().isNull() ||
172 0 : !sscol.antennaId().isDefined(0)) )
173 0 : throw(AipsError("SYSCAL table is incomplete. Cannot proceed."));
174 :
175 0 : if ( !sscol.tsysSpectrum().isNull() && sscol.tsysSpectrum().isDefined(0) )
176 0 : freqDepTsys_ = True;
177 0 : else if ( !sscol.tsys().isNull() && sscol.tsys().isDefined(0) )
178 0 : freqDepTsys_ = False;
179 : else
180 0 : throw(AipsError("SYSCAL table has no Tsys measurements. Cannot proceed."));
181 : }
182 :
183 0 : if (!freqDepTsys_)
184 0 : nChanParList() = 1;
185 :
186 : // Create a new caltable to fill up
187 0 : createMemCalTable();
188 :
189 : // Setup shape of solveAllRPar
190 0 : initSolvePar();
191 :
192 0 : }
193 :
194 :
195 :
196 0 : void StandardTsys::specify(const Record& specify) {
197 :
198 : // Escape if SYSCAL table absent
199 0 : 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 0 : Matrix<Int> tsyscount(nSpw(),nAnt(),0);
204 0 : Matrix<Int> negTsys(nSpw(),nAnt(),0);
205 :
206 0 : Block<String> columns(2);
207 0 : columns[0] = "TIME";
208 0 : columns[1] = "SPECTRAL_WINDOW_ID";
209 0 : Table sysCalTab(sysCalTabName_,Table::Old);
210 0 : TableIterator sysCalIter(sysCalTab,columns);
211 :
212 : // Iterate
213 0 : while (!sysCalIter.pastEnd()) {
214 :
215 : // First extract info from SYSCAL
216 0 : MSSysCal mssc(sysCalIter.table());
217 0 : MSSysCalColumns sccol(mssc);
218 :
219 0 : Int ispw=sccol.spectralWindowId()(0);
220 0 : currSpw()=ispw; // registers everything else!
221 0 : Double timestamp=sccol.time()(0);
222 0 : Double interval=sccol.interval()(0);
223 :
224 0 : Vector<Int> ants;
225 0 : sccol.antennaId().getColumn(ants);
226 :
227 0 : Cube<Float> tsys;
228 0 : if (freqDepTsys_) {
229 0 : sccol.tsysSpectrum().getColumn(tsys);
230 : } else {
231 0 : Array<Float> tmp;
232 0 : sccol.tsys().getColumn(tmp);
233 0 : IPosition shape=tmp.shape();
234 0 : tsys=tmp.reform(IPosition(3, shape(0), 1, shape(1)));
235 : }
236 0 : IPosition tsysshape(tsys.shape());
237 :
238 : // Insist only that channel axis matches
239 0 : 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 0 : if (tsysshape(0)>2)
244 0 : throw(AipsError("Tsys pol axis is implausible"));
245 :
246 : // Now prepare to store in a caltable
247 0 : currSpw()=ispw; // registers everything else!
248 0 : refTime()=timestamp-interval/2.0;
249 0 : currField()=msmc().fieldIdAtTime(refTime());
250 0 : currScan()=msmc().scanNumberAtTime(refTime());
251 :
252 : // Initialize solveAllRPar, etc.
253 0 : solveAllRPar()=0.0;
254 0 : solveAllParOK()=true; // Assume all ok
255 0 : solveAllParErr()=0.1; // what should we use here? ~1/bandwidth?
256 0 : solveAllParSNR()=1.0;
257 :
258 0 : Int npol=tsysshape(0);
259 0 : IPosition blc(3,0,0,0), trc(3,npol-1,nChanPar()-1,0);
260 0 : for (uInt iant=0;iant<ants.nelements();++iant) {
261 0 : Int thisant=ants(iant);
262 0 : blc(2)=trc(2)=thisant; // the MS antenna index (not loop index)
263 0 : Array<Float> currtsys(tsys.xyPlane(iant));
264 0 : solveAllRPar()(blc,trc).nonDegenerate(2)=currtsys;
265 0 : solveAllParOK()(blc,trc)=true;
266 :
267 : // Increment tsys counter
268 0 : ++tsyscount(ispw,thisant);
269 :
270 0 : negTsys(ispw,thisant)+=nfalse(currtsys>=FLT_MIN);
271 :
272 : // Issue warnings for completely bogus Tsys spectra (per pol)
273 0 : for (Int ipol=0;ipol<npol;++ipol) {
274 0 : 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 0 : << " at t=" << MVTime(refTime()/C::day).string(MVTime::YMD,7)
279 : << " are all negative or zero will be entirely flagged."
280 0 : << LogIO::WARN << LogIO::POST;
281 : }
282 : }
283 : // Flag any Tsys<=0.0
284 :
285 0 : LogicalArray mask((solveAllRPar()>=FLT_MIN));
286 0 : MaskedArray<Bool> negs(solveAllParOK(),!mask);
287 0 : negs=false;
288 :
289 0 : Bool uniform=True;
290 0 : if (specify.isDefined("uniform"))
291 0 : uniform=specify.asBool("uniform");
292 :
293 0 : if (uniform)
294 0 : keepNCT();
295 : else
296 0 : keepNCT(ants);
297 :
298 0 : sysCalIter.next();
299 : }
300 :
301 : // Assign scan and fieldid info
302 : // 2016Aug23 (gmoellen): restore ad hoc method, since
303 : // Tsys timestamps are quirky
304 0 : assignCTScanField(*ct_,msName());
305 :
306 0 : logSink() << "Tsys spectra counts per spw for antenna Ids 0-"<<nElem()-1<<" (per pol):" << LogIO::POST;
307 0 : for (Int ispw=0;ispw<nSpw();++ispw) {
308 0 : Vector<Int> tsyscountspw(tsyscount.row(ispw));
309 0 : if (sum(tsyscountspw)>0) {
310 : logSink() << "Spw " << ispw << ": " << tsyscountspw
311 : << " (=" << sum(tsyscountspw) << " spectra;"
312 0 : << " " << nChanParList()(ispw) << " chans per spectra, per pol)"
313 0 : << LogIO::POST;
314 0 : for (Int iant=0;iant<nAnt();++iant) {
315 0 : if (negTsys(ispw,iant)>0)
316 0 : logSink() << " (Found and flagged " << negTsys(ispw,iant)
317 : << " spurious negative (or zero) Tsys channels for ant id="
318 : << iant << " in spw " << ispw << ".)"
319 0 : << LogIO::POST;
320 : }
321 : }
322 : // else
323 : // logSink() << "Spw " << ispw << ": NONE." << LogIO::POST;
324 : }
325 :
326 0 : }
327 :
328 0 : void StandardTsys::keepNCT(const Vector<Int>& ants) {
329 :
330 0 : for (uInt iant=0;iant<ants.nelements();++iant) {
331 0 : Int thisant=ants(iant);
332 :
333 0 : ct_->addRow(1);
334 :
335 0 : CTMainColumns ncmc(*ct_);
336 :
337 : // We are adding to the most-recently added row
338 0 : Int row=ct_->nrow()-1;
339 :
340 : // Meta-info
341 0 : ncmc.time().put(row,refTime());
342 0 : ncmc.fieldId().put(row,currField());
343 0 : ncmc.spwId().put(row,currSpw());
344 0 : ncmc.scanNo().put(row,currScan());
345 0 : ncmc.obsId().put(row,currObs());
346 0 : ncmc.interval().put(row,0.0);
347 0 : ncmc.antenna1().put(row,thisant);
348 0 : ncmc.antenna2().put(row,-1);
349 :
350 : // Params
351 0 : ncmc.fparam().put(row,solveAllRPar().xyPlane(thisant));
352 :
353 : // Stats
354 0 : ncmc.paramerr().put(row,solveAllParErr().xyPlane(thisant));
355 0 : ncmc.snr().put(row,solveAllParErr().xyPlane(thisant));
356 0 : ncmc.flag().put(row,!solveAllParOK().xyPlane(thisant));
357 : }
358 :
359 : // This spw now has some solutions in it
360 0 : spwOK_(currSpw())=True;
361 0 : }
362 :
363 0 : void StandardTsys::correct2(vi::VisBuffer2& vb, Bool trial,
364 : Bool doWtSp, Bool dosync) {
365 :
366 : // Signal channelized weight calibration downstream
367 0 : freqDepCalWt_=doWtSp;
368 :
369 : // Call parent:
370 0 : BJones::correct2(vb,trial,doWtSp,dosync);
371 0 : }
372 :
373 :
374 : // Specialized calcPar that does some sanity checking
375 0 : void StandardTsys::calcPar() {
376 :
377 : // Call parent (this drives generalized interpolation)
378 0 : 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 0 : Cube<Bool> mask(currRPar()<Float(0.0));
383 0 : currParOK()(mask)=false;
384 0 : currRPar()(mask)=0.0; // avoids NaN generation in sqrt, even for flagged points
385 0 : }
386 :
387 :
388 :
389 : // Specialized Jones Matrix calculation
390 0 : void StandardTsys::calcAllJones() {
391 :
392 : // Antenna-based factors are the sqrt(Tsys)
393 0 : convertArray(currJElem(),sqrt(currRPar()));
394 0 : currJElemOK()=currParOK();
395 :
396 0 : }
397 :
398 0 : void StandardTsys::syncWtScale() {
399 :
400 : // cout << "VJ::syncWtScale (" << typeName() << ")" << endl;
401 :
402 0 : 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 0 : switch (this->jonesType()) {
408 0 : case Jones::Scalar:
409 : case Jones::Diagonal: {
410 0 : currWtScale().resize(currRPar().shape()); // same as Tsys cube
411 0 : currWtScale()=0.0;
412 0 : 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 0 : calcWtScale2();
425 :
426 :
427 : // cout << "VJ::syncWtScale: currWtScale() = " << currWtScale() << endl;
428 :
429 :
430 : }
431 : else
432 :
433 : // Just do the single-chan thing
434 0 : BJones::syncWtScale();
435 :
436 : }
437 :
438 :
439 : // Calculate weight update
440 0 : void StandardTsys::calcWtScale() {
441 :
442 : // Initialize to 1.0
443 : // (will only be replaced if a valid calculation is possible)
444 0 : currWtScale().set(1.0);
445 :
446 0 : IPosition ash(currRPar().shape());
447 0 : ash(1)=1; // only on channel weight (so far)
448 0 : Cube<Float> cWS(currWtScale());
449 :
450 : // For each pol and antenna, form 1/mean(Tsys(f))
451 0 : IPosition it3(2,0,2);
452 0 : ArrayIterator<Float> Tsys(currRPar(),it3,false);
453 0 : ArrayIterator<Bool> Tok(currParOK(),it3,false);
454 0 : ArrayIterator<Float> cWSi(cWS,it3,false);
455 :
456 0 : while (!Tsys.pastEnd()) {
457 :
458 : // mask out flagged channels
459 0 : 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 0 : if (Tsysm.nelementsValid()>0) {
463 0 : Float meanTsys=mean(Tsysm);
464 0 : if (meanTsys>0.0)
465 0 : cWSi.array().set(1./meanTsys);
466 : }
467 :
468 0 : Tsys.next();
469 0 : Tok.next();
470 0 : cWSi.next();
471 : }
472 :
473 0 : }
474 :
475 : // Calculate weight update
476 0 : void StandardTsys::calcWtScale2() {
477 :
478 : // 1/Tsys (only where ok)
479 0 : currWtScale()(currParOK())=1.f/currRPar()(currParOK());
480 :
481 0 : }
482 :
483 : } //# NAMESPACE CASA - END
|
