diff --git a/CEP/DP3/DPPP/include/DPPP/GainCal.h b/CEP/DP3/DPPP/include/DPPP/GainCal.h index 0dd60d5c2ea58c2fdd438c3f60041cdd446800be..78bcf6fac682f9923b31c74f3061b43d8a90121e 100644 --- a/CEP/DP3/DPPP/include/DPPP/GainCal.h +++ b/CEP/DP3/DPPP/include/DPPP/GainCal.h @@ -95,14 +95,6 @@ namespace LOFAR { // Invert solution (for applying it) casa::Cube<casa::DComplex> invertSol(const casa::Cube<casa::DComplex>& sol); - // Counts the number of antennas with non-flagged data, - // Set a map for the used antennas in iS, returns the number of antennas - void setAntennaMaps (const casa::Bool* flag, uint freqCell); - - // Remove rows and colums corresponding to antennas with too much - // flagged data from vis and mvis - void removeDeadAntennas (); - // Fills the matrices itsVis and itsMVis void fillMatrices (casa::Complex* model, casa::Complex* data, float* weight, const casa::Bool* flag); @@ -153,11 +145,10 @@ namespace LOFAR { uint itsMaxIter; double itsTolerance; - bool itsPropagateSolutions; // Not used currently, TODO: use this + bool itsPropagateSolutions; uint itsSolInt; // Time cell size uint itsNChan; // Frequency cell size uint itsNFreqCells; - uint itsMinBLperAnt; uint itsTimeSlotsPerParmUpdate; uint itsConverged; diff --git a/CEP/DP3/DPPP/src/GainCal.cc b/CEP/DP3/DPPP/src/GainCal.cc index 143beefe32b6a9068231cdf92d855a5aa032f968..a50fdf948636a9d0324695b79f155620ce3bf5f0 100644 --- a/CEP/DP3/DPPP/src/GainCal.cc +++ b/CEP/DP3/DPPP/src/GainCal.cc @@ -83,7 +83,6 @@ namespace LOFAR { itsSolInt (parset.getInt(prefix + "solint", 1)), itsNChan (parset.getInt(prefix + "nchan", 0)), itsNFreqCells (0), - itsMinBLperAnt (parset.getInt(prefix + "minblperant", 4)), itsTimeSlotsPerParmUpdate (parset.getInt(prefix + "timeslotsperparmupdate", 500)), itsConverged (0), @@ -329,7 +328,7 @@ namespace LOFAR { // Start new solution interval for (uint freqCell=0; freqCell<itsNFreqCells; freqCell++) { - setAntennaMaps(flag, freqCell); + iS[freqCell].clearStationFlagged(); iS[freqCell].resetVis(); } } @@ -475,58 +474,6 @@ namespace LOFAR { } } - void GainCal::setAntennaMaps (const Bool* flag, uint freqCell) { - uint nCr=info().ncorr(); - uint nCh=info().nchan(); - uint nBl=info().nbaselines(); - - iS[freqCell].clearStationFlagged(); - - casa::Vector<casa::uInt> dataPerAntenna; // nAnt - dataPerAntenna.resize(info().antennaNames().size()); - - // TODO: implement smarter graph algorithm here that requires - // less than O(n_ant^3) worst case operations. - bool flaggedAntsAdded=true; - uint loopcount=0; - while (flaggedAntsAdded) { - dataPerAntenna=0; - for (uint bl=0;bl<nBl;++bl) { - uint ant1=info().getAnt1()[bl]; - uint ant2=info().getAnt2()[bl]; - if (ant1==ant2) { - continue; - } - uint chmax=min((freqCell+1)*itsNChan, nCh); - for (uint ch=freqCell*itsNChan;ch<chmax;++ch) { - for (uint cr=0;cr<nCr;++cr) { - if (!(flag[bl*nCr*nCh + ch*nCr + cr] - || iS[freqCell].getStationFlagged()[ant1] - || iS[freqCell].getStationFlagged()[ant2] - )) { - dataPerAntenna(ant1)++; - dataPerAntenna(ant2)++; - } - } - } - } - - flaggedAntsAdded=false; - for (uint ant=0; ant<info().antennaNames().size(); ++ant) { - if (dataPerAntenna(ant)>nCr*itsMinBLperAnt) { - iS[freqCell].getStationFlagged()[ant]=false; // Index in stefcal numbering - } else { // Not enough data - //cout<<"flagging station "<<ant<<", "<<dataPerAntenna(ant)<<endl; - if (!iS[freqCell].getStationFlagged()[ant]) { - iS[freqCell].getStationFlagged()[ant]=true; - flaggedAntsAdded=true; - } - } - } - loopcount++; - } - } - void GainCal::stefcal () { itsTimerSolve.start(); diff --git a/CEP/DP3/DPPP/src/StefCal.cc b/CEP/DP3/DPPP/src/StefCal.cc index be30fc93c0008a06cac96f72ab998813fcceafd6..852cd8acf80574f1d8ffedb195aed1dd61004a3a 100644 --- a/CEP/DP3/DPPP/src/StefCal.cc +++ b/CEP/DP3/DPPP/src/StefCal.cc @@ -285,7 +285,13 @@ namespace LOFAR { vis_p=&_vis(IPosition(6,0,1,time,ch,1,st1)); for (uint st2=0;st2<_nSt;++st2) { t(3) += conj(_z(st2+_nSt*ch+_nSt*_nChan*time,3)) * vis_p[st2]; }// itsVis(IPosition(6,st2,1,time,ch,1,st1)) } } - DComplex invdet= 1./(w(0) * w (3) - w(1)*w(2)); + DComplex invdet= w(0) * w (3) - w(1)*w(2); + if (abs(invdet)==0) { + _stationFlagged[st1] = true; + _g(st1,0) = 0; + continue; + } + invdet= 1./invdet; _g(st1,0) = invdet * ( w(3) * t(0) - w(1) * t(2) ); _g(st1,1) = invdet * ( w(3) * t(1) - w(1) * t(3) ); _g(st1,2) = invdet * ( w(0) * t(2) - w(2) * t(0) ); @@ -333,6 +339,7 @@ namespace LOFAR { //cout<<", w="<<ww<<" "; if (ww==0) { _stationFlagged[st1%_nSt]=true; + _g(st1,0)=0; continue; } _g(st1,0)=tt/ww;