Line data Source code
1 : /*
2 : * SDDoubleCircleGainCal.cpp
3 : *
4 : * Created on: Jun 3, 2016
5 : * Author: nakazato
6 : */
7 :
8 : #include <synthesis/MeasurementComponents/SDDoubleCircleGainCal.h>
9 : #include <synthesis/MeasurementComponents/StandardVisCal.h>
10 : #include <synthesis/MeasurementComponents/SolvableVisCal.h>
11 : #include <synthesis/MeasurementComponents/SDDoubleCircleGainCalImpl.h>
12 : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
13 : #include <synthesis/MeasurementEquations/VisEquation.h>
14 : #include <synthesis/Utilities/PointingDirectionCalculator.h>
15 : #include <synthesis/Utilities/PointingDirectionProjector.h>
16 : #include <synthesis/CalTables/CTIter.h>
17 : #include <msvis/MSVis/VisSet.h>
18 :
19 : #include <casacore/casa/BasicSL/String.h>
20 : #include <casacore/casa/Arrays/Vector.h>
21 : #include <casacore/casa/IO/ArrayIO.h>
22 : #include <casacore/casa/Quanta/Quantum.h>
23 : #include <casacore/casa/Quanta/QuantumHolder.h>
24 : #include <casacore/casa/Utilities/Assert.h>
25 : #include <casacore/tables/TaQL/TableParse.h>
26 : #include <casacore/measures/TableMeasures/ScalarQuantColumn.h>
27 : #include <casacore/measures/TableMeasures/ArrayQuantColumn.h>
28 : #include <casacore/ms/MeasurementSets/MSIter.h>
29 :
30 : #include <iostream>
31 : #include <sstream>
32 : #include <cassert>
33 : #include <map>
34 :
35 : // Debug Message Handling
36 : // if SDGAIN_DEBUG is defined, the macro debuglog and
37 : // debugpost point standard output stream (std::cout and
38 : // std::endl so that debug messages are sent to standard
39 : // output. Otherwise, these macros basically does nothing.
40 : // "Do nothing" behavior is implemented in NullLogger
41 : // and its associating << operator below.
42 : //
43 : // Usage:
44 : // Similar to standard output stream.
45 : //
46 : // debuglog << "Any message" << any_value << debugpost;
47 : //
48 : //#define SDGAINDBLC_DEBUG
49 :
50 : namespace {
51 : struct NullLogger {
52 : };
53 :
54 : template<class T>
55 50064 : inline NullLogger &operator<<(NullLogger &logger, T /*value*/) {
56 50064 : return logger;
57 : }
58 : }
59 :
60 : #ifdef SDGAINDBLC_DEBUG
61 : #define debuglog std::cerr
62 : #define debugpost std::endl
63 : #else
64 : ::NullLogger nulllogger;
65 : #define debuglog nulllogger
66 : #define debugpost 0
67 : #endif
68 :
69 : namespace {
70 0 : inline void fillNChanParList(casacore::String const &msName,
71 : casacore::Vector<casacore::Int> &nChanParList) {
72 0 : casacore::MeasurementSet const ms(msName);
73 0 : casacore::MSSpectralWindow const &msspw = ms.spectralWindow();
74 0 : casacore::ScalarColumn<casacore::Int> nchanCol(msspw, "NUM_CHAN");
75 0 : debuglog << "nchanCol=" << nchanCol.getColumn() << debugpost;
76 0 : nChanParList = nchanCol.getColumn()(
77 0 : casacore::Slice(0, nChanParList.nelements()));
78 0 : }
79 :
80 : template<class T>
81 0 : inline casacore::String toString(casacore::Vector<T> const &v) {
82 0 : std::ostringstream oss;
83 0 : oss << "[";
84 0 : std::string delimiter = "";
85 0 : for (size_t i = 0; i < v.nelements(); ++i) {
86 0 : oss << delimiter << v[i];
87 0 : delimiter = ",";
88 : }
89 0 : oss << "]";
90 0 : return casacore::String(oss);
91 : }
92 :
93 0 : inline casacore::String selectOnSourceAutoCorr(
94 : casacore::MeasurementSet const &ms) {
95 0 : debuglog << "selectOnSource" << debugpost;
96 : casacore::String taqlForState(
97 0 : "SELECT FLAG_ROW FROM $1 WHERE UPPER(OBS_MODE) ~ m/^OBSERVE_TARGET#ON_SOURCE/");
98 0 : casacore::Table stateSel = casacore::tableCommand(taqlForState, ms.state()).table();
99 0 : auto stateIdList = stateSel.rowNumbers();
100 0 : debuglog << "stateIdList = " << stateIdList << debugpost;
101 0 : std::ostringstream oss;
102 0 : oss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && STATE_ID IN "
103 0 : << toString(stateIdList)
104 0 : << " ORDER BY FIELD_ID, ANTENNA1, FEED1, DATA_DESC_ID, TIME";
105 0 : return casacore::String(oss);
106 : }
107 :
108 : class DataColumnAccessor {
109 : public:
110 0 : DataColumnAccessor(casacore::Table const &ms,
111 0 : casacore::String const colName = "DATA") :
112 0 : dataCol_(ms, colName) {
113 0 : }
114 : casacore::Matrix<casacore::Float> operator()(size_t irow) {
115 : return casacore::real(dataCol_(irow));
116 : }
117 0 : casacore::Cube<casacore::Float> getColumn() {
118 0 : return casacore::real(dataCol_.getColumn());
119 : }
120 : private:
121 : DataColumnAccessor() {
122 : }
123 : casacore::ArrayColumn<casacore::Complex> dataCol_;
124 : };
125 :
126 : class FloatDataColumnAccessor {
127 : public:
128 0 : FloatDataColumnAccessor(casacore::Table const &ms) :
129 0 : dataCol_(ms, "FLOAT_DATA") {
130 0 : }
131 : casacore::Matrix<casacore::Float> operator()(size_t irow) {
132 : return dataCol_(irow);
133 : }
134 0 : casacore::Cube<casacore::Float> getColumn() {
135 0 : return dataCol_.getColumn();
136 : }
137 : private:
138 : FloatDataColumnAccessor() {
139 : }
140 : casacore::ArrayColumn<casacore::Float> dataCol_;
141 : };
142 :
143 0 : inline bool isEphemeris(casacore::String const &name) {
144 : // Check if given name is included in MDirection types
145 : casacore::Int nall, nextra;
146 : const casacore::uInt *typ;
147 0 : auto *types = casacore::MDirection::allMyTypes(nall, nextra, typ);
148 0 : auto start_extra = nall - nextra;
149 0 : auto capital_name = name;
150 0 : capital_name.upcase();
151 :
152 0 : for (auto i = start_extra; i < nall; ++i) {
153 0 : if (capital_name == types[i]) {
154 0 : return true;
155 : }
156 : }
157 :
158 0 : return false;
159 : }
160 :
161 0 : inline void updateWeight(casa::NewCalTable &ct) {
162 0 : casa::CTMainColumns ctmc(ct);
163 :
164 : // simply copy FPARAM
165 0 : for (size_t irow = 0; irow < ct.nrow(); ++irow) {
166 0 : ctmc.weight().put(irow, real(ctmc.cparam()(irow)));
167 : }
168 0 : }
169 :
170 28 : casacore::Double rad2arcsec(casacore::Double value_in_rad) {
171 28 : return casacore::Quantity(value_in_rad, "rad").getValue("arcsec");
172 : }
173 : }
174 :
175 : using namespace casacore;
176 : namespace casa {
177 0 : SDDoubleCircleGainCal::SDDoubleCircleGainCal(VisSet &vs) :
178 : VisCal(vs), // virtual base
179 : VisMueller(vs), // virtual base
180 0 : GJones(vs), central_disk_size_(0.0), smooth_(True), currAnt_() {
181 0 : }
182 :
183 10 : SDDoubleCircleGainCal::SDDoubleCircleGainCal(const MSMetaInfoForCal& msmc) :
184 : VisCal(msmc), // virtual base
185 : VisMueller(msmc), // virtual base
186 10 : GJones(msmc), central_disk_size_(0.0), smooth_(True), currAnt_() {
187 10 : }
188 :
189 20 : SDDoubleCircleGainCal::~SDDoubleCircleGainCal() {
190 20 : }
191 :
192 0 : void SDDoubleCircleGainCal::setSolve() {
193 0 : central_disk_size_ = 0.0;
194 0 : smooth_ = True;
195 :
196 : // call parent setSolve
197 0 : SolvableVisCal::setSolve();
198 :
199 : // solint forcibly set to 'int'
200 0 : solint() = "int";
201 0 : }
202 :
203 10 : void SDDoubleCircleGainCal::setSolve(const Record &solve) {
204 : // parameters for double circle gain calibration
205 10 : if (solve.isDefined("smooth")) {
206 10 : smooth_ = solve.asBool("smooth");
207 : }
208 10 : if (solve.isDefined("radius")) {
209 20 : String size_string = solve.asString("radius");
210 20 : QuantumHolder qh;
211 20 : String error;
212 10 : qh.fromString(error, size_string);
213 10 : Quantity q = qh.asQuantity();
214 10 : central_disk_size_ = q.getValue("rad");
215 : }
216 :
217 10 : logSink() << LogIO::DEBUGGING << "smooth=" << smooth_ << LogIO::POST;
218 : logSink() << LogIO::DEBUGGING << "central disk size=" << central_disk_size_
219 : << " rad" << "(" << rad2arcsec(central_disk_size_) << " arcsec)"
220 10 : << LogIO::POST;
221 :
222 10 : if (central_disk_size_ < 0.0) {
223 1 : logSink() << "Negative central disk size is given" << LogIO::EXCEPTION;
224 : }
225 :
226 : // call parent setSolve
227 9 : SolvableVisCal::setSolve(solve);
228 :
229 : // solint forcibly set to 'int'
230 9 : solint() = "int";
231 9 : }
232 :
233 18 : String SDDoubleCircleGainCal::solveinfo() {
234 36 : ostringstream o;
235 36 : o << typeName() << ": " << calTableName() << " smooth="
236 18 : << (smooth_ ? "True" : "False") << endl << " radius="
237 18 : << central_disk_size_;
238 18 : if (central_disk_size_ == 0.0) {
239 2 : o << " (half of primary beam will be used)";
240 : }
241 18 : o << endl;
242 36 : return String(o);
243 : }
244 :
245 9 : void SDDoubleCircleGainCal::globalPostSolveTinker() {
246 :
247 : // apply generic post-solve stuff - is it necessary?
248 9 : SolvableVisCal::globalPostSolveTinker();
249 :
250 : // double-circle gain calibration is implemented here
251 : // assuming that given caltable (on memory) has complete
252 : // set of spectral data required for the calibration
253 :
254 : // setup worker_
255 9 : ROScalarQuantColumn<Double> antennaDiameterColumn(ct_->antenna(),
256 27 : "DISH_DIAMETER");
257 9 : ROArrayQuantColumn<Double> observingFrequencyColumn(ct_->spectralWindow(),
258 27 : "CHAN_FREQ");
259 9 : Int smoothingSize = -1;// use default smoothing size
260 9 : worker_.setCentralRegion(central_disk_size_);
261 9 : if (smooth_) {
262 9 : worker_.setSmoothing(smoothingSize);
263 : } else {
264 0 : worker_.unsetSmoothing();
265 : }
266 :
267 : // sort caltable by TIME
268 27 : NewCalTable sorted(ct_->sort("TIME"));
269 18 : Block<String> col(3);
270 9 : col[0] = "SPECTRAL_WINDOW_ID";
271 9 : col[1] = "FIELD_ID";
272 9 : col[2] = "ANTENNA1";
273 : //col[3] = "FEED1";
274 18 : CTIter ctiter(sorted,col);
275 :
276 18 : Vector<rownr_t> to_be_removed;
277 27 : while (!ctiter.pastEnd()) {
278 18 : Int const thisAntenna = ctiter.thisAntenna1();
279 18 : Quantity antennaDiameterQuant = antennaDiameterColumn(thisAntenna); // nominal
280 18 : worker_.setAntennaDiameter(antennaDiameterQuant.getValue("m"));
281 18 : debuglog<< "antenna diameter = " << worker_.getAntennaDiameter() << "m" << debugpost;
282 18 : Int const thisSpw = ctiter.thisSpw();
283 18 : Vector<Quantity> observingFrequencyQuant = observingFrequencyColumn(thisSpw);
284 18 : Double meanFrequency = 0.0;
285 18 : auto numChan = observingFrequencyQuant.nelements();
286 18 : debuglog<< "numChan = " << numChan << debugpost;
287 18 : assert(numChan > 0);
288 18 : if (numChan % 2 == 0) {
289 0 : meanFrequency = (observingFrequencyQuant[numChan / 2 - 1].getValue("Hz")
290 0 : + observingFrequencyQuant[numChan / 2].getValue("Hz")) / 2.0;
291 : } else {
292 18 : meanFrequency = observingFrequencyQuant[numChan / 2].getValue("Hz");
293 : }
294 : //debuglog << "mean frequency " << meanFrequency.getValue() << " [" << meanFrequency.getFullUnit() << "]" << debugpost;
295 18 : debuglog<< "mean frequency " << meanFrequency << debugpost;
296 18 : worker_.setObservingFrequency(meanFrequency);
297 18 : debuglog<< "observing frequency = " << worker_.getObservingFrequency() / 1e9 << "GHz" << debugpost;
298 18 : Double primaryBeamSize = worker_.getPrimaryBeamSize();
299 18 : debuglog<< "primary beam size = " << rad2arcsec(primaryBeamSize) << " arcsec" << debugpost;
300 :
301 : // get table entry sorted by TIME
302 18 : Vector<Double> time(ctiter.time());
303 18 : Cube<Complex> p(ctiter.cparam());
304 18 : Cube<Bool> fl(ctiter.flag());
305 :
306 : // take real part of CPARAM
307 18 : Cube<Float> preal = real(p);
308 :
309 : // execute double-circle gain calibration
310 18 : worker_.doCalibrate(time, preal, fl);
311 :
312 : // if gain calibration fail (typically due to insufficient
313 : // number of data points), go to next iteration step
314 18 : if (preal.empty()) {
315 : // add row numbers to the "TO-BE-REMOVED" list
316 4 : auto rows = ctiter.table().rowNumbers();
317 2 : size_t nelem = to_be_removed.nelements();
318 2 : size_t nelem_add = rows.nelements();
319 2 : to_be_removed.resize(nelem + nelem_add, True);
320 2 : to_be_removed(Slice(nelem, nelem_add)) = rows;
321 2 : ctiter.next();
322 2 : continue;
323 : }
324 :
325 : // set real part of CPARAM
326 16 : setReal(p, preal);
327 :
328 : // record result
329 16 : ctiter.setcparam(p);
330 16 : ctiter.setflag(fl);
331 16 : ctiter.next();
332 : }
333 :
334 : // remove rows registered to the "TO-BE-REMOVED" list
335 9 : if (to_be_removed.nelements() > 0) {
336 1 : ct_->removeRow(to_be_removed);
337 : }
338 9 : }
339 :
340 1658 : void SDDoubleCircleGainCal::keepNCT() {
341 : // Call parent to do general stuff
342 1658 : GJones::keepNCT();
343 :
344 1658 : if (prtlev()>4)
345 0 : cout << " SDDoubleCircleGainCal::keepNCT" << endl;
346 :
347 : // Set proper antenna id
348 3316 : Vector<Int> a1(nAnt());
349 1658 : a1 = currAnt_;
350 : //indgen(a1);
351 :
352 : // We are adding to the most-recently added rows
353 3316 : RefRows rows(ct_->nrow() - nElem(), ct_->nrow() - 1, 1);
354 :
355 : // Write to table
356 3316 : CTMainColumns ncmc(*ct_);
357 1658 : ncmc.antenna1().putColumnCells(rows, a1);
358 1658 : }
359 :
360 0 : void SDDoubleCircleGainCal::syncWtScale()
361 : {
362 0 : currWtScale().resize(currJElem().shape());
363 :
364 : // We use simple (pre-inversion) square of currJElem
365 0 : Cube<Float> cWS(currWtScale());
366 0 : cWS=real(currJElem()*conj(currJElem()));
367 0 : cWS(!currJElemOK())=1.0;
368 0 : }
369 :
370 0 : void SDDoubleCircleGainCal::selfGatherAndSolve(VisSet& vs,
371 : VisEquation& /* ve */) {
372 0 : SDDoubleCircleGainCalImpl sdcalib;
373 0 : debuglog << "SDDoubleCircleGainCal::selfGatherAndSolve()" << debugpost;
374 :
375 : // TODO: implement pre-application of single dish caltables
376 :
377 0 : debuglog << "nspw = " << nSpw() << debugpost;
378 0 : fillNChanParList(msName(), nChanParList());
379 0 : debuglog << "nChanParList=" << nChanParList() << debugpost;
380 :
381 : // Create a new caltable to fill up
382 0 : createMemCalTable();
383 :
384 : // Setup shape of solveAllRPar
385 0 : nElem() = 1;
386 0 : currAnt_.resize(nElem());
387 0 : currAnt_ = -1;
388 0 : initSolvePar();
389 :
390 : // re-initialize calibration solution to 0.0 and calibration flags to false
391 0 : for (Int ispw=0;ispw<nSpw();++ispw) {
392 0 : currSpw() = ispw;
393 0 : solveAllParOK() = false;
394 0 : solveAllCPar() = Complex(0.0);
395 : }
396 :
397 : // Pick up OFF spectra using STATE_ID
398 0 : auto const msSel = vs.iter().ms();
399 : debuglog << "configure data selection for specific calibration mode"
400 0 : << debugpost;
401 0 : auto const taql = selectOnSourceAutoCorr(msSel);
402 0 : debuglog << "taql = \"" << taql << "\"" << debugpost;
403 0 : MeasurementSet msOnSource(tableCommand(taql, msSel).table());
404 : logSink() << LogIO::DEBUGGING << "msSel.nrow()=" << msSel.nrow()
405 0 : << " msOnSource.nrow()=" << msOnSource.nrow() << LogIO::POST;
406 0 : if (msOnSource.nrow() == 0) {
407 0 : throw AipsError("No reference integration in the data.");
408 : }
409 : String dataColName =
410 0 : (msOnSource.tableDesc().isColumn("FLOAT_DATA")) ? "FLOAT_DATA" : "DATA";
411 : logSink() << LogIO::DEBUGGING << "dataColName = " << dataColName
412 0 : << LogIO::POST;
413 :
414 0 : if (msOnSource.tableDesc().isColumn("FLOAT_DATA")) {
415 0 : executeDoubleCircleGainCal<FloatDataColumnAccessor>(msOnSource);
416 : } else {
417 0 : executeDoubleCircleGainCal<DataColumnAccessor>(msOnSource);
418 : }
419 :
420 : //assignCTScanField(*ct_, msName());
421 :
422 : // update weight
423 0 : updateWeight(*ct_);
424 :
425 : // store caltable
426 0 : storeNCT();
427 0 : }
428 :
429 1658 : void SDDoubleCircleGainCal::selfSolveOne(SDBList &sdbs) {
430 : // do nothing at this moment
431 1658 : auto const nSDB = sdbs.nSDB();
432 1658 : debuglog << "nSDB = " << nSDB << debugpost;
433 3316 : for (Int i = 0; i < nSDB; ++i) {
434 1658 : auto const &sdb = sdbs(i);
435 1658 : debuglog << "SDB" << i << ": ";
436 : debuglog << "fld " << sdb.fieldId()
437 : << " ant " << sdb.antenna1()
438 1658 : << " spw " << sdb.spectralWindow();
439 1658 : auto current_precision = cerr.precision();
440 1658 : cerr.precision(16);
441 1658 : debuglog << " time " << sdb.time() << debugpost;
442 1658 : cerr.precision(current_precision);
443 : }
444 1658 : AlwaysAssert(nSDB == 1, AipsError);
445 : // DebugAssert(nSDB == 1, AipsError);
446 :
447 1658 : auto &sdb = sdbs(0);
448 1658 : debuglog << "spw = " << sdb.spectralWindow()[0] << " antenna1 = "
449 3316 : << sdb.antenna1()[0] << "," << sdb.antenna2()[0] << " nRows = "
450 6632 : << sdb.nRows() << debugpost;
451 :
452 4974 : debuglog << "solveAllCPar().shape() = " << solveAllCPar().shape()
453 1658 : << debugpost;
454 1658 : debuglog << "visCube.shape() = " << sdb.visCubeCorrected().shape()
455 3316 : << debugpost;
456 :
457 1658 : auto const &corrected = sdb.visCubeCorrected();
458 1658 : auto const &flag = sdb.flagCube();
459 :
460 1658 : if (!corrected.conform(solveAllCPar())) {
461 : // resize solution array
462 204 : nElem() = sdb.nRows();
463 204 : currAnt_.resize(nElem());
464 204 : sizeSolveParCurrSpw(sdb.nChannels());
465 : }
466 :
467 1658 : solveAllCPar() = Complex(1.0);
468 1658 : solveAllParOK() = false;
469 1658 : currAnt_ = sdb.antenna1();
470 :
471 1658 : size_t const numCorr = corrected.shape()[0];
472 4770 : for (size_t iCorr = 0; iCorr < numCorr; ++iCorr) {
473 3112 : solveAllCPar().yzPlane(iCorr) = corrected.yzPlane(iCorr);
474 3112 : solveParOK().yzPlane(iCorr) = !flag.yzPlane(iCorr);
475 3112 : solveAllParErr().yzPlane(iCorr) = 0.1; // TODO
476 3112 : solveAllParSNR().yzPlane(iCorr) = 1.0; // TODO
477 : }
478 1658 : }
479 :
480 : template<class Accessor>
481 0 : void SDDoubleCircleGainCal::executeDoubleCircleGainCal(
482 : MeasurementSet const &ms) {
483 0 : logSink() << LogOrigin("SDDoubleCircleGainCal", __FUNCTION__, WHERE);
484 : // setup worker class
485 0 : SDDoubleCircleGainCalImpl worker;
486 :
487 0 : Int smoothingSize = -1;// use default smoothing size
488 0 : worker.setCentralRegion(central_disk_size_);
489 0 : if (smooth_) {
490 0 : worker.setSmoothing(smoothingSize);
491 : } else {
492 0 : worker.unsetSmoothing();
493 : }
494 :
495 : // ArrayColumn<Double> uvwColumn(ms, "UVW");
496 : // Matrix<Double> uvw = uvwColumn.getColumn();
497 : // debuglog<< "uvw.shape " << uvw.shape() << debugpost;
498 : //
499 : // make a map between SOURCE_ID and source NAME
500 0 : auto const &sourceTable = ms.source();
501 0 : ScalarColumn<Int> idCol(sourceTable,
502 : sourceTable.columnName(MSSource::MSSourceEnums::SOURCE_ID));
503 0 : ScalarColumn<String> nameCol(sourceTable,
504 : sourceTable.columnName(MSSource::MSSourceEnums::NAME));
505 0 : std::map<Int, String> sourceMap;
506 0 : for (uInt irow = 0; irow < sourceTable.nrow(); ++irow) {
507 0 : auto sourceId = idCol(irow);
508 0 : if (sourceMap.find(sourceId) == sourceMap.end()) {
509 0 : sourceMap[sourceId] = nameCol(irow);
510 : }
511 : }
512 :
513 : // make a map between FIELD_ID and SOURCE_ID
514 0 : auto const &fieldTable = ms.field();
515 0 : idCol.attach(fieldTable,
516 : fieldTable.columnName(MSField::MSFieldEnums::SOURCE_ID));
517 0 : ROArrayMeasColumn<MDirection> dirCol(fieldTable, "REFERENCE_DIR");
518 0 : std::map<Int, Int> fieldMap;
519 0 : for (uInt irow = 0; irow < fieldTable.nrow(); ++irow) {
520 0 : auto sourceId = idCol(irow);
521 0 : fieldMap[static_cast<Int>(irow)] = sourceId;
522 : }
523 :
524 : // access to subtable columns
525 0 : ROScalarQuantColumn<Double> antennaDiameterColumn(ms.antenna(),
526 : "DISH_DIAMETER");
527 0 : ROArrayQuantColumn<Double> observingFrequencyColumn(ms.spectralWindow(),
528 : "CHAN_FREQ");
529 :
530 : // traverse MS
531 0 : Int cols[] = {MS::FIELD_ID, MS::ANTENNA1, MS::FEED1, MS::DATA_DESC_ID};
532 0 : Int *colsp = cols;
533 0 : Block<Int> sortCols(4, colsp, false);
534 0 : MSIter msIter(ms, sortCols, 0.0, false, false);
535 0 : for (msIter.origin(); msIter.more(); msIter++) {
536 0 : MeasurementSet const currentMS = msIter.table();
537 :
538 0 : uInt nrow = currentMS.nrow();
539 0 : debuglog<< "nrow = " << nrow << debugpost;
540 0 : if (nrow == 0) {
541 0 : debuglog<< "Skip" << debugpost;
542 0 : continue;
543 : }
544 0 : Int ispw = msIter.spectralWindowId();
545 0 : if (nChanParList()[ispw] == 4) {
546 : // Skip WVR
547 : debuglog<< "Skip " << ispw
548 0 : << "(nchan " << nChanParList()[ispw] << ")"
549 0 : << debugpost;
550 0 : continue;
551 : }
552 : logSink() << "Process spw " << ispw
553 0 : << "(nchan " << nChanParList()[ispw] << ")" << LogIO::POST;
554 :
555 0 : Int ifield = msIter.fieldId();
556 0 : ScalarColumn<Int> antennaCol(currentMS, "ANTENNA1");
557 : //currAnt_ = antennaCol(0);
558 0 : Int iantenna = antennaCol(0);
559 0 : ScalarColumn<Int> feedCol(currentMS, "FEED1");
560 : debuglog<< "FIELD_ID " << msIter.fieldId()
561 : << " ANTENNA1 " << iantenna//currAnt_
562 : << " FEED1 " << feedCol(0)
563 : << " DATA_DESC_ID " << msIter.dataDescriptionId()
564 0 : << debugpost;
565 :
566 : // setup PointingDirectionCalculator
567 0 : PointingDirectionCalculator pcalc(currentMS);
568 0 : pcalc.setDirectionColumn("DIRECTION");
569 0 : pcalc.setFrame("J2000");
570 0 : pcalc.setDirectionListMatrixShape(PointingDirectionCalculator::ROW_MAJOR);
571 0 : debuglog<< "SOURCE_ID " << fieldMap[ifield] << " SOURCE_NAME " << sourceMap[fieldMap[ifield]] << debugpost;
572 0 : auto const isEphem = ::isEphemeris(sourceMap[fieldMap[ifield]]);
573 0 : Matrix<Double> offset_direction;
574 0 : if (isEphem) {
575 0 : pcalc.setMovingSource(sourceMap[fieldMap[ifield]]);
576 0 : offset_direction = pcalc.getDirection();
577 : } else {
578 0 : pcalc.unsetMovingSource();
579 0 : Matrix<Double> direction = pcalc.getDirection();
580 :
581 : // absolute coordinate -> offset from center
582 0 : OrthographicProjector projector(1.0f);
583 0 : projector.setDirection(direction);
584 0 : Vector<MDirection> md = dirCol.convert(ifield, MDirection::J2000);
585 : //logSink() << "md.shape() = " << md.shape() << LogIO::POST;
586 0 : auto const qd = md[0].getAngle("rad");
587 0 : auto const d = qd.getValue();
588 0 : auto const lat = d[0];
589 0 : auto const lon = d[1];
590 0 : logSink() << "reference coordinate: lat = " << lat << " lon = " << lon << LogIO::POST;
591 0 : projector.setReferencePixel(0.0, 0.0);
592 0 : projector.setReferenceCoordinate(lat, lon);
593 0 : offset_direction = projector.project();
594 0 : auto const pixel_size = projector.pixel_size();
595 : // convert offset_direction from pixel to radian
596 0 : offset_direction *= pixel_size;
597 : }
598 : // debuglog<< "offset_direction = " << offset_direction << debugpost;
599 : // Double const *direction_p = offset_direction.data();
600 : // for (size_t i = 0; i < 10; ++i) {
601 : // debuglog<< "offset_direction[" << i << "]=" << direction_p[i] << debugpost;
602 : // }
603 :
604 0 : ScalarColumn<Double> timeCol(currentMS, "TIME");
605 0 : Vector<Double> time = timeCol.getColumn();
606 0 : Accessor dataCol(currentMS);
607 0 : Cube<Float> data = dataCol.getColumn();
608 0 : ArrayColumn<Bool> flagCol(currentMS, "FLAG");
609 0 : Cube<Bool> flag = flagCol.getColumn();
610 : // debuglog<< "data = " << data << debugpost;
611 :
612 0 : Vector<Double> gainTime;
613 0 : Cube<Float> gain;
614 0 : Cube<Bool> gain_flag;
615 :
616 : // tell some basic information to worker object
617 0 : Quantity antennaDiameterQuant = antennaDiameterColumn(iantenna);
618 0 : worker.setAntennaDiameter(antennaDiameterQuant.getValue("m"));
619 0 : debuglog<< "antenna diameter = " << worker.getAntennaDiameter() << "m" << debugpost;
620 0 : Vector<Quantity> observingFrequencyQuant = observingFrequencyColumn(ispw);
621 0 : Double meanFrequency = 0.0;
622 0 : auto numChan = observingFrequencyQuant.nelements();
623 0 : debuglog<< "numChan = " << numChan << debugpost;
624 0 : assert(numChan > 0);
625 0 : if (numChan % 2 == 0) {
626 0 : meanFrequency = (observingFrequencyQuant[numChan / 2 - 1].getValue("Hz")
627 0 : + observingFrequencyQuant[numChan / 2].getValue("Hz")) / 2.0;
628 : } else {
629 0 : meanFrequency = observingFrequencyQuant[numChan / 2].getValue("Hz");
630 : }
631 : //debuglog << "mean frequency " << meanFrequency.getValue() << " [" << meanFrequency.getFullUnit() << "]" << debugpost;
632 0 : debuglog<< "mean frequency " << meanFrequency << debugpost;
633 0 : worker.setObservingFrequency(meanFrequency);
634 0 : debuglog<< "observing frequency = " << worker.getObservingFrequency() / 1e9 << "GHz" << debugpost;
635 0 : Double primaryBeamSize = worker.getPrimaryBeamSize();
636 0 : debuglog<< "primary beam size = " << rad2arcsec(primaryBeamSize) << " arcsec" << debugpost;
637 :
638 0 : auto const effective_radius = worker.getRadius();
639 : logSink() << "effective radius of the central region = " << effective_radius << " arcsec"
640 0 : << " (" << rad2arcsec(effective_radius) << " arcsec)"<< LogIO::POST;
641 0 : if (worker.isSmoothingActive()) {
642 0 : auto const effective_smoothing_size = worker.getEffectiveSmoothingSize();
643 0 : logSink() << "smoothing activated. effective size = " << effective_smoothing_size << LogIO::POST;
644 : }
645 : else {
646 0 : logSink() << "smoothing deactivated." << LogIO::POST;
647 : }
648 :
649 : // execute calibration
650 0 : worker.calibrate(data, flag, time, offset_direction, gainTime, gain, gain_flag);
651 : //debuglog<< "gain_time = " << gain_time << debugpost;
652 : //debuglog<<"gain = " << gain << debugpost;
653 0 : size_t numGain = gainTime.nelements();
654 0 : debuglog<< "number of gain " << numGain << debugpost;
655 :
656 0 : currSpw() = ispw;
657 :
658 : // make sure storage and flag for calibration solution allocated
659 : // for the current spw are properly initialized
660 0 : solveAllCPar() = Complex(0.0);
661 0 : solveAllParOK() = false;
662 :
663 0 : currField() = ifield;
664 0 : currAnt_ = iantenna;
665 0 : debuglog << "antenna is " << currAnt_ << debugpost;
666 :
667 0 : size_t numCorr = gain.shape()[0];
668 0 : Slice const chanSlice(0, numChan);
669 0 : for (size_t i = 0; i < numGain; ++i) {
670 0 : refTime() = gainTime[i];
671 0 : Slice const rowSlice(i, 1);
672 0 : for (size_t iCorr = 0; iCorr < numCorr; ++iCorr) {
673 0 : Slice const corrSlice(iCorr, 1);
674 0 : auto cparSlice = solveAllCPar()(corrSlice, chanSlice, Slice(0, 1));
675 0 : convertArray(cparSlice, gain(corrSlice, chanSlice, rowSlice));
676 0 : solveAllParOK()(corrSlice, chanSlice, Slice(0, 1)) = !gain_flag(corrSlice, chanSlice, rowSlice);
677 0 : solveAllParErr().yzPlane(iCorr) = 0.1; // TODO
678 0 : solveAllParSNR().yzPlane(iCorr) = 1.0; // TODO
679 : }
680 :
681 0 : keepNCT();
682 : }
683 :
684 : }
685 0 : }
686 : }
|
