Line data Source code
1 : //# CalSet.cc: Implementation of Calibration parameter cache
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 <algorithm>
28 : #include <synthesis/CalTables/CalSet.h>
29 :
30 : #include <synthesis/CalTables/CalTable2.h>
31 : //#include <synthesis/CalTables/CalTable.h>
32 : #include <synthesis/CalTables/CalDescColumns2.h>
33 : #include <synthesis/CalTables/SolvableVJTable.h>
34 : #include <synthesis/CalTables/CalMainColumns2.h>
35 :
36 : #include <synthesis/CalTables/TimeVarVJDesc.h>
37 : #include <synthesis/CalTables/SolvableVJDesc.h>
38 : #include <synthesis/CalTables/SolvableVJMRec.h>
39 : #include <synthesis/CalTables/SolvableVJMCol.h>
40 :
41 : #include <casacore/tables/Tables/TableDesc.h>
42 : #include <casacore/tables/Tables/SetupNewTab.h>
43 : #include <casacore/tables/Tables/Table.h>
44 : #include <casacore/tables/Tables/ScaColDesc.h>
45 : #include <casacore/tables/Tables/ArrColDesc.h>
46 : #include <casacore/tables/Tables/ScalarColumn.h>
47 : #include <casacore/tables/Tables/ArrayColumn.h>
48 : #include <casacore/tables/Tables/RefRows.h>
49 :
50 : #include <casacore/casa/Arrays.h>
51 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
52 : #include <casacore/casa/BasicSL/String.h>
53 : #include <casacore/casa/Utilities/Assert.h>
54 : #include <casacore/casa/OS/Timer.h>
55 : #include <casacore/casa/Exceptions/Error.h>
56 : #include <casacore/casa/OS/Path.h>
57 : #include <casacore/casa/sstream.h>
58 :
59 : #include <casacore/casa/Logging/LogMessage.h>
60 : #include <casacore/casa/Logging/LogSink.h>
61 :
62 : namespace casa { //# NAMESPACE CASA - BEGIN
63 :
64 : // ------------------------------------------------------------------
65 :
66 : // From shape only, this is the solve context
67 0 : template<class T> CalSet<T>::CalSet(const casacore::Int& nSpw) :
68 : calTableName_(""),
69 : nSpw_(nSpw),
70 : nPar_(0),
71 : nChan_(0),
72 : nElem_(0),
73 : nTime_(0),
74 0 : spwOK_(nSpw_,false),
75 : startChan_(nSpw_,0),
76 0 : freq_(nSpw_,NULL),
77 0 : MJDStart_(nSpw_,NULL),
78 0 : MJDStop_(nSpw_,NULL),
79 0 : MJDTimeStamp_(nSpw_,NULL),
80 0 : fieldId_(nSpw_,NULL),
81 0 : fieldName_(nSpw_,NULL),
82 0 : sourceName_(nSpw_,NULL),
83 0 : par_(nSpw_,NULL),
84 0 : parOK_(nSpw_,NULL),
85 0 : parErr_(nSpw_,NULL),
86 0 : parSNR_(nSpw_,NULL),
87 : // iSolutionOK_(nSpw_,NULL),
88 0 : iFit_(nSpw_,NULL),
89 0 : iFitwt_(nSpw_,NULL),
90 0 : solutionOK_(nSpw_,NULL),
91 0 : fit_(nSpw_,NULL),
92 0 : fitwt_(nSpw_,NULL)
93 : {
94 0 : calTabDesc_=NULL;
95 0 : calTab_=NULL;
96 0 : svjmcol_=NULL;
97 0 : };
98 :
99 : // From shape only, this is the solve context
100 : template<class T> CalSet<T>::CalSet(const casacore::Int& nSpw,
101 : const casacore::Int& nPar,
102 : const casacore::Vector<casacore::Int>& nChan,
103 : const casacore::Int& nElem,
104 : const casacore::Vector<casacore::Int>& nTime) :
105 : calTableName_(""),
106 : nSpw_(nSpw),
107 : nPar_(nPar),
108 : nChan_(nChan),
109 : nElem_(nElem),
110 : nTime_(nTime),
111 : spwOK_(nSpw_,false),
112 : startChan_(nSpw_,0),
113 : freq_(nSpw_,NULL),
114 : MJDStart_(nSpw_,NULL),
115 : MJDStop_(nSpw_,NULL),
116 : MJDTimeStamp_(nSpw_,NULL),
117 : fieldId_(nSpw_,NULL),
118 : fieldName_(nSpw_,NULL),
119 : sourceName_(nSpw_,NULL),
120 : par_(nSpw_,NULL),
121 : parOK_(nSpw_,NULL),
122 : parErr_(nSpw_,NULL),
123 : parSNR_(nSpw_,NULL),
124 : // iSolutionOK_(nSpw_,NULL),
125 : iFit_(nSpw_,NULL),
126 : iFitwt_(nSpw_,NULL),
127 : solutionOK_(nSpw_,NULL),
128 : fit_(nSpw_,NULL),
129 : fitwt_(nSpw_,NULL)
130 : {
131 : calTabDesc_=NULL;
132 : calTab_=NULL;
133 : svjmcol_=NULL;
134 : // Resize caches
135 : inflate();
136 : };
137 :
138 :
139 : // From existing CalTable name, this is apply context
140 0 : template<class T> CalSet<T>::CalSet(const casacore::String& calTableName,
141 : const casacore::String& select,
142 : const casacore::Int& nSpw,
143 : const casacore::Int& nPar,
144 : const casacore::Int& nElem) :
145 : calTableName_(calTableName),
146 : nSpw_(nSpw),
147 : nPar_(nPar),
148 : nChan_(nSpw,0),
149 : nElem_(nElem),
150 : nTime_(nSpw,0),
151 0 : spwOK_(nSpw_,false),
152 : startChan_(nSpw_,0),
153 0 : freq_(nSpw_,NULL),
154 0 : MJDStart_(nSpw_,NULL),
155 0 : MJDStop_(nSpw_,NULL),
156 0 : MJDTimeStamp_(nSpw_,NULL),
157 0 : fieldId_(nSpw_,NULL),
158 0 : fieldName_(nSpw_,NULL),
159 0 : sourceName_(nSpw_,NULL),
160 0 : par_(nSpw_,NULL),
161 0 : parOK_(nSpw_,NULL),
162 0 : parErr_(nSpw_,NULL),
163 0 : parSNR_(nSpw_,NULL),
164 : // iSolutionOK_(nSpw_,NULL),
165 0 : iFit_(nSpw_,NULL),
166 0 : iFitwt_(nSpw_,NULL),
167 0 : solutionOK_(nSpw_,NULL),
168 0 : fit_(nSpw_,NULL),
169 0 : fitwt_(nSpw_,NULL)
170 : {
171 0 : calTabDesc_=NULL;
172 0 : calTab_=NULL;
173 0 : svjmcol_=NULL;
174 :
175 : // Fill from table
176 0 : load(calTableName,select);
177 :
178 : // Set spwOK_ according to available calibration
179 0 : setSpwOK();
180 :
181 0 : }
182 :
183 0 : template<class T> CalSet<T>::~CalSet() {
184 0 : deflate();
185 0 : }
186 :
187 0 : template<class T> void CalSet<T>::resize(const casacore::Int& nPar,
188 : const casacore::Vector<casacore::Int>& nChan,
189 : const casacore::Int& nElem,
190 : const casacore::Vector<casacore::Int>& nTime) {
191 0 : nPar_=nPar;
192 0 : nChan_=nChan;
193 0 : nElem_=nElem;
194 0 : nTime_=nTime;
195 :
196 0 : inflate();
197 :
198 0 : }
199 :
200 :
201 :
202 : // Inflate cache to proper size
203 0 : template<class T> void CalSet<T>::inflate() {
204 :
205 : // Construct shaped pointed-to objects in cache
206 :
207 : // TODO:
208 : // Consider initialization value (per type)
209 :
210 : // Delete exiting cache
211 0 : deflate();
212 :
213 0 : for (casacore::Int ispw=0; ispw<nSpw_; ispw++) {
214 0 : casacore::uInt ntime=nTime_(ispw);
215 0 : if (ntime > 0) {
216 :
217 0 : freq_[ispw] = new casacore::Vector<casacore::Double>(nChan_(ispw),0.0);
218 :
219 0 : MJDStart_[ispw] = new casacore::Vector<casacore::Double>(ntime,0.0);
220 0 : MJDStop_[ispw] = new casacore::Vector<casacore::Double>(ntime,0.0);
221 0 : MJDTimeStamp_[ispw] = new casacore::Vector<casacore::Double>(ntime,0.0);
222 0 : fieldName_[ispw] = new casacore::Vector<casacore::String>(ntime,"");
223 0 : sourceName_[ispw] = new casacore::Vector<casacore::String>(ntime,"");
224 0 : fieldId_[ispw] = new casacore::Vector<casacore::Int>(ntime,-1);
225 :
226 0 : casacore::IPosition parshape(4,nPar_,nChan_(ispw),nElem_,ntime);
227 0 : par_[ispw] = new casacore::Array<T>(parshape,(T)1.0);
228 0 : parOK_[ispw] = new casacore::Array<casacore::Bool>(parshape,false);
229 0 : parErr_[ispw] = new casacore::Array<casacore::Float>(parshape,0.0f);
230 0 : parSNR_[ispw] = new casacore::Array<casacore::Float>(parshape,0.0f);
231 :
232 : // iSolutionOK_[ispw] = new casacore::Matrix<casacore::Bool>(nElem_,ntime,false);
233 0 : iFit_[ispw] = new casacore::Matrix<casacore::Float>(nElem_,ntime,0.0f);
234 0 : iFitwt_[ispw] = new casacore::Matrix<casacore::Float>(nElem_,ntime,0.0f);
235 0 : solutionOK_[ispw] = new casacore::Vector<casacore::Bool>(ntime,false);
236 0 : fit_[ispw] = new casacore::Vector<casacore::Float>(ntime,0.0f);
237 0 : fitwt_[ispw] = new casacore::Vector<casacore::Float>(ntime,0.0f);
238 : }
239 : }
240 :
241 0 : }
242 :
243 :
244 :
245 :
246 0 : template<class T> void CalSet<T>::deflate() {
247 :
248 : // Delete parameter memory
249 :
250 0 : for (casacore::Int ispw=0; ispw<nSpw_; ispw++) {
251 0 : if (MJDStart_[ispw]) delete MJDStart_[ispw];
252 0 : if (MJDStop_[ispw]) delete MJDStop_[ispw];
253 0 : if (MJDTimeStamp_[ispw]) delete MJDTimeStamp_[ispw];
254 0 : if (fieldName_[ispw]) delete fieldName_[ispw];
255 0 : if (sourceName_[ispw]) delete sourceName_[ispw];
256 0 : if (fieldId_[ispw]) delete fieldId_[ispw];
257 0 : if (par_[ispw]) delete par_[ispw];
258 0 : if (parOK_[ispw]) delete parOK_[ispw];
259 0 : if (parErr_[ispw]) delete parErr_[ispw];
260 0 : if (parSNR_[ispw]) delete parSNR_[ispw];
261 : // if (iSolutionOK_[ispw]) delete iSolutionOK_[ispw];
262 0 : if (iFit_[ispw]) delete iFit_[ispw];
263 0 : if (iFitwt_[ispw]) delete iFitwt_[ispw];
264 0 : if (solutionOK_[ispw]) delete solutionOK_[ispw];
265 0 : if (fit_[ispw]) delete fit_[ispw];
266 0 : if (fitwt_[ispw]) delete fitwt_[ispw];
267 0 : MJDStart_[ispw]=NULL;
268 0 : MJDStop_[ispw]=NULL;
269 0 : MJDTimeStamp_[ispw]=NULL;
270 0 : fieldName_[ispw]=NULL;
271 0 : sourceName_[ispw]=NULL;
272 0 : fieldId_[ispw]=NULL;
273 0 : par_[ispw]=NULL;
274 0 : parOK_[ispw]=NULL;
275 0 : parErr_[ispw]=NULL;
276 0 : parSNR_[ispw]=NULL;
277 : // iSolutionOK_[ispw]=NULL;
278 0 : iFit_[ispw]=NULL;
279 0 : iFitwt_[ispw]=NULL;
280 0 : solutionOK_[ispw]=NULL;
281 0 : fit_[ispw]=NULL;
282 0 : fitwt_[ispw]=NULL;
283 : }
284 0 : }
285 :
286 :
287 :
288 0 : template<class T> void CalSet<T>::load (const casacore::String& file,
289 : const casacore::String& select)
290 : {
291 : // Load data from a calibration table
292 : // casacore::Input:
293 : // file const casacore::String& Cal table name
294 : // select const casacore::String& Selection string
295 : //
296 :
297 :
298 : // cout << "CalSet::load...(nSpw_=" << nSpw_ << ")" << endl;
299 :
300 0 : casacore::Timer timer;
301 0 : timer.mark();
302 :
303 0 : casacore::LogMessage message(casacore::LogOrigin("CalSet","load"));
304 :
305 : // Decode the Jones matrix type
306 : /*
307 : casacore::Int jonesType = 0;
308 : if (nPar_ == 1) jonesType = 1;
309 : if (nPar_ == 2) jonesType = 2;
310 : if (nPar_ == 4) jonesType = 3;
311 : */
312 :
313 : // Open, select, sort the calibration table
314 0 : SolvableVisJonesTable svjtab(file);
315 0 : svjtab.select2(select);
316 :
317 : // Get no. of antennas and time slots
318 0 : casacore::Int numberAnt = svjtab.maxAntenna() + 1;
319 :
320 : // cout << "Initial selection: " << timer.all_usec()/1.0e6 << endl;
321 0 : timer.mark();
322 :
323 0 : AlwaysAssert(numberAnt==nElem_,casacore::AipsError)
324 :
325 0 : casacore::Int nDesc=svjtab.nRowDesc();
326 0 : casacore::Vector<casacore::Int> spwmap(nDesc,-1);
327 :
328 0 : casacore::Vector<casacore::Int> nRowPerCDI;
329 0 : svjtab.rowsPerCalDescId(nRowPerCDI);
330 :
331 : // cout << "Optimized CDI count: " << timer.all_usec()/1.0e6 << endl;
332 0 : timer.mark();
333 :
334 0 : for (casacore::Int idesc=0;idesc<nDesc;idesc++) {
335 :
336 : // This cal desc
337 0 : CalDescRecord* calDescRec = new CalDescRecord (svjtab.getRowDesc(idesc));
338 :
339 : // Get this spw ID
340 0 : casacore::Vector<casacore::Int> spwlist;
341 0 : calDescRec->getSpwId(spwlist);
342 0 : casacore::Int nSpw; spwlist.shape(nSpw);
343 0 : if (nSpw > 1) {}; // ERROR!!! Should only be one spw per cal desc!
344 0 : spwmap(idesc)=spwlist(0);
345 :
346 : // Trap spwids that do not occur in the MS
347 : // (Since we rely on the casacore::MS meta info for the calibration solutions,
348 : // we cannot identify spws that do not occur in the MS.)
349 0 : if (spwlist(0)>nSpw_-1)
350 : throw(casacore::AipsError("Caltable '"+file+"' contains spw = "+
351 0 : casacore::String::toString(spwlist(0))+
352 0 : " which does not occur in the MS. Cannot proceed."));
353 :
354 : // In next few steps, need to watch for repeat spws in new cal descs!!
355 :
356 : // Get number of channels this spw
357 0 : casacore::Vector<casacore::Int> nchanlist;
358 :
359 0 : calDescRec->getNumChan(nchanlist);
360 0 : nChan_(spwmap(idesc))=nchanlist(0);
361 :
362 : // Get channel range / start channel
363 0 : casacore::Cube<casacore::Int> chanRange;
364 0 : calDescRec->getChanRange(chanRange);
365 0 : startChan_(spwmap(idesc))=chanRange(0,0,0);
366 :
367 : // Get slot count for this desc's spw
368 0 : nTime_(spwmap(idesc))=nRowPerCDI(idesc)/numberAnt;
369 :
370 0 : delete calDescRec;
371 : }
372 :
373 : // cout << "CalDesc sizeup: " << timer.all_usec()/1.0e6 << endl;
374 : // cout << "nTime_ = " << nTime_ << endl;
375 :
376 0 : timer.mark();
377 :
378 : // At this point, we know how big our slot-dep caches must be
379 : // (in private data), so initialize them
380 0 : inflate();
381 :
382 : // Remember if we found and filled any solutions
383 0 : casacore::Bool solfillok(false);
384 :
385 : // cout << "CalSet inflated: " << timer.all_usec()/1.0e6 << endl;
386 :
387 :
388 : // Fill per caldesc
389 0 : casacore::Double ttime(0.0);
390 0 : for (casacore::Int idesc=0;idesc<nDesc;idesc++) {
391 :
392 : // cout << "CDI = " << idesc << " "<< flush;
393 :
394 0 : timer.mark();
395 :
396 0 : casacore::Int thisSpw=spwmap(idesc);
397 :
398 : // Reopen and globally select caltable
399 : //SolvableVisJonesTable svjtabspw(file);
400 0 : CalTable2 svjtabspw(file);
401 0 : svjtabspw.select2(select);
402 :
403 : // cout << " Sel: " << timer.all_usec()/1.0e6 << flush;
404 :
405 : // isolate this caldesc:
406 0 : std::ostringstream selectstr;
407 0 : selectstr << "CAL_DESC_ID == " << idesc;
408 0 : casacore::String caldescsel; caldescsel = selectstr.str();
409 0 : svjtabspw.select2(caldescsel);
410 :
411 : // cout << " CDIsel: " << timer.all_usec()/1.0e6 << flush;
412 :
413 : // cout << "Sorting..." << endl;
414 0 : casacore::Block<casacore::String> scol(1);
415 0 : scol[0]="TIME";
416 0 : svjtabspw.sort2(scol);
417 :
418 : // cout << " casacore::Sort: " << timer.all_usec()/1.0e6 << flush;
419 :
420 0 : casacore::Int nrow = svjtabspw.nRowMain();
421 0 : casacore::IPosition out(3,0,0,0); // par, chan, row
422 0 : casacore::IPosition in(4,0,0,0,0); // par, chan, ant, slot
423 0 : if (nrow>0) {
424 :
425 : // Found some solutions to fill
426 0 : solfillok=true;
427 :
428 : // Ensure sorted on time
429 0 : casacore::Block <casacore::String> sortCol(1,"TIME");
430 0 : svjtabspw.sort2(sortCol);
431 :
432 : // Extract the gain table columns
433 : //ROSolvableVisJonesMCol svjmcol(svjtabspw);
434 : ROSolvableCalSetMCol<T> *svjmcol;
435 0 : svjmcol= new ROSolvableCalSetMCol<T>(svjtabspw);
436 :
437 0 : casacore::Vector<casacore::Int> calDescId; (*svjmcol).calDescId().getColumn(calDescId);
438 0 : casacore::Vector<casacore::Double> time; (*svjmcol).time().getColumn(time);
439 0 : casacore::Vector<casacore::Double> interval; (*svjmcol).interval().getColumn(interval);
440 0 : casacore::Vector<casacore::Int> antenna1; (*svjmcol).antenna1().getColumn(antenna1);
441 : // cout << "antennas = " << antenna1 << endl;
442 0 : casacore::Vector<casacore::Int> fieldId; (*svjmcol).fieldId().getColumn(fieldId);
443 : // cout << "fieldId = " << fieldId << endl;
444 0 : casacore::Vector<casacore::String> fieldName; (*svjmcol).fieldName().getColumn(fieldName);
445 0 : casacore::Vector<casacore::String> sourceName; (*svjmcol).sourceName().getColumn(sourceName);
446 0 : casacore::Vector<casacore::Bool> totalSolOk; (*svjmcol).totalSolnOk().getColumn(totalSolOk);
447 0 : casacore::Vector<casacore::Float> totalFit; (*svjmcol).totalFit().getColumn(totalFit);
448 0 : casacore::Vector<casacore::Float> totalFitWt; (*svjmcol).totalFitWgt().getColumn(totalFitWt);
449 0 : casacore::Array<T> gain; (*svjmcol).gain().getColumn(gain);
450 0 : casacore::Cube<casacore::Bool> solOk; (*svjmcol).solnOk().getColumn(solOk);
451 0 : casacore::Cube<casacore::Float> fit; (*svjmcol).fit().getColumn(fit);
452 0 : casacore::Cube<casacore::Float> fitWt; (*svjmcol).fitWgt().getColumn(fitWt);
453 0 : casacore::Cube<casacore::Bool> flag; (*svjmcol).flag().getColumn(flag);
454 0 : casacore::Cube<casacore::Float> snr; (*svjmcol).snr().getColumn(snr);
455 :
456 : // Read the calibration information
457 0 : casacore::Double /*lastTime(-1.0),*/ thisTime(0.0), thisInterval(0.0);
458 0 : casacore::Int islot(-1);
459 : casacore::Int iant;
460 :
461 0 : for (casacore::Int irow=0; irow<nrow; irow++) {
462 0 : out(2)=irow;
463 :
464 0 : thisTime=time(irow);
465 :
466 : // If this is a new solution
467 0 : if (irow%numberAnt==0) {
468 :
469 0 : islot++;
470 0 : in(3)=islot;
471 :
472 0 : thisInterval=interval(irow);
473 0 : (*MJDTimeStamp_[thisSpw])(islot) = thisTime;
474 0 : (*MJDStart_[thisSpw])(islot) = thisTime - thisInterval / 2.0;
475 0 : (*MJDStop_[thisSpw])(islot) = thisTime + thisInterval / 2.0;
476 0 : (*fieldId_[thisSpw])(islot) = fieldId(irow);
477 0 : (*fieldName_[thisSpw])(islot) = fieldName(irow);
478 0 : (*sourceName_[thisSpw])(islot) = sourceName(irow);
479 :
480 0 : (*solutionOK_[thisSpw])(islot) = totalSolOk(irow);
481 0 : (*fit_[thisSpw])(islot) = totalFit(irow);
482 0 : (*fitwt_[thisSpw])(islot) = totalFitWt(irow);
483 :
484 : //lastTime = thisTime;
485 : };
486 :
487 0 : iant=antenna1(irow);
488 0 : in(2)=iant;
489 :
490 : // (*iSolutionOK_[thisSpw])(iant,islot) = solOk(0,0,irow);
491 0 : (*iFit_[thisSpw])(iant,islot) = fit(0,0,irow);
492 0 : (*iFitwt_[thisSpw])(iant,islot) = fitWt(0,0,irow);
493 :
494 0 : for (casacore::Int ichan=0; ichan<nChan_(thisSpw); ichan++) {
495 :
496 :
497 0 : out(1)=in(1)=ichan;
498 :
499 0 : for (casacore::Int ipar=0; ipar<nPar_; ipar++) {
500 0 : in(0)=out(0)=ipar;
501 0 : (*par_[thisSpw])(in)=gain(out);
502 0 : (*parOK_[thisSpw])(in) = (solOk(out) && !flag(out));
503 0 : (*parSNR_[thisSpw])(in) = snr(out);
504 : }
505 : }
506 :
507 : } // irow
508 : } // nrow>0
509 :
510 : // cout << "cs fields = " << *fieldId_[thisSpw] << endl;
511 :
512 0 : casacore::Double itime=timer.all_usec()/1.0e6;
513 0 : ttime+=itime;
514 :
515 : // cout << " Totals: " << itime << " " << ttime << endl;
516 :
517 : } // idesc
518 :
519 : // If we found no solutions in selected table, abort:
520 0 : if (!solfillok) {
521 0 : throw(casacore::AipsError(" Specified cal table selection selects no solutions in this table. Please review setapply settings."));
522 : }
523 :
524 :
525 0 : };
526 :
527 0 : template<class T> void CalSet<T>::initCalTableDesc(const casacore::String& type, const casacore::Int& parType)
528 : {
529 0 : if (calTabDesc_) {delete calTabDesc_;calTabDesc_=NULL;}
530 0 : calTabDesc_=new CalTableDesc2(type,parType);
531 0 : }
532 :
533 0 : template<class T> void CalSet<T>::attach()
534 : {
535 0 : }
536 :
537 :
538 0 : template<class T> void CalSet<T>::store (const casacore::String& file,
539 : const casacore::String& type,
540 : const casacore::Bool& append,
541 : const casacore::String& msname)
542 : {
543 : // Write the solutions to an output calibration table
544 : // casacore::Input:
545 : // file casacore::String Cal table name
546 : // append casacore::Bool Append if true, else overwrite
547 : //
548 :
549 : // total rows to be written per Spw
550 0 : casacore::Vector<casacore::Int> nRow(nSpw_,0);
551 0 : for (casacore::Int iSpw=0;iSpw<nSpw_;iSpw++)
552 0 : if (solutionOK_[iSpw])
553 0 : nRow(iSpw)=nElem_*ntrue(*(solutionOK_[iSpw]));
554 :
555 : // Escape if nothing to write
556 0 : if (sum(nRow)==0)
557 0 : throw(casacore::AipsError("No valid calibration solutions; no table written."));
558 :
559 : // Initialization:
560 : // No. of rows in cal_main, cal_desc and cal_history
561 0 : casacore::Int nMain = 0;
562 0 : casacore::Int nDesc = 0;
563 : //casacore::Int nHist = 0;
564 :
565 0 : if (calTabDesc_ == NULL)
566 : {
567 0 : std::ostringstream str;
568 : str << "CalSet::store(): Perhaps CalSet.initCalTableDesc() was not called "
569 : << "before calling CalSet::store()?"
570 0 : << " Jones = " << type << " File = " << file;
571 :
572 0 : throw(casacore::AipsError(str.str()));
573 : }
574 :
575 0 : calTabDesc_->addDesc(calTabDesc_->defaultFitDesc(),calTabDesc_->calMainDesc());
576 :
577 : // Calibration table
578 : // SolvableVisJonesTable *tab;
579 :
580 : // Open the output file if it already exists and is being appended to.
581 0 : if (calTab_) delete calTab_; calTab_=NULL;
582 :
583 0 : if (append && casacore::Table::isWritable (file)) {
584 : // tab = new SolvableVisJonesTable (file, casacore::Table::Update);
585 0 : calTab_ = new CalTable2 (file, *calTabDesc_, casacore::Table::Update);
586 0 : nMain = calTab_->nRowMain();
587 0 : nDesc = calTab_->nRowDesc();
588 : //nHist = calTab_->nRowHistory();
589 : } else {
590 : // Create a new calibration table
591 0 : casacore::Table::TableOption access = casacore::Table::New;
592 : // tab = new SolvableVisJonesTable (file, type, access);
593 0 : calTab_ = new CalTable2 (file, *calTabDesc_, access);
594 : };
595 :
596 : // Write every spw w/ max number of channels
597 : // (eventually, CalTable should permit variable-shape cols)
598 0 : casacore::Int maxNumChan(1);
599 0 : for (casacore::Int iSpw=0; iSpw<nSpw_; iSpw++)
600 0 : if (par_[iSpw]!=NULL)
601 0 : maxNumChan=std::max(maxNumChan,nChan_(iSpw));
602 :
603 : // Some default values
604 0 : casacore::Double dzero = 0;
605 0 : casacore::IPosition ip(2,1,maxNumChan);
606 :
607 : // CalDesc Sub-table records
608 : CalDescRecord* descRec;
609 0 : casacore::Vector<casacore::Int> calDescNum(nSpw_); calDescNum=-1;
610 0 : for (casacore::Int iSpw=0; iSpw<nSpw_; iSpw++) {
611 :
612 : // Write a CalDesc for each spw which has solutions
613 : // Note: CalDesc index != SpwId, in general
614 :
615 0 : if (nRow(iSpw)>0 && par_[iSpw]!=NULL) {
616 :
617 : // Access to existing CAL_DESC columns:
618 0 : CalDescColumns2 cd(*calTab_);
619 :
620 : // Check if this spw already in CAL_DESC
621 : // cout << "spwCol = " << cd.spwId().getColumn() << endl;
622 :
623 0 : casacore::Bool newCD(true);
624 0 : for (casacore::Int iCD=0;iCD<nDesc;iCD++) {
625 :
626 0 : casacore::IPosition iCDip(1,0);
627 0 : if ( iSpw==(cd.spwId()(iCD))(iCDip) ) {
628 : // Don't need new CAL_DESC entry
629 0 : newCD=false;
630 0 : calDescNum(iSpw)=iCD;
631 0 : break;
632 : }
633 : }
634 :
635 0 : if (newCD) {
636 :
637 : // Cal_desc fields
638 0 : casacore::Vector <casacore::Int> spwId(1,iSpw);
639 0 : casacore::Matrix <casacore::Double> chanFreq(ip, dzero);
640 0 : casacore::Matrix <casacore::Double> chanWidth(ip, dzero);
641 0 : casacore::Array <casacore::String> polznType(ip, casacore::String(""));
642 0 : casacore::Cube <casacore::Int> chanRange(casacore::IPosition(3,2,1,maxNumChan), 0);
643 0 : casacore::Vector <casacore::Int> numChan(1,nChan_(iSpw));
644 0 : for (casacore::Int ichan=0; ichan<nChan_(iSpw); ichan++) {
645 0 : chanRange(0,0,ichan)=startChan_(iSpw);
646 0 : chanRange(1,0,ichan)=startChan_(iSpw) + nChan_(iSpw) -1;
647 : }
648 :
649 : // Fill the cal_desc record
650 0 : descRec = new CalDescRecord;
651 0 : descRec->defineNumSpw (1);
652 0 : descRec->defineNumChan (numChan);
653 0 : descRec->defineNumReceptors (2);
654 0 : descRec->defineNJones (2);
655 0 : descRec->defineSpwId (spwId);
656 0 : descRec->defineChanFreq (chanFreq);
657 0 : descRec->defineChanWidth (chanWidth);
658 0 : descRec->defineChanRange (chanRange);
659 0 : descRec->definePolznType (polznType);
660 0 : descRec->defineJonesType ("full");
661 0 : descRec->defineMSName (casacore::Path(msname).baseName());
662 : // descRec->defineMSName ("");
663 :
664 : // Write the cal_desc record
665 :
666 0 : calTab_->putRowDesc (nDesc, *descRec);
667 0 : delete descRec;
668 :
669 : // This spw will have this calDesc index in main table
670 0 : calDescNum(iSpw) = nDesc;
671 0 : nDesc++;
672 : }
673 :
674 : }
675 :
676 : }
677 :
678 :
679 : // Now write MAIN table in column-wise fashion
680 :
681 : // Starting row in this slot
682 :
683 0 : for (casacore::Int iSpw=0; iSpw<nSpw_; iSpw++) {
684 :
685 : // Write table for spws which have solutions
686 0 : if (par_[iSpw]!=NULL) {
687 :
688 : // Create references to cal data for this spw
689 0 : casacore::Vector<casacore::Bool> thisSolOK; thisSolOK.reference(*(solutionOK_[iSpw]));
690 0 : casacore::Vector<casacore::Double> thisMJDTimeStamp; thisMJDTimeStamp.reference(*(MJDTimeStamp_[iSpw]));
691 0 : casacore::Vector<casacore::Double> thisMJDStart; thisMJDStart.reference(*(MJDStart_[iSpw]));
692 0 : casacore::Vector<casacore::Double> thisMJDStop; thisMJDStop.reference(*(MJDStop_[iSpw]));
693 0 : casacore::Vector<casacore::Int> thisFieldId; thisFieldId.reference(*(fieldId_[iSpw]));
694 0 : casacore::Vector<casacore::String> thisFieldName; thisFieldName.reference(*(fieldName_[iSpw]));
695 0 : casacore::Vector<casacore::String> thisSourceName; thisSourceName.reference(*(sourceName_[iSpw]));
696 0 : casacore::Vector<casacore::Float> thisFit; thisFit.reference(*(fit_[iSpw]));
697 0 : casacore::Vector<casacore::Float> thisFitwt; thisFitwt.reference(*(fitwt_[iSpw]));
698 : // casacore::Array<casacore::Complex> thisAntGain; thisAntGain.reference(*(par_[iSpw]));
699 0 : casacore::Array<T> thisAntGain; thisAntGain.reference(*(par_[iSpw]));
700 : // casacore::Matrix<casacore::Bool> thisISolutionOK; thisISolutionOK.reference(*(iSolutionOK_[iSpw]));
701 0 : casacore::Matrix<casacore::Float> thisIFit; thisIFit.reference(*(iFit_[iSpw]));
702 0 : casacore::Matrix<casacore::Float> thisIFitwt; thisIFitwt.reference(*(iFitwt_[iSpw]));
703 :
704 0 : casacore::Int thisnRow=nRow(iSpw);
705 :
706 : // Only if there are rows to add to table
707 0 : if (thisnRow > 0) {
708 :
709 : // These are constant columns (with boring values, currently)
710 0 : casacore::Vector<casacore::Double> timeEP(thisnRow,0.0);
711 0 : casacore::Vector<casacore::Int> feed1(thisnRow,0);
712 0 : casacore::Vector<casacore::Int> arrayId(thisnRow,0);
713 0 : casacore::Vector<casacore::Int> obsId(thisnRow,0);
714 0 : casacore::Vector<casacore::Int> scanNum(thisnRow,0);
715 0 : casacore::Vector<casacore::Int> procId(thisnRow,0);
716 0 : casacore::Vector<casacore::Int> stateId(thisnRow,0);
717 0 : casacore::Vector<casacore::Int> phaseId(thisnRow,0);
718 0 : casacore::Vector<casacore::Int> pulsarBin(thisnRow,0);
719 0 : casacore::Vector<casacore::Int> pulsarGateId(thisnRow,0);
720 0 : casacore::Vector<casacore::Int> freqGroup(thisnRow,0);
721 0 : casacore::Vector<casacore::Int> calHistId(thisnRow,0);
722 :
723 : // This is constant
724 0 : casacore::Vector<casacore::Int> calDescId(thisnRow,calDescNum(iSpw));
725 :
726 : // These are constant per slot
727 : // (these cols should be incremental)
728 0 : casacore::Vector<casacore::Double> time(thisnRow,0.0);
729 0 : casacore::Vector<casacore::Double> interval(thisnRow,0.0);
730 0 : casacore::Vector<casacore::Int> fieldId(thisnRow,0);
731 0 : casacore::Vector<casacore::String> fieldName(thisnRow,"");
732 0 : casacore::Vector<casacore::String> sourceName(thisnRow,"");
733 0 : casacore::Vector<casacore::Bool> totalSolOk(thisnRow,false);
734 0 : casacore::Vector<casacore::Float> totalFit(thisnRow,0.0);
735 0 : casacore::Vector<casacore::Float> totalFitWt(thisnRow,0.0);
736 :
737 : // These vary
738 0 : casacore::Vector<casacore::Int> antenna1(thisnRow,0);
739 0 : casacore::Cube<T> gain(casacore::IPosition(3,nPar(),maxNumChan,thisnRow),T(0.0));
740 0 : casacore::Cube<casacore::Bool> solOk(nPar(),maxNumChan,thisnRow,false);
741 0 : casacore::Cube<casacore::Float> fit(1,maxNumChan,thisnRow,0.0);
742 0 : casacore::Cube<casacore::Float> fitWt(1,maxNumChan,thisnRow,0.0);
743 0 : casacore::Cube<casacore::Bool> flag(nPar(),maxNumChan,thisnRow,false);
744 0 : casacore::Cube<casacore::Float> snr(nPar(),maxNumChan,thisnRow,false);
745 :
746 0 : casacore::IPosition out(3,0,0,0); // par, chan, row
747 0 : casacore::IPosition in(4,0,0,0,0); // par, chan, ant, slot
748 0 : casacore::Int thisRow(0);
749 0 : for (casacore::Int islot = 0; islot < nTime_(iSpw); islot++) {
750 0 : in(3)=islot;
751 0 : if (thisSolOK(islot)) {
752 :
753 : // Fill slot-constant cols:
754 0 : casacore::Slice thisSlice(thisRow,nElem_);
755 0 : time(thisSlice)=thisMJDTimeStamp(islot);
756 0 : casacore::Double dt=(thisMJDStop(islot) - thisMJDStart(islot));
757 0 : if (dt<0.0) dt=1.0;
758 0 : interval(thisSlice)=dt;
759 0 : fieldId(thisSlice)=thisFieldId(islot);
760 0 : fieldName(thisSlice)=thisFieldName(islot);
761 0 : sourceName(thisSlice)=thisSourceName(islot);
762 0 : totalSolOk(thisSlice)=thisSolOK(islot);
763 0 : totalFit(thisSlice)=thisFit(islot);
764 0 : totalFitWt(thisSlice)=thisFitwt(islot);
765 :
766 : // Loop over the number of antennas
767 0 : for (casacore::Int iant = 0; iant < nElem_; iant++) {
768 0 : out(2)=thisRow;
769 0 : in(2)=iant;
770 : // Antenna index
771 0 : antenna1(thisRow)=iant;
772 :
773 0 : gain.xyPlane(thisRow)(casacore::IPosition(2,0,0),casacore::IPosition(2,nPar()-1,nChan_(iSpw)-1))=
774 : thisAntGain(casacore::IPosition(4,0,0,iant,islot),
775 0 : casacore::IPosition(4,nPar()-1,nChan_(iSpw)-1,iant,islot)).nonDegenerate(2);
776 :
777 : // Per-chan fit pars
778 0 : for (casacore::Int ichan=0; ichan<nChan_(iSpw); ichan++) {
779 : // Gain stats (slot constant, per spw?)
780 : //solOk(0,ichan,thisRow) = thisISolutionOK(iant,islot);
781 0 : fit(0,ichan,thisRow) = thisIFit(iant,islot);
782 0 : fitWt(0,ichan,thisRow) = thisIFitwt(iant,islot);
783 :
784 :
785 0 : for (casacore::Int ipar=0; ipar<nPar(); ++ipar) {
786 0 : solOk(ipar,ichan,thisRow) = parOK(iSpw)(casacore::IPosition(4,ipar,ichan,iant,islot));
787 0 : flag(ipar,ichan,thisRow) = !solOk(ipar,ichan,thisRow);
788 0 : snr(ipar,ichan,thisRow) = parSNR(iSpw)(casacore::IPosition(4,ipar,ichan,iant,islot));
789 : } // ipar
790 : } // ichan
791 :
792 : // next time round is next row
793 0 : thisRow++;
794 : }; // iant
795 :
796 :
797 : }; // thisSolOK(islot)
798 : }; // islot
799 :
800 : // Now push everything to the disk table
801 0 : calTab_->addRowMain(thisnRow);
802 : // SolvableVisJonesMCol svjmcol(*tab);
803 0 : if (svjmcol_) {delete svjmcol_; svjmcol_=NULL;}
804 0 : svjmcol_ = new SolvableCalSetMCol<T>(*calTab_);
805 :
806 0 : casacore::RefRows refRows(nMain,nMain+thisnRow-1);
807 0 : svjmcol_->time().putColumnCells(refRows,time);
808 0 : svjmcol_->timeEP().putColumnCells(refRows,timeEP);
809 0 : svjmcol_->interval().putColumnCells(refRows,interval);
810 0 : svjmcol_->antenna1().putColumnCells(refRows,antenna1);
811 0 : svjmcol_->feed1().putColumnCells(refRows,feed1);
812 0 : svjmcol_->fieldId().putColumnCells(refRows,fieldId);
813 0 : svjmcol_->arrayId().putColumnCells(refRows,arrayId);
814 0 : svjmcol_->obsId().putColumnCells(refRows,obsId);
815 0 : svjmcol_->scanNo().putColumnCells(refRows,scanNum);
816 0 : svjmcol_->processorId().putColumnCells(refRows,procId);
817 0 : svjmcol_->stateId().putColumnCells(refRows,stateId);
818 0 : svjmcol_->phaseId().putColumnCells(refRows,phaseId);
819 0 : svjmcol_->pulsarBin().putColumnCells(refRows,pulsarBin);
820 0 : svjmcol_->pulsarGateId().putColumnCells(refRows,pulsarGateId);
821 0 : svjmcol_->freqGrp().putColumnCells(refRows,freqGroup);
822 0 : svjmcol_->fieldName().putColumnCells(refRows,fieldName);
823 0 : svjmcol_->sourceName().putColumnCells(refRows,sourceName);
824 0 : svjmcol_->gain().putColumnCells(refRows,gain);
825 0 : svjmcol_->totalSolnOk().putColumnCells(refRows,totalSolOk);
826 0 : svjmcol_->totalFit().putColumnCells(refRows,totalFit);
827 0 : svjmcol_->totalFitWgt().putColumnCells(refRows,totalFitWt);
828 0 : svjmcol_->solnOk().putColumnCells(refRows,solOk);
829 0 : svjmcol_->fit().putColumnCells(refRows,fit);
830 0 : svjmcol_->fitWgt().putColumnCells(refRows,fitWt);
831 0 : svjmcol_->calDescId().putColumnCells(refRows,calDescId);
832 0 : svjmcol_->calHistoryId().putColumnCells(refRows,calHistId);
833 0 : svjmcol_->flag().putColumnCells(refRows,flag);
834 0 : svjmcol_->snr().putColumnCells(refRows,snr);
835 :
836 0 : nMain = calTab_->nRowMain();
837 :
838 : } // thisnRow > 0
839 : } // par_[iSpw]!=NULL
840 : } // iSpw
841 :
842 0 : delete calTab_;
843 :
844 0 : };
845 : /*
846 : template<class T> void CalSet<T>::store (const casacore::String& file,
847 : const casacore::String& type,
848 : const casacore::String& msname,
849 : const casacore::Bool& append)
850 : {
851 : // Write the solutions to an output calibration table
852 : // casacore::Input:
853 : // file casacore::String Cal table name
854 : // append casacore::Bool Append if true, else overwrite
855 : //
856 :
857 : // total rows to be written per Spw
858 : casacore::Vector<casacore::Int> nRow(nSpw_,0);
859 : for (casacore::Int iSpw=0;iSpw<nSpw_;iSpw++)
860 : if (solutionOK_[iSpw])
861 : nRow(iSpw)=nElem_*ntrue(*(solutionOK_[iSpw]));
862 :
863 : // Escape if nothing to write
864 : if (sum(nRow)==0)
865 : throw(casacore::AipsError("No valid calibration solutions; no table written."));
866 :
867 : // Initialization:
868 : // No. of rows in cal_main, cal_desc and cal_history
869 : casacore::Int nMain = 0;
870 : casacore::Int nDesc = 0;
871 : //casacore::Int nHist = 0;
872 :
873 : // Calibration table
874 : SolvableVisJonesTable *tab;
875 :
876 : // Open the output file if it already exists and is being appended to.
877 : if (append && casacore::Table::isWritable (file)) {
878 : tab = new SolvableVisJonesTable (file, casacore::Table::Update);
879 : nMain = tab->nRowMain();
880 : nDesc = tab->nRowDesc();
881 : //nHist = tab->nRowHistory();
882 : } else {
883 : // Create a new calibration table
884 : casacore::Table::TableOption access = casacore::Table::New;
885 : tab = new SolvableVisJonesTable (file, type, access);
886 : };
887 :
888 : // Write every spw w/ max number of channels
889 : // (eventually, CalTable should permit variable-shape cols)
890 : casacore::Int maxNumChan(1);
891 : for (casacore::Int iSpw=0; iSpw<nSpw_; iSpw++)
892 : if (par_[iSpw]!=NULL)
893 : maxNumChan=max(maxNumChan,nChan_(iSpw));
894 :
895 : // Some default values
896 : casacore::Double dzero = 0;
897 : casacore::IPosition ip(2,1,maxNumChan);
898 :
899 : // CalDesc Sub-table records
900 : CalDescRecord* descRec;
901 : casacore::Vector<casacore::Int> calDescNum(nSpw_); calDescNum=-1;
902 : for (casacore::Int iSpw=0; iSpw<nSpw_; iSpw++) {
903 :
904 : // Write a CalDesc for each spw which has solutions
905 : // Note: CalDesc index != SpwId, in general
906 :
907 : if (nRow(iSpw)>0 && par_[iSpw]!=NULL) {
908 :
909 : // Access to existing CAL_DESC columns:
910 : CalDescColumns cd(*tab);
911 :
912 : // Check if this spw already in CAL_DESC
913 : // cout << "spwCol = " << cd.spwId().getColumn() << endl;
914 :
915 : casacore::Bool newCD(true);
916 : for (casacore::Int iCD=0;iCD<nDesc;iCD++) {
917 :
918 : casacore::IPosition iCDip(1,0);
919 : if ( iSpw==(cd.spwId()(iCD))(iCDip) ) {
920 : // Don't need new CAL_DESC entry
921 : newCD=false;
922 : calDescNum(iSpw)=iCD;
923 : break;
924 : }
925 : }
926 :
927 : if (newCD) {
928 :
929 : // Cal_desc fields
930 : casacore::Vector <casacore::Int> spwId(1,iSpw);
931 : casacore::Matrix <casacore::Double> chanFreq(ip, dzero);
932 : casacore::Matrix <casacore::Double> chanWidth(ip, dzero);
933 : casacore::Array <casacore::String> polznType(ip, "");
934 : casacore::Cube <casacore::Int> chanRange(casacore::IPosition(3,2,1,maxNumChan), 0);
935 : casacore::Vector <casacore::Int> numChan(1,nChan_(iSpw));
936 : for (casacore::Int ichan=0; ichan<nChan_(iSpw); ichan++) {
937 : chanRange(0,0,ichan)=startChan_(iSpw);
938 : chanRange(1,0,ichan)=startChan_(iSpw) + nChan_(iSpw) -1;
939 : }
940 :
941 : // Fill the cal_desc record
942 : descRec = new CalDescRecord;
943 : descRec->defineNumSpw (1);
944 : descRec->defineNumChan (numChan);
945 : descRec->defineNumReceptors (2);
946 : descRec->defineNJones (2);
947 : descRec->defineSpwId (spwId);
948 : descRec->defineChanFreq (chanFreq);
949 : descRec->defineChanWidth (chanWidth);
950 : descRec->defineChanRange (chanRange);
951 : descRec->definePolznType (polznType);
952 : descRec->defineJonesType ("full");
953 : descRec->defineMSName (casacore::Path(msname).baseName());
954 :
955 : // Write the cal_desc record
956 :
957 : tab->putRowDesc (nDesc, *descRec);
958 : delete descRec;
959 :
960 : // This spw will have this calDesc index in main table
961 : calDescNum(iSpw) = nDesc;
962 : nDesc++;
963 : }
964 :
965 : }
966 :
967 : }
968 :
969 :
970 : // Now write MAIN table in column-wise fashion
971 :
972 : // Starting row in this slot
973 :
974 : for (casacore::Int iSpw=0; iSpw<nSpw_; iSpw++) {
975 :
976 : // Write table for spws which have solutions
977 : if (par_[iSpw]!=NULL) {
978 :
979 : // Create references to cal data for this spw
980 : casacore::Vector<casacore::Bool> thisSolOK; thisSolOK.reference(*(solutionOK_[iSpw]));
981 : casacore::Vector<casacore::Double> thisMJDTimeStamp; thisMJDTimeStamp.reference(*(MJDTimeStamp_[iSpw]));
982 : casacore::Vector<casacore::Double> thisMJDStart; thisMJDStart.reference(*(MJDStart_[iSpw]));
983 : casacore::Vector<casacore::Double> thisMJDStop; thisMJDStop.reference(*(MJDStop_[iSpw]));
984 : casacore::Vector<casacore::Int> thisFieldId; thisFieldId.reference(*(fieldId_[iSpw]));
985 : casacore::Vector<casacore::String> thisFieldName; thisFieldName.reference(*(fieldName_[iSpw]));
986 : casacore::Vector<casacore::String> thisSourceName; thisSourceName.reference(*(sourceName_[iSpw]));
987 : casacore::Vector<casacore::Float> thisFit; thisFit.reference(*(fit_[iSpw]));
988 : casacore::Vector<casacore::Float> thisFitwt; thisFitwt.reference(*(fitwt_[iSpw]));
989 : casacore::Array<casacore::Complex> thisAntGain; thisAntGain.reference(*(par_[iSpw]));
990 : // casacore::Matrix<casacore::Bool> thisISolutionOK; thisISolutionOK.reference(*(iSolutionOK_[iSpw]));
991 : casacore::Matrix<casacore::Float> thisIFit; thisIFit.reference(*(iFit_[iSpw]));
992 : casacore::Matrix<casacore::Float> thisIFitwt; thisIFitwt.reference(*(iFitwt_[iSpw]));
993 :
994 : casacore::Int thisnRow=nRow(iSpw);
995 :
996 : // Only if there are rows to add to table
997 : if (thisnRow > 0) {
998 :
999 : // These are constant columns (with boring values, currently)
1000 : casacore::Vector<casacore::Double> timeEP(thisnRow,0.0);
1001 : casacore::Vector<casacore::Int> feed1(thisnRow,0);
1002 : casacore::Vector<casacore::Int> arrayId(thisnRow,0);
1003 : casacore::Vector<casacore::Int> obsId(thisnRow,0);
1004 : casacore::Vector<casacore::Int> scanNum(thisnRow,0);
1005 : casacore::Vector<casacore::Int> procId(thisnRow,0);
1006 : casacore::Vector<casacore::Int> stateId(thisnRow,0);
1007 : casacore::Vector<casacore::Int> phaseId(thisnRow,0);
1008 : casacore::Vector<casacore::Int> pulsarBin(thisnRow,0);
1009 : casacore::Vector<casacore::Int> pulsarGateId(thisnRow,0);
1010 : casacore::Vector<casacore::Int> freqGroup(thisnRow,0);
1011 : casacore::Vector<casacore::Int> calHistId(thisnRow,0);
1012 :
1013 : // This is constant
1014 : casacore::Vector<casacore::Int> calDescId(thisnRow,calDescNum(iSpw));
1015 :
1016 : // These are constant per slot
1017 : // (these cols should be incremental)
1018 : casacore::Vector<casacore::Double> time(thisnRow,0.0);
1019 : casacore::Vector<casacore::Double> interval(thisnRow,0.0);
1020 : casacore::Vector<casacore::Int> fieldId(thisnRow,0);
1021 : casacore::Vector<casacore::String> fieldName(thisnRow,"");
1022 : casacore::Vector<casacore::String> sourceName(thisnRow,"");
1023 : casacore::Vector<casacore::Bool> totalSolOk(thisnRow,false);
1024 : casacore::Vector<casacore::Float> totalFit(thisnRow,0.0);
1025 : casacore::Vector<casacore::Float> totalFitWt(thisnRow,0.0);
1026 :
1027 : // These vary
1028 : casacore::Vector<casacore::Int> antenna1(thisnRow,0);
1029 : casacore::Cube<casacore::Complex> gain(casacore::IPosition(3,nPar(),maxNumChan,thisnRow),casacore::Complex(0.0,0.0));
1030 : casacore::Cube<casacore::Bool> solOk(nPar(),maxNumChan,thisnRow,false);
1031 : casacore::Cube<casacore::Float> fit(1,maxNumChan,thisnRow,0.0);
1032 : casacore::Cube<casacore::Float> fitWt(1,maxNumChan,thisnRow,0.0);
1033 : casacore::Cube<casacore::Bool> flag(nPar(),maxNumChan,thisnRow,false);
1034 : casacore::Cube<casacore::Float> snr(nPar(),maxNumChan,thisnRow,false);
1035 :
1036 : casacore::IPosition out(3,0,0,0); // par, chan, row
1037 : casacore::IPosition in(4,0,0,0,0); // par, chan, ant, slot
1038 : casacore::Int thisRow(0);
1039 : for (casacore::Int islot = 0; islot < nTime_(iSpw); islot++) {
1040 : in(3)=islot;
1041 : if (thisSolOK(islot)) {
1042 :
1043 : // Fill slot-constant cols:
1044 : casacore::Slice thisSlice(thisRow,nElem_);
1045 : time(thisSlice)=thisMJDTimeStamp(islot);
1046 : casacore::Double dt=(thisMJDStop(islot) - thisMJDStart(islot));
1047 : if (dt<0.0) dt=1.0;
1048 : interval(thisSlice)=dt;
1049 : fieldId(thisSlice)=thisFieldId(islot);
1050 : fieldName(thisSlice)=thisFieldName(islot);
1051 : sourceName(thisSlice)=thisSourceName(islot);
1052 : totalSolOk(thisSlice)=thisSolOK(islot);
1053 : totalFit(thisSlice)=thisFit(islot);
1054 : totalFitWt(thisSlice)=thisFitwt(islot);
1055 :
1056 : // Loop over the number of antennas
1057 : for (casacore::Int iant = 0; iant < nElem_; iant++) {
1058 : out(2)=thisRow;
1059 : in(2)=iant;
1060 : // Antenna index
1061 : antenna1(thisRow)=iant;
1062 :
1063 :
1064 : gain.xyPlane(thisRow)(casacore::IPosition(2,0,0),casacore::IPosition(2,nPar()-1,nChan_(iSpw)-1))=
1065 : thisAntGain(casacore::IPosition(4,0,0,iant,islot),casacore::IPosition(4,nPar()-1,nChan_(iSpw)-1,iant,islot)).nonDegenerate(2);
1066 :
1067 :
1068 : // Per-chan fit pars
1069 : for (casacore::Int ichan=0; ichan<nChan_(iSpw); ichan++) {
1070 : // Gain stats (slot constant, per spw?)
1071 : //solOk(0,ichan,thisRow) = thisISolutionOK(iant,islot);
1072 : fit(0,ichan,thisRow) = thisIFit(iant,islot);
1073 : fitWt(0,ichan,thisRow) = thisIFitwt(iant,islot);
1074 :
1075 :
1076 : for (casacore::Int ipar=0; ipar<nPar(); ++ipar) {
1077 : solOk(ipar,ichan,thisRow) = parOK(iSpw)(casacore::IPosition(4,ipar,ichan,iant,islot));
1078 : flag(ipar,ichan,thisRow) = !solOk(ipar,ichan,thisRow);
1079 : snr(ipar,ichan,thisRow) = parSNR(iSpw)(casacore::IPosition(4,ipar,ichan,iant,islot));
1080 : }
1081 : }
1082 :
1083 : // next time round is next row
1084 : thisRow++;
1085 : };
1086 :
1087 : };
1088 : };
1089 :
1090 : // Now push everything to the disk table
1091 : tab->addRowMain(thisnRow);
1092 : SolvableVisJonesMCol svjmcol(*tab);
1093 :
1094 : casacore::RefRows refRows(nMain,nMain+thisnRow-1);
1095 : svjmcol.time().putColumnCells(refRows,time);
1096 : svjmcol.timeEP().putColumnCells(refRows,timeEP);
1097 : svjmcol.interval().putColumnCells(refRows,interval);
1098 : svjmcol.antenna1().putColumnCells(refRows,antenna1);
1099 : svjmcol.feed1().putColumnCells(refRows,feed1);
1100 : svjmcol.fieldId().putColumnCells(refRows,fieldId);
1101 : svjmcol.arrayId().putColumnCells(refRows,arrayId);
1102 : svjmcol.obsId().putColumnCells(refRows,obsId);
1103 : svjmcol.scanNo().putColumnCells(refRows,scanNum);
1104 : svjmcol.processorId().putColumnCells(refRows,procId);
1105 : svjmcol.stateId().putColumnCells(refRows,stateId);
1106 : svjmcol.phaseId().putColumnCells(refRows,phaseId);
1107 : svjmcol.pulsarBin().putColumnCells(refRows,pulsarBin);
1108 : svjmcol.pulsarGateId().putColumnCells(refRows,pulsarGateId);
1109 : svjmcol.freqGrp().putColumnCells(refRows,freqGroup);
1110 : svjmcol.fieldName().putColumnCells(refRows,fieldName);
1111 : svjmcol.sourceName().putColumnCells(refRows,sourceName);
1112 : svjmcol.gain().putColumnCells(refRows,gain);
1113 : svjmcol.totalSolnOk().putColumnCells(refRows,totalSolOk);
1114 : svjmcol.totalFit().putColumnCells(refRows,totalFit);
1115 : svjmcol.totalFitWgt().putColumnCells(refRows,totalFitWt);
1116 : svjmcol.solnOk().putColumnCells(refRows,solOk);
1117 : svjmcol.fit().putColumnCells(refRows,fit);
1118 : svjmcol.fitWgt().putColumnCells(refRows,fitWt);
1119 : svjmcol.calDescId().putColumnCells(refRows,calDescId);
1120 : svjmcol.calHistoryId().putColumnCells(refRows,calHistId);
1121 : svjmcol.flag().putColumnCells(refRows,flag);
1122 : svjmcol.snr().putColumnCells(refRows,snr);
1123 :
1124 :
1125 : nMain = tab->nRowMain();
1126 :
1127 : }
1128 : }
1129 : }
1130 :
1131 : delete tab;
1132 :
1133 : };
1134 : */
1135 : } //# NAMESPACE CASA - END
1136 :
|
