diff --git a/CEP/DP3/DPPP/include/DPPP/DPBuffer.h b/CEP/DP3/DPPP/include/DPPP/DPBuffer.h
index df9d330b52f11f6fa16ecc201d9b9a71f7ccb25f..d8c3a1cd870388ed5b0ba4c7cfa010bb3b357c0d 100644
--- a/CEP/DP3/DPPP/include/DPPP/DPBuffer.h
+++ b/CEP/DP3/DPPP/include/DPPP/DPBuffer.h
@@ -165,7 +165,7 @@ namespace LOFAR {
{ return itsTime; }
// Get or set the row numbers used by the DPInput class.
- // It can be empty (e.g. when MSReader inserted an dummy time slot).
+ // It can be empty (e.g. when MSReader inserted a dummy time slot).
void setRowNrs (const casa::Vector<uint>& rownrs)
{ itsRowNrs.reference (rownrs); }
const casa::Vector<uint>& getRowNrs() const
diff --git a/CEP/DP3/DPPP/include/DPPP/DPInput.h b/CEP/DP3/DPPP/include/DPPP/DPInput.h
index c40f49b148c52e9f599590c9098b94c3108fc05c..58e4f816b28e269593d820cf95e585b3c8007604 100644
--- a/CEP/DP3/DPPP/include/DPPP/DPInput.h
+++ b/CEP/DP3/DPPP/include/DPPP/DPInput.h
@@ -98,6 +98,10 @@ namespace LOFAR {
const casa::Vector<int>& getAnt2() const
{ return itsAnt2; }
+ // Get the baseline table index of the autocorrelations.
+ // A negative value means there are no autocorrelations for that antenna.
+ const vector<int>& getAutoCorrIndex() const;
+
// Get the lengths of the baselines (in meters).
const vector<double>& getBaselineLengths() const;
@@ -153,14 +157,16 @@ namespace LOFAR {
uint itsNrChan;
uint itsNrCorr;
uint itsNrBl;
- casa::Vector<int> itsAnt1; //# ant1 of all baselines
- casa::Vector<int> itsAnt2; //# ant2 of all baselines
- mutable vector<double> itsBLength; //# baseline lengths
+ casa::Vector<casa::Int> itsAnt1; //# ant1 of all baselines
+ casa::Vector<casa::Int> itsAnt2; //# ant2 of all baselines
+ mutable vector<double> itsBLength; //# baseline lengths
+ mutable vector<int> itsAutoCorrIndex; //# autocorr index per ant
casa::Vector<casa::String> itsAntNames;
vector<casa::MPosition> itsAntPos;
casa::MPosition itsArrayPos;
casa::MDirection itsPhaseCenter;
casa::Vector<double> itsChanFreqs;
+ casa::Vector<double> itsChanWidths;
};
} //# end namespace
diff --git a/CEP/DP3/DPPP/include/DPPP/MSReader.h b/CEP/DP3/DPPP/include/DPPP/MSReader.h
index ad1bf4c58ed9c0dc54da77fa7b62d0bbf099b9a3..513f5a0c91b463a8c0b4afcb55209d47fb5dd488 100644
--- a/CEP/DP3/DPPP/include/DPPP/MSReader.h
+++ b/CEP/DP3/DPPP/include/DPPP/MSReader.h
@@ -52,6 +52,7 @@ namespace LOFAR {
// Currently the following can be given:
// <ul>
// <li> msin: name of the MS
+ // <li> msin.autoweight: calculate weights from autocorrelations? [no]
// <li> msin.startchan: first channel to use [0]
// <li> msin.nchan: number of channels to use [all]
// <li> msin.useflag: use the existing flags? [yes]
@@ -101,6 +102,8 @@ namespace LOFAR {
// WEIGHT_SPECTRUM is used if present and containing valid data,
// otherwise column WEIGHT is used. The weights of an inserted
// time slot are set to 0.
+ // If autoweight is on, the autocorrelations are used to
+ // calculate proper weights.
// </td>
// </tr>
// <tr>
@@ -211,10 +214,14 @@ namespace LOFAR {
// Calculate the UVWs for a missing time slot.
void calcUVW();
+ // Calculate the weights from the autocorrelations.
+ void autoWeight (casa::Cube<float>& weights);
+
//# Data members.
casa::Table itsMS;
casa::TableIterator itsIter;
casa::String itsDataColName;
+ bool itsAutoWeight; //# calculate weights from autocorr?
bool itsHasWeightSpectrum;
bool itsUseFlags;
bool itsUseAllChan; //# all channels (i.e. no slicer)?
diff --git a/CEP/DP3/DPPP/src/Averager.cc b/CEP/DP3/DPPP/src/Averager.cc
index d639fd8f7cdd0df68a36e688e238cc3f1b7962c0..1dba3f824238fff0c6288876c7b7426f1e3fc006 100644
--- a/CEP/DP3/DPPP/src/Averager.cc
+++ b/CEP/DP3/DPPP/src/Averager.cc
@@ -213,7 +213,7 @@ namespace LOFAR {
float sumw=0;
uint np = 0;
for (uint j=0; j<nch; ++j) {
- sumd += indata[inxi];
+ sumd += indata[inxi]; // Note: weight is accounted for in process
sumw += inwght[inxi];
np += innp[inxi];
inxi += ncorr;
diff --git a/CEP/DP3/DPPP/src/DPInput.cc b/CEP/DP3/DPPP/src/DPInput.cc
index a3e673f721e89bd83c9750151f1ccbcd71a08564..d01bb472847dd14d59e22521f60b377601879eb4 100644
--- a/CEP/DP3/DPPP/src/DPInput.cc
+++ b/CEP/DP3/DPPP/src/DPInput.cc
@@ -78,6 +78,22 @@ namespace LOFAR {
return itsBLength;
}
+ const vector<int>& DPInput::getAutoCorrIndex() const
+ {
+ if (itsAutoCorrIndex.empty()) {
+ int nant = 1 + std::max(max(itsAnt1), max(itsAnt2));
+ itsAutoCorrIndex.resize (nant);
+ std::fill (itsAutoCorrIndex.begin(), itsAutoCorrIndex.end(), -1);
+ // Keep the baseline table index for the autocorrelations.
+ for (uint i=0; i<itsAnt1.size(); ++i) {
+ if (itsAnt1[i] == itsAnt2[i]) {
+ itsAutoCorrIndex[itsAnt1[i]] = i;
+ }
+ }
+ }
+ return itsAutoCorrIndex;
+ }
+
Cube<bool> DPInput::fetchFullResFlags (const DPBuffer& buf,
const RefRows& rowNrs,
bool merge)
diff --git a/CEP/DP3/DPPP/src/MSReader.cc b/CEP/DP3/DPPP/src/MSReader.cc
index a319dd17596af11da5d1eddf12bd0d8d2c83c963..d7e3235f581a75a4d58c40e9e37b867a16a9b3bc 100644
--- a/CEP/DP3/DPPP/src/MSReader.cc
+++ b/CEP/DP3/DPPP/src/MSReader.cc
@@ -58,6 +58,7 @@ namespace LOFAR {
string endTimeStr = parset.getString (prefix+"endtime", "");
itsUseFlags = parset.getBool (prefix+"useflag", true);
itsDataColName = parset.getString (prefix+"datacolumn", "DATA");
+ itsAutoWeight = parset.getBool (prefix+"autoweight", False);
// Prepare the MS access and get time info.
double startTime, endTime;
prepare (startTime, endTime, itsInterval);
@@ -262,6 +263,7 @@ namespace LOFAR {
os << " ntimes: " << itsMS.nrow() / itsNrBl << std::endl;
os << " time interval: " << itsInterval << std::endl;
os << " DATA column: " << itsDataColName << std::endl;
+ os << " autoweight: " << itsAutoWeight << std::endl;
}
void MSReader::showCounts (std::ostream& os) const
@@ -342,8 +344,8 @@ namespace LOFAR {
ASSERTSTR (sortab.nrow() == itsNrBl,
"The MS appears to have multiple subbands");
// Get the baselines.
- ROScalarColumn<int>(itsIter.table(), "ANTENNA1").getColumn (itsAnt1);
- ROScalarColumn<int>(itsIter.table(), "ANTENNA2").getColumn (itsAnt2);
+ ROScalarColumn<Int>(itsIter.table(), "ANTENNA1").getColumn (itsAnt1);
+ ROScalarColumn<Int>(itsIter.table(), "ANTENNA2").getColumn (itsAnt2);
// Keep the row numbers of the first part to be used for the meta info
// of possibly missing time slots.
itsBaseRowNrs = itsIter.table().rowNumbers(itsMS);
@@ -367,8 +369,10 @@ namespace LOFAR {
// Read the center frequencies of all channels.
Table spwtab(itsMS.keywordSet().asTable("SPECTRAL_WINDOW"));
ROArrayColumn<double> freqCol (spwtab, "CHAN_FREQ");
+ ROArrayColumn<double> widthCol (spwtab, "CHAN_WIDTH");
// Take only the channels used in the input.
- itsChanFreqs = freqCol(0);
+ itsChanFreqs = freqCol(0);
+ itsChanWidths = widthCol(0);
// Get the array position using the telescope name from the OBSERVATION
// subtable.
Table obstab (itsMS.keywordSet().asTable ("OBSERVATION"));
@@ -446,35 +450,79 @@ namespace LOFAR {
weights = 0;
return weights;
}
- // Get weights for entire spectrum if pesent.
+ Cube<float> weights;
+ // Get weights for entire spectrum if present.
if (itsHasWeightSpectrum) {
ROArrayColumn<float> wsCol(itsMS, "WEIGHT_SPECTRUM");
// Using getColumnCells(rowNrs,itsColSlicer) fails for LofarStMan.
// Hence work around it.
- Cube<float> weights = wsCol.getColumnCells (rowNrs);
- return (itsUseAllChan ? weights : weights(itsArrSlicer));
- }
- // No spectrum present, so get global weights and assign to each channel.
- ROArrayColumn<float> wCol(itsMS, "WEIGHT");
- Matrix<float> inArr = wCol.getColumnCells (rowNrs);
- Cube<float> outArr(itsNrCorr, itsNrChan, itsNrBl);
- float* inPtr = inArr.data();
- float* outPtr = outArr.data();
- for (uint i=0; i<itsNrBl; ++i) {
- // If global weights are zero, set them to 1. Some old MSs need that.
- for (uint k=0; k<itsNrCorr; ++k) {
- if (inPtr[k] == 0.) {
- inPtr[k] = 1.;
- }
+ weights.reference (wsCol.getColumnCells (rowNrs));
+ if (!itsUseAllChan) {
+ weights.reference (weights(itsArrSlicer));
}
- for (uint j=0; j<itsNrChan; ++j) {
+ } else {
+ // No spectrum present; get global weights and assign to each channel.
+ ROArrayColumn<float> wCol(itsMS, "WEIGHT");
+ Matrix<float> inArr = wCol.getColumnCells (rowNrs);
+ Cube<float> outArr(itsNrCorr, itsNrChan, itsNrBl);
+ float* inPtr = inArr.data();
+ float* outPtr = outArr.data();
+ for (uint i=0; i<itsNrBl; ++i) {
+ // If global weights are zero, set them to 1. Some old MSs need that.
for (uint k=0; k<itsNrCorr; ++k) {
- *outPtr++ = inPtr[k];
+ if (inPtr[k] == 0.) {
+ inPtr[k] = 1.;
+ }
+ }
+ for (uint j=0; j<itsNrChan; ++j) {
+ for (uint k=0; k<itsNrCorr; ++k) {
+ *outPtr++ = inPtr[k];
+ }
+ }
+ inPtr += itsNrCorr;
+ }
+ weights.reference (outArr);
+ }
+ if (itsAutoWeight) {
+ // Adapt weights using autocorrelations.
+ autoWeight (weights);
+ }
+ return weights;
+ }
+
+ void MSReader::autoWeight (Cube<float>& weights)
+ {
+ const double* chanWidths = itsChanWidths.data();
+ uint npol = weights.shape()[0];
+ uint nchan = weights.shape()[1];
+ uint nbl = weights.shape()[2];
+ // Get the autocorrelations indices.
+ const vector<int>& autoInx = getAutoCorrIndex();
+ // Calculate the weight for each cross-correlation data point.
+ const Complex* data = itsBuffer.getData().data();
+ float* weight = weights.data();
+ for (uint i=0; i<nbl; ++i) {
+ // Can only be done if both autocorrelations are present.
+ if (itsAnt1[i] != itsAnt2[i] &&
+ autoInx[itsAnt1[i]] >= 0 && autoInx[itsAnt2[i]] >= 0) {
+ // Get offset of both autocorr in data array.
+ const Complex* auto1 = data + autoInx[itsAnt1[i]]*nchan*npol;
+ const Complex* auto2 = data + autoInx[itsAnt2[i]]*nchan*npol;
+ for (uint j=0; j<nchan; ++j) {
+ double w = chanWidths[j] * itsInterval;
+ *weight++ *= w / (auto1[0].real() * auto2[0].real()); // XX
+ if (npol == 4) {
+ *weight++ *= w / (auto1[0].real() * auto2[1].real()); // XY
+ *weight++ *= w / (auto1[1].real() * auto2[0].real()); // YX
+ *weight++ *= w / (auto1[1].real() * auto2[1].real()); // YY
+ } else if (npol == 2) {
+ *weight++ *= w / (auto1[1].real() * auto2[1].real()); // YY
+ }
}
+ } else {
+ weight += nchan*npol;
}
- inPtr += itsNrCorr;
}
- return outArr;
}
Cube<bool> MSReader::getFullResFlags (const RefRows& rowNrs)
diff --git a/CEP/DP3/DPPP/test/tAverager.cc b/CEP/DP3/DPPP/test/tAverager.cc
index 3329568362eab81fab1ab5a8add862b35a4697e1..9edf29153d1543926fb8bc5224d04ffedcafced0 100644
--- a/CEP/DP3/DPPP/test/tAverager.cc
+++ b/CEP/DP3/DPPP/test/tAverager.cc
@@ -117,9 +117,7 @@ private:
for (int j=itsCount*itsNAvgTime; j<itsCount*itsNAvgTime+navgtime; ++j) {
for (int i=0; i<int(data.size()); ++i) {
data.data()[i] += Complex(i+j*10,i-1000+j*6);
- if (!itsFlag) {
- weights.data()[i] += float(1);
- }
+ weights.data()[i] += float(1);
}
}
fullResFlags(Slicer(IPosition(3,0,0,0),