Commits
1423 1423 | } |
1424 1424 | } |
1425 1425 | gsl_matrix_free(lu); |
1426 1426 | gsl_matrix_free(inv_hess); |
1427 1427 | // SNR[i], according to aips, is 1/sqrt(sigma2*hess(i1,i1)*0.5); |
1428 1428 | // Note that in AIPS i1 ranges over 1..NANT |
1429 1429 | // We use 1 as a success code. |
1430 1430 | return 1; |
1431 1431 | } |
1432 1432 | |
1433 + | |
1434 + | Int |
1435 + | findRefAntWithData(SDBList& sdbs, Vector<Int>& refAntList, Int prtlev) { |
1436 + | std::set<Int> activeAntennas; |
1437 + | for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) { |
1438 + | SolveDataBuffer& s(sdbs(ibuf)); |
1439 + | if (!s.Ok()) |
1440 + | continue; |
1441 + | Cube<Bool> fl = s.flagCube(); |
1442 + | for (Int irow=0; irow!=s.nRows(); irow++) { |
1443 + | if (s.flagRow()(irow)) |
1444 + | continue; |
1445 + | Int a1(s.antenna1()(irow)); |
1446 + | Int a2(s.antenna2()(irow)); |
1447 + | // Not using irow |
1448 + | Matrix<Bool> flr = fl.xyPlane(irow); |
1449 + | if (!allTrue(flr)) { |
1450 + | activeAntennas.insert(a1); |
1451 + | activeAntennas.insert(a2); |
1452 + | } |
1453 + | } |
1454 + | } |
1455 + | if (prtlev > 2) { |
1456 + | cout << "[FringeJones.cc::findRefAntWithData] activeAntennas: "; |
1457 + | std::copy( |
1458 + | activeAntennas.begin(), |
1459 + | activeAntennas.end(), |
1460 + | std::ostream_iterator<Int>(std::cout, " ") |
1461 + | ); |
1462 + | cout << endl; |
1463 + | } |
1464 + | Int refAnt = -1; |
1465 + | for (Vector<Int>::ConstIteratorSTL a = refAntList.begin(); a != refAntList.end(); a++) { |
1466 + | if (activeAntennas.find(*a) != activeAntennas.end()) { |
1467 + | if (prtlev > 2) |
1468 + | cout << "[FringeJones.cc::findRefAntWithData] We are choosing refant " << *a << endl; |
1469 + | refAnt = *a; |
1470 + | break; |
1471 + | } else { |
1472 + | if (prtlev > 2) |
1473 + | cout << "[FringeJones.cc::findRefAntWithData] No data for refant " << *a << endl; |
1474 + | } |
1475 + | } |
1476 + | return refAntList(0); |
1477 + | } |
1478 + | |
1479 + | |
1433 1480 | // Stolen from SolveDataBuffer |
1434 1481 | void |
1435 1482 | aggregateTimeCentroid(SDBList& sdbs, Int refAnt, std::map<Int, Double>& aggregateTime) { |
1436 1483 | // Weighted average of SDBs' timeCentroids |
1437 1484 | std::map<Int, Double> aggregateWeight; |
1438 1485 | for (Int i=0; i < sdbs.nSDB(); ++i) { |
1439 1486 | SolveDataBuffer& s = sdbs(i); |
1440 1487 | Vector<Double> wtvD; |
1441 1488 | // Sum over correlations and channels to get a vector of weights for each row |
1442 1489 | Vector<Float> wtv(partialSums(s.weightSpectrum(), IPosition(2, 0, 1))); |
1886 1933 | for (uInt ispw=0;ispw<spwMap().nelements();++ispw) |
1887 1934 | if (spwMap()(ispw)>-1) |
1888 1935 | KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw)); |
1889 1936 | } |
1890 1937 | |
1891 1938 | } |
1892 1939 | |
1893 1940 | void FringeJones::setSolve(const Record& solve) { |
1894 1941 | |
1895 1942 | // Call parent to do conventional things |
1943 + | if (prtlev() > 2) { |
1944 + | cout << "Before GJones::setSolve" << endl |
1945 + | << "FringeJones::setSolve()" <<endl |
1946 + | << "FringeJones::refant() = "<< refant() <<endl |
1947 + | << "FringeJones::refantlist() = "<< refantlist() <<endl; |
1948 + | } |
1896 1949 | GJones::setSolve(solve); |
1897 1950 | |
1898 1951 | // if (!ct_) |
1899 1952 | // throw(AipsError("No calibration table specified")); |
1900 1953 | // cerr << "setSolve here, ct_: "<< ct_ << endl; |
1901 1954 | |
1902 1955 | // Trap unspecified refant: |
1903 1956 | if (refant()<0) |
1904 1957 | throw(AipsError("Please specify a good reference antenna (refant) explicitly.")); |
1905 1958 | if (solve.isDefined("zerorates")) { |
2037 2090 | Double t0 = sdbs(0).time()(0); |
2038 2091 | Double dt0 = refTime() - t0; |
2039 2092 | //Double df0 = ref_freq - sdbs.freqs()(0); |
2040 2093 | //Double df0 = 0; |
2041 2094 | Double df0 = centroidFreq - sdbs(0).freqs()(0); // global center to global edge |
2042 2095 | |
2043 2096 | logSink() << "Solving for fringes for spw=" << currSpw() << " at t=" |
2044 2097 | << MVTime(refTime()/C::day).string(MVTime::YMD,7) << LogIO::POST; |
2045 2098 | |
2046 2099 | std::map<Int, Double> aggregateTime; |
2100 + | // Set the refant to the first choice that has data! |
2101 + | refant() = findRefAntWithData(sdbs, refantlist(), prtlev()); |
2102 + | if (refant()<0) |
2103 + | throw(AipsError("No valid reference antenna supplied.")); |
2104 + | |
2047 2105 | aggregateTimeCentroid(sdbs, refant(), aggregateTime); |
2048 2106 | |
2049 2107 | if (DEVDEBUG) { |
2050 2108 | std::cerr << "Weighted time centroids" << endl; |
2051 2109 | for (auto it=aggregateTime.begin(); it!=aggregateTime.end(); ++it) |
2052 2110 | std::cerr << it->first << " => " << it->second - t0 << std::endl; |
2053 2111 | } |
2054 2112 | |
2055 2113 | DelayRateFFT drf(sdbs, refant(), delayWindow(), rateWindow()); |
2056 2114 | drf.FFT(); |