diff --git a/Appl/CEP/CS1/CS1_pp_lib/scripts/CS1_Offline_pipeline_control.py b/Appl/CEP/CS1/CS1_pp_lib/scripts/CS1_Offline_pipeline_control.py
new file mode 100644
index 0000000000000000000000000000000000000000..6f630b8c56d41cafecc3fec32543d01c90a65d17
--- /dev/null
+++ b/Appl/CEP/CS1/CS1_pp_lib/scripts/CS1_Offline_pipeline_control.py
@@ -0,0 +1,127 @@
+from optparse import OptionParser
+from LOFAR_Parset import Parset
+import sys, os, time
+
+parser = OptionParser(usage='%prog [options] -m<Obs name>')
+parser.add_option("-o", "--obs", dest="obs", default="",
+ help="Observation name")
+
+options, args = parser.parse_args()
+
+Observation = options.obs
+
+##settings
+node_dir = '/local/renting'
+run_id = 'test01'
+nodes = [('lioff001','SB0.MS'),('lioff002','SB1.MS'),('lioff003','SB2.MS'),
+('lioff004','SB3.MS'),('lioff005','SB4.MS'),('lioff006','SB5.MS'),
+('lioff007','SB6.MS'),('lioff008','SB7.MS'),('lioff009','SB8.MS'),
+('lioff010','SB9.MS'),('lioff011','SB10.MS'),('lioff012','SB11.MS'),
+('lioff013','SB12.MS'),('lioff014','SB13.MS'),('lioff015','SB14.MS'),
+('lioff016','SB15.MS'),('lioff017','SB16.MS'),('lioff018','SB17.MS')
+]
+
+
+if Observation == "":
+ print "No Observation set given, exiting"
+ sys.exit(1)
+else:
+ if not os.path.exists(Observation):
+ print "Observation does not exist, exiting"
+ sys.exit(2)
+
+try:
+ os.environ['AIPSPATH']
+except:
+ print "AIPSPATH not set, exiting"
+ sys.exit(3)
+
+try:
+ os.environ['SSH_AGENT_PID']
+except:
+ print "SSH-AGENT not active, exiting"
+ sys.exit(3)
+
+pidlist = []
+
+def spawn(host, MS):
+ try:
+ pid = os.fork()
+ except OSError, e:
+ print 'Unable to fork:' + str(e)
+ os._exit(1)
+ if pid == 0: ## the new client
+ ##Copy node script so we are sure we are running the right version
+ os.system('scp CS1_Offline_pipeline_node.py ' + host + ':' + node_dir +
+ '/CS1_Offline_pipeline_node.py')
+ ##Copy the subband MS
+ os.system('scp -r ' + Observation + '/' + MS +
+ ' ' + host + ':' + node_dir + '/' + MS)
+ ##Run the node script on every node
+ os.system('ssh -t ' + host +
+ ' "setenv PYTHONPATH /app/LOFAR/stable;' +
+ 'source /app/scripts/doUnstableAIPS++;cd '+node_dir + ';python ' + node_dir +
+ '/CS1_Offline_pipeline_node.py -m' +
+ node_dir + '/' + MS +
+ ' >> ' + run_id + '.log"')
+ os._exit(0)
+ else: ## parent process
+ pidlist.append(pid)
+
+for node in nodes:
+ spawn(node[0], node[1])
+
+time.sleep(10) ## wait for childs to spawn
+for p in pidlist: ## wait for the clients to finish
+ (pid, status) = os.waitpid(p, 0)
+ print 'result: ' + str(nodes[pidlist.index(pid)]) + ' ' + str(pid) + ' ' + str(status)
+
+MSlist = []
+for node in nodes:
+ ## get the squashed measurementsets back
+ os.system('scp -r ' + node[0] + ':' + node_dir + '/s' + node[1] + ' .')
+ MSlist.append('s' + node[1])
+
+MS = 'result.MS'
+
+Combiner_parset = Parset()
+Combiner_parset['inms'] = MSlist
+Combiner_parset['out'] = MS
+Combiner_parset.writeToFile("CS1_SPWCombine.parset")
+fd = open("CS1_SPWComine.debug", 'w')
+fd.write("Global 20\n")
+fd.write("Everything 20\n")
+fd.close()
+
+## Unsure if these are usefull imager parameters, probably onlt around 60 MHz
+Imager_parset = Parset()
+Imager_parset['ms'] = MS
+Imager_parset['compress'] = "False"
+Imager_parset['datamode'] = "channel"
+Imager_parset['imagemode'] = "mfs"
+Imager_parset['spwid'] = [0,1]
+Imager_parset['nchan'] = 256
+Imager_parset['start'] = 0
+Imager_parset['step'] = 1
+Imager_parset['nx'] = 512
+Imager_parset['ny'] = 512
+Imager_parset['cellx'] = 750
+Imager_parset['celly'] = 750
+Imager_parset['stokes'] = "I"
+Imager_parset['weighttype'] = "natural"
+Imager_parset['weightnpixels'] = 1024
+Imager_parset['tile'] = 32
+Imager_parset['padding'] = 1.0
+Imager_parset['gridfunction'] = "SF"
+Imager_parset['imagetype'] = "observed"
+Imager_parset['imagename'] = "observed.image"
+Flagger_parset.writeToFile("CS1_Imager.parset")
+fd = open("CS1_Imager.debug", 'w')
+fd.write("Global 20\n")
+fd.write("Everything 20\n")
+fd.close()
+
+## some processing on the final MS
+os.system("/app/LOFAR/stable/CS1_SPWCombine")
+os.system("glish -l /app/LOFAR/stable/flag_auto.g " + MS)
+#os.system("/app/LOFAR/stable/CS1_Imager")
diff --git a/Appl/CEP/CS1/CS1_pp_lib/scripts/CS1_Offline_pipeline_node.py b/Appl/CEP/CS1/CS1_pp_lib/scripts/CS1_Offline_pipeline_node.py
new file mode 100644
index 0000000000000000000000000000000000000000..529a3a960c992a57e7604dd0510803ec24a6f49e
--- /dev/null
+++ b/Appl/CEP/CS1/CS1_pp_lib/scripts/CS1_Offline_pipeline_node.py
@@ -0,0 +1,65 @@
+from optparse import OptionParser
+from LOFAR_Parset import Parset
+import sys, os
+
+parser = OptionParser(usage='%prog [options] -m<MS name>')
+parser.add_option("-m", "--msin", dest="msin", default="",
+ help="Measurementset name")
+
+options, args = parser.parse_args()
+
+MS = options.msin
+
+if MS == "":
+ print "No Measurement set given, exiting"
+ sys.exit(1)
+else:
+ if not os.path.exists(MS):
+ print "MS does not exist, exiting"
+ sys.exit(2)
+
+try:
+ os.environ['AIPSPATH']
+except:
+ print "AIPSPATH not set, exiting"
+ sys.exit(3)
+
+Bandpass_parset = Parset()
+Bandpass_parset['ms'] = MS
+Bandpass_parset['fixed'] = 5
+Bandpass_parset['window'] = 1
+Bandpass_parset.writeToFile("CS1_BandpassCorrector.parset")
+fd = open("CS1_BandpassCorrector.debug", 'w')
+fd.write("Global 20\n")
+fd.write("Everything 20\n")
+fd.close()
+
+Flagger_parset = Parset()
+Flagger_parset['ms'] = MS
+Flagger_parset['existing'] = False
+Flagger_parset['threshold']= 1.1
+Flagger_parset.writeToFile("CS1_frequencyFlagger.parset")
+fd = open("CS1_FrequencyFlagger.debug", 'w')
+fd.write("Global 20\n")
+fd.write("Everything 20\n")
+fd.close()
+
+Squasher_parset = Parset()
+(head, tail) = os.path.split(MS)
+Squasher_parset['inms'] = MS
+Squasher_parset['outms'] = head + "/s" + tail
+Squasher_parset['start'] = 0
+Squasher_parset['step'] = 64
+Squasher_parset['nchan'] = 256
+Squasher_parset['threshold'] = 0
+Squasher_parset['useflags'] = True
+Squasher_parset['allcolumns'] = True
+Squasher_parset.writeToFile("CS1_DataSquasher.parset")
+fd = open("CS1_DataSquasher.debug", 'w')
+fd.write("Global 20\n")
+fd.write("Everything 20\n")
+fd.close()
+
+os.system("/app/LOFAR/stable/CS1_BandpassCorrector")
+os.system("/app/LOFAR/stable/CS1_FrequencyFlagger")
+os.system("/app/LOFAR/stable/CS1_DataSquasher")
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/FlaggerStatistics.cc b/Appl/CEP/CS1/CS1_pp_lib/src/FlaggerStatistics.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a32d7e89fb8bfa37c5b1355d1d069e5f6eebf6cf
--- /dev/null
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/FlaggerStatistics.cc
@@ -0,0 +1,565 @@
+/***************************************************************************
+ * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * renting@astron.nl *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#include <tables/Tables.h>
+#include <tables/Tables/TableIter.h>
+#include "ComplexMedianFlagger.h"
+#include <casa/Quanta/MVEpoch.h>
+
+using namespace WSRT; //should be nfra ?
+using namespace casa;
+
+enum CorrelationTypes {None=0,I=1,Q=2,U=3,V=4,RR=5,RL=6,LR=7,LL=8,XX=9,XY=10,YX=11,YY=12}; //found somewhere in AIPS++, don't remember where
+
+//===============>>> itoa <<<===============
+//ANSI C++ doesn't seem to have a decent function for this, or I'm not aware of it. Need to rename it to IntToStr(), to avoid confusion
+std::string itoa(int value, int base=10)
+{ // found on http://www.jb.man.ac.uk/~slowe/cpp/itoa.html
+ //maybe copyright Robert Jan Schaper, Ray-Yuan Sheu, Rodrigo de Salvo Braz, Wes Garland and John Maloney
+ enum { kMaxDigits = 35 };
+ std::string buf;
+ buf.reserve( kMaxDigits ); // Pre-allocate enough space.
+ // check that the base if valid
+ if (base < 2 || base > 16) return buf;
+ int quotient = value;
+ // Translating number to string with base:
+ do {
+ buf += "0123456789abcdef"[ std::abs( quotient % base ) ];
+ quotient /= base;
+ } while ( quotient );
+ // Append the negative sign for base 10
+ if ( value < 0 && base == 10) buf += '-';
+ std::reverse( buf.begin(), buf.end() );
+ return buf;
+}
+
+//===============>>> ComplexMedianFlagger::ComplexMedianFlagger <<<===============
+/* initialize some meta data and get the datastorage the right size. */
+ComplexMedianFlagger::ComplexMedianFlagger(MS_File* InputMSfile,
+ int InputWindowSize,
+ bool UseOnlyXpolarizations,
+ double InputMinThreshold,
+ double InputMaxThreshold)
+{
+ MSfile = InputMSfile;
+ WindowSize = InputWindowSize;
+ MinThreshold = InputMinThreshold;
+ MaxThreshold = InputMaxThreshold;
+ NumAntennae = (*MSfile).itsNumAntennae;
+ NumPairs = (*MSfile).itsNumPairs;
+ NumBands = (*MSfile).itsNumBands;
+ NumChannels = (*MSfile).itsNumChannels;
+ NumPolarizations = (*MSfile).itsNumPolarizations;
+ NoiseLevel = (*MSfile).itsNoiseLevel;
+ NumTimeslots = (*MSfile).itsNumTimeslots;
+ AntennaNames = (*MSfile).itsAntennaNames;
+
+ PairsIndex.resize(NumPairs);
+ Statistics = Cube<int>(NumBands, NumAntennae, NumAntennae, 0);
+ PolarizationsToCheck.resize(NumPolarizations);
+ DeterminePolarizationsToCheck(UseOnlyXpolarizations);
+
+ TimeslotData.resize(NumPairs*NumBands);
+ for (int i = 0; i < NumPairs*NumBands; i++)
+ { TimeslotData[i].resize(NumPolarizations, NumChannels, WindowSize);
+ }
+
+ int index = 0;
+ for (int i = 0; i < NumAntennae; i++)
+ { for(int j = i; j < NumAntennae; j++)
+ { PairsIndex[index] = pairii(i, j);
+ BaselineIndex[pairii(i, j)] = index++;
+ }
+ }
+ ComputeBaselineLengths();
+}
+
+//===============>>> ComplexMedianFlagger::~ComplexMedianFlagger <<<===============
+
+ComplexMedianFlagger::~ComplexMedianFlagger()
+{
+}
+
+//===============>>> ComplexMedianFlagger::DetermineCorrelationsToCheck <<<===============
+/* create a list of polarizations we want to check, maybe we only want to to XY, YX */
+void ComplexMedianFlagger::DeterminePolarizationsToCheck(bool UseOnlyXpolarizations)
+{
+ if (UseOnlyXpolarizations)
+ {
+ bool noCorrError = true;
+ for (int i = 0; i < NumPolarizations; i++)
+ {
+ switch((*MSfile).itsPolarizations[i])
+ {
+ case None:
+ case I:
+ case RR:
+ case LL:
+ case XX:
+ case YY:
+ if (UseOnlyXpolarizations)
+ PolarizationsToCheck[i] = false;
+ break;
+ case Q:
+ case U:
+ case V:
+ case RL:
+ case LR:
+ case XY:
+ case YX:
+ noCorrError = false;
+ if (UseOnlyXpolarizations)
+ PolarizationsToCheck[i] = true;
+ break;
+ }
+ }
+ if (noCorrError)
+ {
+ cout << "There are no crosspolarizations to flag!";
+ exit(1);
+ }
+ }
+ else
+ {
+ for (int i = 0; i < NumPolarizations; i++)
+ { PolarizationsToCheck[i] = true;
+ }
+ }
+}
+
+//===============>>> ComplexMedianFlagger::ComputeBaselineLengths <<<===============
+/* compute baseline lengths, and determine the longest one.*/
+void ComplexMedianFlagger::ComputeBaselineLengths()
+{
+ MaxBaselineLength = 0.0;
+ BaselineLengths.resize(NumPairs);
+ //Antenna positions
+ MSAntenna antenna = (*MSfile).antenna();
+ ROArrayColumn<Double> position(antenna, "POSITION");
+ for (int i = 0; i < NumAntennae; i++ )
+ {
+ for (int j = i; j < NumAntennae; j++)
+ {
+ Vector<Double> p(position(i) - position(j));
+ double temp = sqrt(p(0)*p(0) + p(1)*p(1) + p(2)*p(2));
+ BaselineLengths[BaselineIndex[pairii(i, j)]] = temp;
+ if (temp > MaxBaselineLength && temp < 3000000) //radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243
+ { MaxBaselineLength = temp; // non-existent antenna's can have position (0,0,0)
+ }
+ }
+ }
+}
+
+//===============>>> ComplexMedianFlagger::FlagDataOrBaselines <<<===============
+/*This function outputs the gathered statistics.*/
+void ComplexMedianFlagger::ProcessStatistics()
+{
+ vector<int> bands(NumBands);
+ vector<int> antennae(NumAntennae);
+ unsigned int namelength = 6;
+ for(int i = 0; i < NumAntennae; i++)
+ {
+ if (namelength < AntennaNames[i].size())
+ { namelength = AntennaNames[i].size();
+ }
+ }
+ for (int i = 0; i < NumBands; i++)
+ {
+ cout << "Band: " << i+1 << endl;
+ cout << string(namelength+1,' ');
+ for(int j = 0; j < NumAntennae; j++)
+ {
+ string out = AntennaNames[j];
+ out.resize(namelength+1,' ');
+ cout << out;
+ }
+ cout << endl;
+ for(int j = 0; j < NumAntennae; j++)
+ {
+ string out = AntennaNames[j];
+ out.resize(namelength+1,' ');
+ cout << out;
+ for(int k = 0; k < NumAntennae; k++)
+ {
+ if (k < j) //We print a complete array, but we have inspected only those where k >= j
+ {
+ int val = 100 * Statistics(i,k,j) / (NumChannels * NumTimeslots * NumPolarizations);
+ bands[i] += val;
+ antennae[j] += val;
+ antennae[k] += val;
+ string out = itoa(val) + "%";
+ out.resize(namelength+1,' ');
+ cout << out;
+ }
+ else
+ {
+ int val = 100 * Statistics(i,j,k) / (NumChannels * NumTimeslots * NumPolarizations);
+ bands[i] += val;
+ antennae[j] += val;
+ antennae[k] += val;
+ string out = itoa(val) + "%";
+ out.resize(namelength+1,' ');
+ cout << out;
+ }
+ /*{ cout << rms << "Faulty baseline detected: Antenna " << j+1
+ << ":" << k+1 << " SpectralWindow: "<< i+1 << endl;
+ }
+ }*/
+ }
+ cout << endl;
+ }
+ }
+ cout << "Bands (flagged %): ";
+ for (int i = 0; i < NumBands; i++)
+ {
+ string out = string("IVC") + itoa(i);
+ out.resize(namelength+1,' ');
+ cout << out;
+ }
+ cout << endl << " ";
+ for (int i = 0; i < NumBands; i++)
+ {
+ string out = itoa(bands[i] / (NumAntennae*NumAntennae)) + "%";
+ out.resize(namelength+1,' ');
+ cout << out;
+ }
+ cout << endl << "Antennae (flagged %): " ;
+ for(int j = 0; j < NumAntennae; j++)
+ {
+ string out = AntennaNames[j];
+ out.resize(namelength+1,' ');
+ cout << out;
+ }
+ cout << endl << " ";
+ for (int i = 0; i < NumAntennae; i++)
+ {
+ string out = itoa(antennae[i] / (NumBands*NumAntennae*2)) + "%";
+ out.resize(namelength+1,' ');
+ cout << out;
+ }
+ cout << endl;
+}
+//===============>>> ComplexMedianFlagger::FlagTimeslot <<<===============
+/* This function inspects each visibility in a cetain baseline-band
+and flags on complexe distance, then determines to flag the entire baseline-band
+based on the RMS of the points it didn't flag.*/
+bool ComplexMedianFlagger::FlagBaselineBand(Matrix<Bool>* Flags,
+ Cube<Complex>* Timeslots,
+ vector<double>* rms,
+ int* flagCounter,
+ double FlagThreshold,
+ int Position,
+ bool flagRMS)
+{
+ Vector<Float> Reals(WindowSize);
+ Vector<Float> Imags(WindowSize);
+ vector<double> RMS(NumPolarizations);
+ vector<int> RMSCounter(NumPolarizations);
+ bool FlagCompleteRow = true;
+ int flagcount = 0;
+ for (int i = NumChannels-1; i >= 0; i--)
+ {
+ bool FlagAllPolarizations = false;
+ for (int j = NumPolarizations-1; j >= 0; j--)
+ { //we need to loop twice, once to determine FlagAllCorrelations
+ if (!FlagAllPolarizations && PolarizationsToCheck[j])
+ {
+ for (int k = 0; k < WindowSize; k++)
+ { //This might be faster in some other way ?
+ Reals[k] = (*Timeslots)(j, i, k).real();
+ Imags[k] = (*Timeslots)(j, i, k).imag();
+ }
+ float real = medianInPlace(Reals, false);
+ float imag = medianInPlace(Imags, false);
+ FlagAllPolarizations |= FlagThreshold < abs((*Timeslots)(j, i, Position)
+ - Complex(real, imag));
+ }
+ }
+ for (int j = NumPolarizations-1; j >= 0; j--)
+ { //the second loop we set the flags or calculate RMS
+ if (FlagAllPolarizations)
+ { (*Flags)(j, i) = true;
+ flagcount++;
+ }
+ else
+ {
+ FlagCompleteRow = false;
+ (*Flags)(j, i) = false; //could check existing flag ?
+ double value = abs((*Timeslots)(j, i, Position));
+ RMS[j] += value*value;
+ (RMSCounter[j])++;
+ }
+ }
+ }
+ //these need to be separated out into a different function for clarity
+ if (flagcount > 0.9 * NumChannels * NumPolarizations) //more as 90% bad
+ {
+ FlagCompleteRow = true;
+ (*Flags) = true;
+ flagcount = NumChannels * NumPolarizations;
+ }
+ if (!FlagCompleteRow && flagRMS) //RMSCounter > 0
+ {
+ for (int j = NumPolarizations-1; j >= 0; j--)
+ {
+ (*rms)[j] = sqrt( RMS[j] / RMSCounter[j] );
+// if (PolarizationsToCheck[j] && ((*rms)[j] > 1.5 || (*rms)[j] <= 0.5 * NoiseLevel)) //we think this RMS value is non-physical
+ if (PolarizationsToCheck[j] && ((*rms)[j] > 25 || (*rms)[j] <= 0.5 * NoiseLevel)) //we think this RMS value is non-physical
+ {
+ FlagCompleteRow = true;
+ (*Flags) = true;
+ flagcount = NumChannels * NumPolarizations;
+ }
+ }
+ }
+ (*flagCounter) += flagcount;
+ return FlagCompleteRow;
+}
+//===============>>> ComplexMedianFlagger::FlagBaseline <<<===============
+/* This function iterates over baseline and band and uses FlagBaselineBand() to determine
+ for each one if it needs to be flagged. It treats the autocorrelations separately,
+ to detect entire malfunctioning telescopes. Finally it writes the flags.
+*/
+void ComplexMedianFlagger::FlagTimeslot(TableIterator* flag_iter,
+ bool flagDatapoints,
+ bool flagRMS,
+ bool ExistingFlags,
+ Int Position)
+{
+ Cube<bool> Flags(NumPolarizations, NumChannels, NumPairs*NumBands, false);
+ vector<bool> FlagCompleteRow(NumPairs*NumBands);
+ vector<int> stats(NumPairs*NumBands);
+ for (int i = 0; i < NumBands; i++)
+ {
+ vector<double> rms(NumPolarizations);
+ vector< vector<double> > AutoCorrRMS(NumPolarizations, NumAntennae);
+ for (int j = 0; j < NumPolarizations; j++) // this needs to be written more elegantly
+ {
+ if ((*MSfile).itsPolarizations[j] == YX)
+ { PolarizationsToCheck[j] = false; // We're not checking YX because WSRT writes 0.0 into them
+ }
+ }
+ for(int j = 0; j < NumAntennae; j++)
+ { //check auto correlations
+ int index = i*NumPairs + BaselineIndex[pairii(j, j)];
+ Matrix<Bool> flags = Flags.xyPlane(index);
+ FlagCompleteRow[index] = FlagBaselineBand(&flags,
+ &(TimeslotData[index]),
+ &rms, &(stats[index]),
+ 3000 * NoiseLevel,
+ Position, flagRMS);
+ for(int k = 0; k < NumPolarizations; k++)
+ { AutoCorrRMS[k][j] = rms[k];
+ }
+ }
+ for(int k = 0; k < NumPolarizations; k++) //check on median
+ {
+ if (PolarizationsToCheck[k])
+ {
+ double median = medianInPlace(Vector<double>(AutoCorrRMS[k]), false);
+ for(int j = 0; j < NumAntennae; j++)
+ {
+ int index = i*NumPairs + BaselineIndex[pairii(j, j)];
+ if (AutoCorrRMS[k][j] < median/10 || AutoCorrRMS[k][j] > median*10)
+ { FlagCompleteRow[index] = true;
+ }
+ }
+ }
+ }
+ for (int j = 0; j < NumPolarizations; j++) // this needs to be written more elegantly
+ {
+ if ((*MSfile).itsPolarizations[j] == YX)
+ { PolarizationsToCheck[j] = true; // We're not checking YX
+ }
+ }
+ // this code coud be separated into different functions
+ vector< vector<double> > CrossCorrRMS(NumPolarizations, NumPairs - NumAntennae);
+ for(int j = 0; j < NumAntennae; j++)
+ {
+ for(int k = j+1; k < NumAntennae; k++)
+ {
+ int index = i*NumPairs + BaselineIndex[pairii(j, k)];
+ Matrix<Bool> flags = Flags.xyPlane(index);
+ if ((BaselineLengths[BaselineIndex[pairii(j, k)]] < 3000000))//radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243 m
+ { //we skip the non-existent telescopes
+ FlagCompleteRow[index] = FlagBaselineBand(&flags,
+ &(TimeslotData[index]),
+ &rms, &(stats[index]),
+ (MinThreshold + (MaxThreshold - MinThreshold)
+ * BaselineLengths[BaselineIndex[pairii(j, k)]]
+ / MaxBaselineLength
+ ) * NoiseLevel,
+ Position, flagRMS);
+ for(int l = 0; l < NumPolarizations; l++)
+ { CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] = rms[l]; //a very complex way of indexing, could be simplified?
+ }
+ if ( ( FlagCompleteRow[i*NumPairs + BaselineIndex[pairii(j, j)]]
+ ||FlagCompleteRow[i*NumPairs + BaselineIndex[pairii(k, k)]])
+ &&flagRMS && !FlagCompleteRow[index])
+ { //We flag it if the autocorrelation is flagged
+ FlagCompleteRow[index] = true;
+ }
+ }
+ }
+ }
+ for(int l = 0; l < NumPolarizations; l++) // check on median
+ {
+ if (PolarizationsToCheck[l])
+ {
+ double median = medianInPlace(Vector<double>(CrossCorrRMS[l]), false);
+ for(int j = 0; j < NumAntennae; j++)
+ {
+ for(int k = j+1; k < NumAntennae; k++)
+ {
+ int index = i*NumPairs + BaselineIndex[pairii(j, k)];
+ if (!FlagCompleteRow[index])
+ { if( CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] < median/10 //a very complex way of indexing, could be simplified?
+ || CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] > median*10)
+ { FlagCompleteRow[index] = true;
+ }
+ }
+ }
+ }
+ }
+ }
+ //this code could be separated in to different functions
+ for(int j = 0; j < NumAntennae; j++) //write the data
+ {
+ for(int k = j; k < NumAntennae; k++)
+ {
+ int index = i*NumPairs + BaselineIndex[pairii(j, k)];
+ Matrix<Bool> flags = Flags.xyPlane(index);
+ if (FlagCompleteRow[index])
+ {
+ Matrix<Bool> flags = Flags.xyPlane(index);
+ flags = true;
+ Statistics(i,j,k) += NumChannels * NumPolarizations;
+ }
+ else
+ { Statistics(i,j,k) += stats[index];
+ }
+ (*MSfile).WriteDataPointFlags(flag_iter, &flags, FlagCompleteRow[index], ExistingFlags);
+ (*flag_iter)++;
+ }
+ }
+ }
+}
+
+//===============>>> ComplexMedianFlagger::UpdateTimeslotData <<<===============
+/* This function reads the visibility data for one timeslot and checks if
+ a mosaicing mode is active*/
+/* The datastructure Timeslotdata is rather complex because it flattens
+ half-filled antenna x antenna x bands matrix into a vector of NumPairs X bands length. */
+bool ComplexMedianFlagger::UpdateTimeslotData(vector<int>* OldFields,
+ vector<int>* OldBands,
+ int* TimeCounter,
+ Table* TimeslotTable,
+ double* Time)
+{
+ int rowcount = (*TimeslotTable).nrow();
+ ROTableVector<Int> antenna1((*TimeslotTable), "ANTENNA1");
+ ROTableVector<Int> antenna2((*TimeslotTable), "ANTENNA2");
+ ROTableVector<Int> bandnr ((*TimeslotTable), "DATA_DESC_ID");
+ ROTableVector<Int> fieldid ((*TimeslotTable), "FIELD_ID");
+ ROTableVector<Double> time ((*TimeslotTable), "TIME_CENTROID");//for testing purposes
+ ROArrayColumn<Complex> data ((*TimeslotTable), "DATA");
+ Cube<Complex> Data(NumPolarizations, NumChannels, rowcount);
+
+ data.getColumn(Data); //We're not checking Data.nrow() Data.ncolumn(), assuming all data is the same size.
+ (*Time) = time(0);//for testing purposes, might be useful in the future
+
+ bool NewField_or_Frequency = false;
+
+ for (int i = 0; i < rowcount; i++)
+ {
+ int bi = BaselineIndex[pairii(antenna1(i), antenna2(i))];
+ int field = fieldid(i);
+ int band = bandnr(i);
+ int index = (band % NumBands) * NumPairs + bi;
+ if ((*TimeCounter) > WindowSize - 1)
+ {
+ if (field != (*OldFields)[bi]) //pointing mosaicing
+ { NewField_or_Frequency = true;
+ }
+ if (band != (*OldBands)[index]) //frequency mosaicing
+ { NewField_or_Frequency = true;
+ }
+ }
+ (*OldFields)[bi] = field;
+ (*OldBands)[index] = band;
+
+ TimeslotData[index].xyPlane((*TimeCounter) % WindowSize) = Data.xyPlane(i); //TimeslotData is only WindowSize size!
+ }
+ (*TimeCounter)++; //we want to reset if we detect a gap in DATA_DESC_ID or FIELD. Maybe TIME too???
+ return NewField_or_Frequency; //assuming everybody is changing field or frequency at the same time!
+}
+
+//===============>>> ComplexMedianFlagger::FlagDataOrBaselines <<<===============
+/* This function iterates over the data per timeslot and uses Flagtimeslot()
+ to actually flag datapoints (if flagDatapoints), and entire baselines (if flagRMS)*/
+void ComplexMedianFlagger::FlagDataOrBaselines(bool flagDatapoints,
+ bool flagRMS,
+ bool ExistingFlags)
+{
+ TableIterator timeslot_iter = (*MSfile).TimeslotIterator();
+ TableIterator flag_iter = (*MSfile).TimeAntennaIterator();
+ int TimeCounter = 0;
+ double Time = 0.0;//for testing purposes
+ vector<int> OldFields(NumPairs); //to check on multipointing and mosaicing
+ vector<int> OldBands(NumPairs*NumBands); //to check on multifrequency and freq mosaicing
+
+ int step = NumTimeslots / 10 + 1; //not exact but it'll do
+ int row = 0;
+
+ while (!timeslot_iter.pastEnd())
+ {
+ Table TimeslotTable = timeslot_iter.table();
+ bool NewFieldorFreq = UpdateTimeslotData(&OldFields,
+ &OldBands,
+ &TimeCounter,
+ &TimeslotTable,
+ &Time);
+ //cout << "Processing: " << MVTime(Time/(24*3600)).string(MVTime::YMD) << endl; //for testing purposes
+
+ if (TimeCounter == WindowSize - 1)
+ { //We have filled WindowSize timeslots and need to flag the first WindowSize/2 timeslots
+ for (int position = 0; position < WindowSize/2; position++)
+ { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags, position);
+ }
+ }
+ if (TimeCounter > WindowSize - 1) //nothing special just falg a timeslot
+ { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags,
+ (TimeCounter + WindowSize/2) % WindowSize);
+ }
+ timeslot_iter++;
+ if (row++ % step == 0) // to tell the user how much % we have processed,
+ { cout << 10*(row/step) << "%" << endl; //not very accurate for low numbers of timeslots, but it'll do for now
+ }
+ if (timeslot_iter.pastEnd() || NewFieldorFreq)
+ { //We still need to flag the last WindowSize/2 timeslots
+ for (int position = WindowSize/2 + 1; position < WindowSize; position++)
+ { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags,
+ (TimeCounter + position) % WindowSize);
+ }
+ TimeCounter = 0; //reset because we have changed to a new Field or frequency
+ }
+ }
+ ProcessStatistics(); //is there a baseline that should be flagged?
+}
+//===============>>> ComplexMedianFlagger <<<===============
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/FlaggerStatistics.h b/Appl/CEP/CS1/CS1_pp_lib/src/FlaggerStatistics.h
new file mode 100644
index 0000000000000000000000000000000000000000..0619007dc1b81bd000b9c0dc48536c5bc891836e
--- /dev/null
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/FlaggerStatistics.h
@@ -0,0 +1,103 @@
+/***************************************************************************
+ * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * renting@astron.nl *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#ifndef __FLAGGER_COMPLEXMEDIANFLAGGER_H__
+#define __FLAGGER_COMPLEXMEDIANFLAGGER_H__
+
+#include <casa/Arrays.h>
+#include <utility>
+#include <vector>
+#include <list>
+#include <map>
+#include "MS_File.h"
+
+namespace WSRT
+{
+ using casa::Cube;
+ using casa::Matrix;
+ using casa::Complex;
+ using casa::Int;
+ using casa::Double;
+ using casa::Bool;
+ using casa::Table;
+ using std::list;
+ using std::vector;
+
+ typedef pair<int, int> pairii;
+
+ class ComplexMedianFlagger
+ {
+ public:
+ ComplexMedianFlagger(MS_File* MSfile,
+ int InputWindowSize,
+ bool UseOnlyXpolarizations,
+ double InputMinThreshold,
+ double InputMaxThreshold);
+ ~ComplexMedianFlagger();
+
+ void FlagDataOrBaselines(bool flagDatapoints,
+ bool flagRMS,
+ bool ExistingFlags);
+
+ protected:
+ int NumAntennae;
+ int NumPairs;
+ int NumBands;
+ int NumChannels;
+ int NumPolarizations;
+ int WindowSize;
+ int NumTimeslots;
+ double MinThreshold;
+ double MaxThreshold;
+ double MaxBaselineLength;
+ double NoiseLevel;
+ vector<double> BaselineLengths;
+ vector<pairii> PairsIndex;
+ map<pairii, int> BaselineIndex;
+ vector< Cube<Complex> > TimeslotData;
+ vector< bool > PolarizationsToCheck;
+ vector<casa::String> AntennaNames;
+ Cube< int > Statistics;
+ MS_File* MSfile;
+
+ void DeterminePolarizationsToCheck(Bool UseOnlyXpolarizations);
+ void ComputeBaselineLengths();
+ void ProcessStatistics();
+ void FlagTimeslot(casa::TableIterator* flag_iter,
+ bool flagDatapoints,
+ bool flagRMS,
+ bool ExistingFlags,
+ int Position);
+ bool FlagBaselineBand(Matrix<Bool>* Flags,
+ Cube<Complex>* Timeslots,
+ vector<double>* rms,
+ int* flagCounter,
+ double FlagThreshold,
+ int Position,
+ bool flagRMS);
+ bool UpdateTimeslotData(vector<int>* OldFields,
+ vector<int>* OldBands,
+ int* TimeCounter,
+ Table* TimeslotTable,
+ double* Time);
+ private:
+ }; // ComplexMedianFlagger
+}; // namespace WSRT
+
+#endif // __FLAGGER_COMPLEXMEDIANFLAGGER_H__
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/MsInfo.cc b/Appl/CEP/CS1/CS1_pp_lib/src/MsInfo.cc
new file mode 100644
index 0000000000000000000000000000000000000000..926f0f63c874a35244e993183ff98ac39e60beca
--- /dev/null
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/MsInfo.cc
@@ -0,0 +1,156 @@
+/***************************************************************************
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
+ * renting@astron.nl *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#include <casa/BasicMath/Math.h>
+#include <casa/Arrays.h>
+
+#include <iostream>
+
+#include "MsInfo.h"
+
+using namespace LOFAR;
+using namespace casa;
+
+//===============>>> Ms_Info::Ms_Info <<<===============
+
+MsInfo::MsInfo(const std::string& msname)
+{
+ MsName = msname;
+ MaxBaselineLength = 0.0;
+ NumSamples = 0;
+ NumAntennae = 0;
+ NumFields = 0;
+ NumBands = 0;
+ NumChannels = 0;
+ NumPolarizations = 0;
+ NumPairs = 0;
+ NumTimeslots = 0;
+ NoiseLevel = 0.0;
+}
+
+//===============>>> Ms_Info::~Ms_Info <<<===============
+
+MsInfo::~MsInfo()
+{
+}
+
+//===============>>> Ms_Info::update <<<===============
+
+void MsInfo::Update(void)
+{
+ MeasurementSet MS(MsName);
+
+ //Number of samples
+ NumSamples = MS.nrow();
+ //Number of Fields
+ MSField fields = MS.field();
+ NumFields = fields.nrow();
+
+ //Number of Antennae
+ MSAntenna antennae = MS.antenna();
+ NumAntennae = antennae.nrow();
+
+ //Antenna Names
+ ROScalarColumn<String> ANT_NAME_col(antennae, "NAME");
+ Vector<String> ant_names = ANT_NAME_col.getColumn();
+ ant_names.tovector(AntennaNames);
+
+ //Number of channels in the Band
+ MSSpectralWindow spectral_window = MS.spectralWindow();
+ ROScalarColumn<Int> NUM_CHAN_col(spectral_window, "NUM_CHAN");
+ NumChannels = NUM_CHAN_col(0);
+
+ //Number of polarizations
+ MSPolarization polarization = MS.polarization();
+ ROScalarColumn<Int> NUM_CORR_col(polarization, "NUM_CORR");
+ NumPolarizations = NUM_CORR_col(0);
+ ROArrayColumn<Int> CORR_TYPE_col(polarization, "CORR_TYPE");
+ Polarizations.resize(NumPolarizations);
+ CORR_TYPE_col.get(0, Polarizations);
+
+ //calculate theoretical noise level
+ ROScalarColumn<Double> EXPOSURE_col(MS, "EXPOSURE");
+ Double exposure = EXPOSURE_col(0);
+
+ ROScalarColumn<Double> TOTAL_BANDWIDTH_col(spectral_window, "TOTAL_BANDWIDTH");
+ Double bandwidth = TOTAL_BANDWIDTH_col(0) / NumChannels;
+
+ NoiseLevel = 1.0 / sqrt(bandwidth * exposure);
+
+ //calculate number of timeslots
+ ROScalarColumn<Double> INTERVAL_col(MS, "INTERVAL");
+ Double Interval = INTERVAL_col(0);
+
+ //Number of timeslots
+ ROScalarColumn<Double> TIME_CENTROID_col(MS, "TIME_CENTROID");
+ Double firstdate = TIME_CENTROID_col(0);
+ Double lastdate = TIME_CENTROID_col(NumSamples-1);
+
+ NumTimeslots = (int)((lastdate-firstdate)/Interval) + 1;
+
+ //calculate number of baselines.
+ NumPairs = (NumAntennae) * (NumAntennae + 1) / 2; //Triangular numbers formula
+
+ //calculate number of Bands
+ NumBands = NumSamples / (NumPairs * NumTimeslots);
+
+ ComputeBaselineLengths(MS);
+}
+
+//===============>>> MS_Info::PrintInfo <<<===============
+
+void MsInfo::PrintInfo(void)
+{
+ std::cout << "NumSamples: " << NumSamples << std::endl;
+ std::cout << "NumFields: " << NumFields << std::endl;
+ std::cout << "NumAntennae: " << NumAntennae << std::endl;
+ std::cout << "NumChannels: " << NumChannels << std::endl;
+ std::cout << "NumPolarizations: " << NumPolarizations << std::endl;
+ std::cout << "Noiselevel: " << NoiseLevel << std::endl;
+ std::cout << "Interval: " << Interval << std::endl;
+ std::cout << "Numtimeslots: " << NumTimeslots << std::endl;
+ std::cout << "NumPairs: " << NumPairs << std::endl;
+ std::cout << "NumBands: " << NumBands << std::endl;
+}
+
+//===============>>> MS_Info::ComputeBaselineLengths <<<===============
+/* compute baseline lengths, and determine the longest one.*/
+
+void MsInfo::ComputeBaselineLengths(casa::MeasurementSet& MS)
+{
+ BaselineLengths.resize(NumPairs);
+ //Antenna positions
+ MSAntenna antenna = MS.antenna();
+ ROArrayColumn<Double> position(antenna, "POSITION");
+ for (int i = 0; i < NumAntennae; i++ )
+ {
+ for (int j = i; j < NumAntennae; j++)
+ {
+ Vector<Double> p(position(i) - position(j));
+ double temp = sqrt(p(0)*p(0) + p(1)*p(1) + p(2)*p(2));
+ BaselineLengths[BaselineIndex[baseline_t(i, j)]] = temp;
+ if (temp > MaxBaselineLength && temp < 3000000) //radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243
+ { MaxBaselineLength = temp; // non-existent antenna's can have position (0,0,0)
+ }
+ }
+ }
+}
+
+//===============>>> MsInfo <<<===============
+
diff --git a/Appl/CEP/CS1/CS1_pp_lib/src/MsInfo.h b/Appl/CEP/CS1/CS1_pp_lib/src/MsInfo.h
new file mode 100644
index 0000000000000000000000000000000000000000..e3ca327e2d5c5fd95b9b2f03820c4b749764ba2f
--- /dev/null
+++ b/Appl/CEP/CS1/CS1_pp_lib/src/MsInfo.h
@@ -0,0 +1,67 @@
+/***************************************************************************
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
+ * renting@astron.nl *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#ifndef __CS1_PP_MS_INFO_H__
+#define __CS1_PP_MS_INFO_H__
+
+#include <ms/MeasurementSets.h>
+#include <tables/Tables.h>
+#include <utility>
+#include <vector>
+#include <string>
+#include <map>
+
+namespace LOFAR
+{
+ typedef std::pair<int, int> baseline_t;
+
+ class MsInfo
+ {
+ public:
+ MsInfo(const std::string& msname);
+ ~MsInfo();
+
+ int NumSamples;
+ int NumAntennae;
+ int NumFields;
+ int NumBands;
+ int NumChannels;
+ int NumPolarizations;
+ int NumPairs;
+ int NumTimeslots;
+ double NoiseLevel;
+ std::vector<casa::String> AntennaNames;
+ casa::Vector<casa::Int> Polarizations;
+ double MaxBaselineLength;
+ std::vector<double> BaselineLengths;
+ std::vector<baseline_t> PairsIndex;
+ std::map<baseline_t, int> BaselineIndex;
+
+ void Update(void);
+ void PrintInfo(void);
+
+ protected:
+ private:
+ std::string MsName;
+ double Interval;
+ void ComputeBaselineLengths(casa::MeasurementSet& MS);
+ }; // MsInfo
+}; // namespace LOFAR
+
+#endif // __CS1_PP_MS_INFO_H__