diff --git a/.gitattributes b/.gitattributes
index 4c7633cb43f4816537bac0fa90a5fc6972308339..1bb14e74b0092b366e0c04e3dd8d4a598882a4cd 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -709,6 +709,7 @@ CEP/DP3/DPPP/test/CS1_IDPPP.parset -text
CEP/DP3/DPPP/test/tNDPPP.in_MS.tgz -text svneol=unset#application/x-compressed-tar
CEP/DP3/DPPP/test/tmwflagger.in_cd -text
CEP/DP3/DPPP/test/tmwflagger.in_vd -text
+CEP/Imager/ImageLofar/src/addImagingInfo -text
CEP/Imager/LofarFT/TODO.txt -text
CEP/Imager/LofarFT/include/LofarFT/LofarATerm.h -text
CEP/Imager/LofarFT/include/LofarFT/LofarCFDefs.h -text
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/gui/quality/histogrampage.h b/CEP/DP3/AOFlagger/include/AOFlagger/gui/quality/histogrampage.h
index 9d2cbd6035eeea957bc1d0588bae2029cd0e1df5..365b08d9b3574e790790f7c1b7cff1dbbd02dada 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/gui/quality/histogrampage.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/gui/quality/histogrampage.h
@@ -62,11 +62,11 @@ class HistogramPage : public Gtk::HBox {
void plotPolarization(const HistogramCollection &histogramCollection, unsigned polarization);
void plotPolarization(const class LogHistogram &totalHistogram, const class LogHistogram &rfiHistogram);
void plotFit(const class LogHistogram &histogram, const std::string &title);
- void plotSlope(const class LogHistogram &histogram, const std::string &title);
+ void plotSlope(const class LogHistogram &histogram, const std::string &title, bool useLowerLimit2);
void onPlotPropertiesClicked();
void onDataExportClicked();
void readFromFile();
- void updateSlopeFrame();
+ void updateSlopeFrame(const LogHistogram &histogram);
void addSlopeText(std::stringstream &str, const LogHistogram &histogram, bool updateRange);
void updateDataWindow();
@@ -112,7 +112,7 @@ class HistogramPage : public Gtk::HBox {
Gtk::Frame _slopeFrame;
Gtk::VBox _slopeBox;
Gtk::TextView _slopeTextView;
- Gtk::CheckButton _drawSlopeButton;
+ Gtk::CheckButton _drawSlopeButton, _drawSlope2Button;
Gtk::CheckButton _slopeAutoRangeButton;
Gtk::Entry _slopeStartEntry, _slopeEndEntry, _slopeRFIRatio;
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/quality/histogramtablesformatter.h b/CEP/DP3/AOFlagger/include/AOFlagger/quality/histogramtablesformatter.h
index 3fa65cff496931f014f7acf8257d3023c25e9d67..96624903d8cd1470a0be2d4c2e2a3b7bc6ce6dda 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/quality/histogramtablesformatter.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/quality/histogramtablesformatter.h
@@ -141,8 +141,10 @@ class HistogramTablesFormatter {
{
Close();
openMainTable(true);
- _measurementSet->rwKeywordSet().removeField(TableName(table));
- casa::Table::deleteTable(TableFilename(table));
+ if(_measurementSet->keywordSet().isDefined(TableName(table)))
+ _measurementSet->rwKeywordSet().removeField(TableName(table));
+ if(_measurementSet->isReadable(TableFilename(table)))
+ casa::Table::deleteTable(TableFilename(table));
}
}
@@ -165,6 +167,11 @@ class HistogramTablesFormatter {
{
return TableExists(HistogramCountTable) && TableExists(HistogramTypeTable);
}
+ void RemoveAll()
+ {
+ RemoveTable(HistogramCountTable);
+ RemoveTable(HistogramTypeTable);
+ }
private:
HistogramTablesFormatter(const HistogramTablesFormatter &) { } // don't allow copies
void operator=(const HistogramTablesFormatter &) { } // don't allow assignment
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h b/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h
index 4f335323d7aeca79ed57eaafc953174ea186ef23..3daf1efc9fc8ce4b1624a802b962b57f80ff26f5 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h
@@ -263,6 +263,34 @@ class LogHistogram : public Serializable
return (double) (sumOffset/(long double) n);
}
+ double NormalizedSlopeStdDev(double startAmplitude, double endAmplitude, double slope, double offset) const
+ {
+ long double sqErrorSum = 0.0, xsqErrorSum = 0.0, xSum = 0.0;
+ unsigned long n = 0;
+ // determine the 'average' x
+ for(const_iterator i=begin();i!=end();++i)
+ {
+ if(i.value() >= startAmplitude && i.value() < endAmplitude)
+ {
+ xSum += log10(i.value());
+ ++n;
+ }
+ }
+ const long double avgX = xSum / n;
+ for(const_iterator i=begin();i!=end();++i)
+ {
+ if(i.value() >= startAmplitude && i.value() < endAmplitude)
+ {
+ long double y = log10(i.normalizedCount());
+ long double x = log10(i.value());
+ long double ySlope = x*slope + offset;
+ sqErrorSum += (y-ySlope)*(y-ySlope);
+ xsqErrorSum += (x - avgX)*(x - avgX);
+ }
+ }
+ return (double) sqrtl(sqErrorSum/(xsqErrorSum * (long double) (n-2)));
+ }
+
double PowerLawUpperLimit(double constrainingAmplitude, double exponent, double factor) const
{
const double count = NormalizedCountAbove(constrainingAmplitude);
@@ -278,6 +306,14 @@ class LogHistogram : public Serializable
return pow(term, 1.0/(exponent+1.0));
}
+ double PowerLawLowerLimit2(double constrainingAmplitude, double exponent, double factor, double rfiRatio) const
+ {
+ const double countPart = NormalizedCountAbove(constrainingAmplitude);
+ const double countTotal = NormalizedTotalCount() * rfiRatio;
+ const double term = (countPart - countTotal) * (exponent+1.0)/factor + pow(constrainingAmplitude, exponent+1.0);
+ return pow(term/-exponent, 1.0/(exponent+1.0));
+ }
+
void GetRFIRegion(double &start, double &end) const
{
double sigmaEstimate = AmplitudeWithMaxNormalizedCount();
diff --git a/CEP/DP3/AOFlagger/src/aoquality.cpp b/CEP/DP3/AOFlagger/src/aoquality.cpp
index 3712941f5a6a37a305faaf8575aa06a0ffbc4c8b..0bfe8f39939c54501e42f7356384979438604490 100644
--- a/CEP/DP3/AOFlagger/src/aoquality.cpp
+++ b/CEP/DP3/AOFlagger/src/aoquality.cpp
@@ -427,6 +427,23 @@ void actionRemove(const std::string &filename)
formatter.RemoveAllQualityTables();
}
+void printRFISlopeForHistogram(const std::map<HistogramCollection::AntennaPair, LogHistogram*> &histogramMap, char polarizationSymbol, const AntennaInfo *antennae)
+{
+ for(std::map<HistogramCollection::AntennaPair, LogHistogram*>::const_iterator i=histogramMap.begin(); i!=histogramMap.end();++i)
+ {
+ const unsigned a1 = i->first.first, a2 = i->first.second;
+ Baseline baseline(antennae[a1], antennae[a2]);
+ double length = baseline.Distance();
+ const LogHistogram &histogram = *i->second;
+ double start, end;
+ histogram.GetRFIRegion(start, end);
+ double slope = histogram.NormalizedSlope(start, end);
+ double offset = histogram.NormalizedSlopeOffset(start, end, slope);
+ double stddev = histogram.NormalizedSlopeStdDev(start, end, slope, offset);
+ std::cout << polarizationSymbol << '\t' << a1 << '\t' << a2 << '\t' << length << '\t' << slope << '\t' << stddev << '\n';
+ }
+}
+
void actionHistogram(const std::string &filename, const std::string &query)
{
HistogramTablesFormatter histogramFormatter(filename);
@@ -454,19 +471,21 @@ void actionHistogram(const std::string &filename, const std::string &query)
for(size_t a=0;a<antennaCount;++a)
antennae[a] = set.GetAntennaInfo(a);
+ HistogramCollection *summedCollection = collection.CreateSummedPolarizationCollection();
+ const std::map<HistogramCollection::AntennaPair, LogHistogram*> &histogramMap = summedCollection->GetRFIHistogram(0);
+ printRFISlopeForHistogram(histogramMap, '*', antennae);
+ delete summedCollection;
for(unsigned p=0;p<polarizationCount;++p)
{
const std::map<HistogramCollection::AntennaPair, LogHistogram*> &histogramMap = collection.GetRFIHistogram(p);
- for(std::map<HistogramCollection::AntennaPair, LogHistogram*>::const_iterator i=histogramMap.begin(); i!=histogramMap.end();++i)
- {
- const unsigned a1 = i->first.first, a2 = i->first.second;
- Baseline baseline(antennae[a1], antennae[a2]);
- double length = baseline.Distance();
- const LogHistogram &histogram = *i->second;
- double slope = histogram.NormalizedSlopeInRFIRegion();
- std::cout << p << '\t' << a1 << '\t' << a2 << '\t' << length << '\t' << slope << '\n';
- }
+ printRFISlopeForHistogram(histogramMap, '0' + p, antennae);
}
+ } else if(query == "remove")
+ {
+ histogramFormatter.RemoveAll();
+ } else
+ {
+ std::cerr << "Unknown histogram command: " << query << "\n";
}
}
diff --git a/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp b/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp
index 3053bae58f951ca74fc852ca080f7b8c65018342..cd60769e8260caa43b43f37a09b8c712d315eb58 100644
--- a/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp
+++ b/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp
@@ -52,6 +52,7 @@ HistogramPage::HistogramPage() :
_dataExportButton("Data"),
_slopeFrame("Slope"),
_drawSlopeButton("Draw"),
+ _drawSlope2Button("Draw2"),
_slopeAutoRangeButton("Auto range"),
_plotPropertiesWindow(0),
_histograms(0),
@@ -131,6 +132,8 @@ HistogramPage::HistogramPage() :
_slopeBox.pack_start(_slopeTextView, Gtk::PACK_SHRINK);
_drawSlopeButton.signal_clicked().connect(sigc::mem_fun(*this, &HistogramPage::updatePlot));
_slopeBox.pack_start(_drawSlopeButton, Gtk::PACK_SHRINK);
+ _drawSlope2Button.signal_clicked().connect(sigc::mem_fun(*this, &HistogramPage::updatePlot));
+ _slopeBox.pack_start(_drawSlope2Button, Gtk::PACK_SHRINK);
_slopeBox.pack_start(_slopeAutoRangeButton, Gtk::PACK_SHRINK);
_slopeAutoRangeButton.set_active(true);
@@ -224,7 +227,6 @@ void HistogramPage::updatePlot()
plotPolarization(*_summedPolarizationHistograms, 0);
_plotWidget.Update();
- updateSlopeFrame();
updateDataWindow();
}
}
@@ -248,6 +250,15 @@ void HistogramPage::plotPolarization(const LogHistogram &totalHistogram, const L
{
plotFit(totalHistogram, "Fit to total");
}
+ if(_drawSlopeButton.get_active())
+ {
+ plotSlope(totalHistogram, "Fitted slope", false);
+ }
+ if(_drawSlope2Button.get_active())
+ {
+ plotSlope(totalHistogram, "Fitted slope", true);
+ }
+ updateSlopeFrame(totalHistogram);
}
if(_rfiHistogramButton.get_active())
@@ -259,10 +270,15 @@ void HistogramPage::plotPolarization(const LogHistogram &totalHistogram, const L
{
plotFit(rfiHistogram, "Fit to RFI");
}
- }
- if(_drawSlopeButton.get_active())
- {
- plotSlope(rfiHistogram, "Fitted slope");
+ updateSlopeFrame(rfiHistogram);
+ if(_drawSlopeButton.get_active())
+ {
+ plotSlope(rfiHistogram, "Fitted slope", false);
+ }
+ if(_drawSlope2Button.get_active())
+ {
+ plotSlope(rfiHistogram, "Fitted slope", true);
+ }
}
if(_notRFIHistogramButton.get_active())
@@ -375,7 +391,7 @@ void HistogramPage::addRayleighDifferenceToPlot(const LogHistogram &histogram, d
}
}
-void HistogramPage::plotSlope(const LogHistogram &histogram, const std::string &title)
+void HistogramPage::plotSlope(const LogHistogram &histogram, const std::string &title, bool useLowerLimit2)
{
double start, end;
if(_slopeAutoRangeButton.get_active())
@@ -385,13 +401,30 @@ void HistogramPage::plotSlope(const LogHistogram &histogram, const std::string &
start = atof(_slopeStartEntry.get_text().c_str());
end = atof(_slopeEndEntry.get_text().c_str());
}
- double slope = histogram.NormalizedSlope(start, end);
- double offset = histogram.NormalizedSlopeOffset(start, end, slope);
+ double
+ xMin = log10(histogram.MinPositiveAmplitude()),
+ rfiRatio = atof(_slopeRFIRatio.get_text().c_str()),
+ slope = histogram.NormalizedSlope(start, end),
+ offset = histogram.NormalizedSlopeOffset(start, end, slope),
+ upperLimit = log10(histogram.PowerLawUpperLimit(start, slope, pow10(offset))),
+ lowerLimit = useLowerLimit2 ?
+ log10(histogram.PowerLawLowerLimit2(start, slope, pow10(offset), rfiRatio)) :
+ log10(histogram.PowerLawLowerLimit(start, slope, pow10(offset), rfiRatio));
+ double xStart, xEnd;
+ if(std::isfinite(lowerLimit))
+ xStart = lowerLimit;
+ else
+ xStart = log10(start) - 1.0;
+ if(std::isfinite(upperLimit))
+ xEnd = upperLimit;
+ else
+ xEnd = log10(histogram.MaxAmplitude());
+ double
+ yStart = xStart*slope + offset,
+ yEnd = xEnd*slope + offset;
_plot.StartLine(title, "Amplitude in arbitrary units (log)", "Frequency (log)");
- double xStart = log10(start / 10.0);
- double xEnd = log10(histogram.MaxAmplitude());
- double yStart = xStart*slope + offset;
- double yEnd = xEnd*slope + offset;
+ if(useLowerLimit2 && std::isfinite(xMin))
+ _plot.PushDataPoint(xMin, yStart);
_plot.PushDataPoint(xStart, yStart);
_plot.PushDataPoint(xEnd, yEnd);
}
@@ -415,21 +448,13 @@ void HistogramPage::onDataExportClicked()
updateDataWindow();
}
-void HistogramPage::updateSlopeFrame()
+void HistogramPage::updateSlopeFrame(const LogHistogram &histogram)
{
std::stringstream str;
str << "Slopes:";
- LogHistogram summedHistogram;
- _summedPolarizationHistograms->GetRFIHistogramForCrossCorrelations(0, summedHistogram);
- addSlopeText(str, summedHistogram, true);
+ addSlopeText(str, histogram, true);
- for(size_t p=0;p<_histograms->PolarizationCount();++p)
- {
- LogHistogram histogram;
- _histograms->GetRFIHistogramForCrossCorrelations(p, histogram);
- addSlopeText(str, histogram, false);
- }
_slopeTextView.get_buffer()->set_text(str.str());
}
@@ -453,11 +478,17 @@ void HistogramPage::addSlopeText(std::stringstream &str, const LogHistogram &his
}
double rfiRatio = atof(_slopeRFIRatio.get_text().c_str());
- double slope = histogram.NormalizedSlope(minRange, maxRange);
- double offset = histogram.NormalizedSlopeOffset(minRange, maxRange, slope);
- double upperLimit = histogram.PowerLawUpperLimit(minRange, slope, pow10(offset));
- double lowerLimit = histogram.PowerLawLowerLimit(minRange, slope, pow10(offset), rfiRatio);
- str << '\n' << slope << '[' << log10(lowerLimit) << ';' << log10(upperLimit) << ']';
+ const double
+ slope = histogram.NormalizedSlope(minRange, maxRange),
+ offset = histogram.NormalizedSlopeOffset(minRange, maxRange, slope),
+ error = histogram.NormalizedSlopeStdDev(minRange, maxRange, slope, offset),
+ upperLimit = histogram.PowerLawUpperLimit(minRange, slope, pow10(offset)),
+ lowerLimit = histogram.PowerLawLowerLimit(minRange, slope, pow10(offset), rfiRatio),
+ lowerError = fabs(lowerLimit - histogram.PowerLawLowerLimit(minRange, slope - error, pow10(offset), rfiRatio)),
+ lowerLimit2 = histogram.PowerLawLowerLimit2(minRange, slope, pow10(offset), rfiRatio);
+ str << '\n' << slope << "±" << error << "\n["
+ << log10(lowerLimit) << "±" << lowerError << ';' << log10(upperLimit) << ']' << '\n'
+ << log10(lowerLimit2);
}
void HistogramPage::updateDataWindow()
diff --git a/CEP/Imager/ImageLofar/CMakeLists.txt b/CEP/Imager/ImageLofar/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5c8f01efa2dd9479e09d633d91687aa8ea4fe0cd
--- /dev/null
+++ b/CEP/Imager/ImageLofar/CMakeLists.txt
@@ -0,0 +1,9 @@
+# $Id$
+
+lofar_package(ImageLofar 0.1 DEPENDS Common)
+
+lofar_find_package(Casacore REQUIRED COMPONENTS images)
+
+add_subdirectory(include/ImageLofar)
+add_subdirectory(src)
+add_subdirectory(test)
diff --git a/CEP/Imager/ImageLofar/include/ImageLofar/CMakeLists.txt b/CEP/Imager/ImageLofar/include/ImageLofar/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4b42128ae8df3c4f115b7984101f140e0e4e82b2
--- /dev/null
+++ b/CEP/Imager/ImageLofar/include/ImageLofar/CMakeLists.txt
@@ -0,0 +1,14 @@
+# $Id$
+
+# List of header files that will be installed.
+set(inst_HEADERS
+)
+
+# Create symbolic link to include directory.
+execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_BINARY_DIR}/include/${PACKAGE_NAME})
+
+# Install header files.
+install(FILES ${inst_HEADERS} DESTINATION include/${PACKAGE_NAME})
+
diff --git a/CEP/Imager/ImageLofar/src/CMakeLists.txt b/CEP/Imager/ImageLofar/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fedd4639e3a871d82413a4ca043bfadfdeaf3c9c
--- /dev/null
+++ b/CEP/Imager/ImageLofar/src/CMakeLists.txt
@@ -0,0 +1,11 @@
+# $Id$
+
+include(LofarPackageVersion)
+
+lofar_add_library(lofarimage
+ Package__Version.cc
+)
+
+install(PROGRAMS
+ addImagingInfo
+ DESTINATION bin)
diff --git a/CEP/Imager/ImageLofar/src/addImagingInfo b/CEP/Imager/ImageLofar/src/addImagingInfo
new file mode 100755
index 0000000000000000000000000000000000000000..70e110c9b703cd7fbadc7e18f8ea5ded28b3d842
--- /dev/null
+++ b/CEP/Imager/ImageLofar/src/addImagingInfo
@@ -0,0 +1,219 @@
+#!/usr/bin/env python
+
+import pyrap.tables as pt
+import lofar.parmdb as pdb
+
+""" Add a subtable of an MS to the image """
+def addSubTable (image, msName, subName, removeColumns=[]):
+ # Make a selection of all rows/columns of the MS subtable
+ sel = pt.taql ("select * from '" + msName + "/" + subName + "'")
+ # Remove the required columns.
+ if len(removeColumns) > 0:
+ sel.removecols (removeColumns)
+ # Strip LOFAR_ from column names
+ for col in sel.colnames():
+ if len(col) > 6 and col[:6] == "LOFAR_":
+ sel.renamecol (col, col[6:])
+ # Copy the subtable to the image and add it as a subtable.
+ # Always prefix subtable name with LOFAR_.
+ subNameOut = subName;
+ if len(subNameOut) < 6 or subNameOut[:6] != "LOFAR_":
+ subNameOut = "LOFAR_" + subNameOut
+ subtab = sel.copy (image.name() + "/" + subNameOut, deep=True)
+ image.putkeyword (subNameOut, subtab)
+ print "Added subtable", subNameOut, "containing", subtab.nrows(), "rows"
+ subtab.close()
+ sel.close()
+
+""" Create the empty LOFAR_QUALITY subtable """
+def addQualityTable (image):
+ # Create the table using TaQL.
+ tab = pt.taql ("create table '" + image.name() + "/LOFAR_QUALITY' " +
+ "QUALITY_MEASURE string, VALUE string, FLAG_ROW bool")
+ image.putkeyword ("LOFAR_QUALITY", tab)
+ tab.close()
+ print "Added subtable LOFAR_QUALITY containing 0 rows"
+
+""" Create the LOFAR_ORIGIN subtable and fill from all MSs """
+def addOriginTable (image, msNames):
+ # Concatenate the OBSERVATION subtables of all MSs.
+ obsNames = [name + "/OBSERVATION" for name in msNames]
+ obstab = pt.table(obsNames, ack=False)
+ # Select and rename the required columns.
+ # Some columns are not in the LOFAR_OBSERVATION table.
+ # Create them by selecting a similarly typed column and fill them later.
+ selstr = "LOFAR_OBSERVATION_ID as OBSERVATION_ID"
+ selstr += ",LOFAR_SUB_ARRAY_POINTING as SUB_ARRAY_POINTING"
+ selstr += ",LOFAR_SUB_ARRAY_POINTING as SUBBAND"
+ selstr += ",LOFAR_SUB_ARRAY_POINTING as NUM_CHAN"
+ selstr += ",LOFAR_OBSERVATION_FREQUENCY_MIN as CHANNEL_WIDTH"
+ selstr += ",LOFAR_OBSERVATION_FREQUENCY_MIN as EXPOSURE"
+ selstr += ",LOFAR_OBSERVATION_FREQUENCY_MIN as FREQUENCY_MIN"
+ selstr += ",LOFAR_OBSERVATION_FREQUENCY_MAX as FREQUENCY_MAX"
+ selstr += ",LOFAR_OBSERVATION_FREQUENCY_CENTER as FREQUENCY_CENTER"
+ selstr += ",LOFAR_OBSERVATION_START as START"
+ selstr += ",LOFAR_OBSERVATION_END as END"
+ selstr += ",FLAG_ROW"
+ sel = obstab.select(selstr)
+ # Copy the subtable to the image and add it as a subtable.
+ subtab = sel.copy (image.name() + "/" + "LOFAR_ORIGIN", deep=True)
+ subtab = pt.table (image.name() + "/" + "LOFAR_ORIGIN", readonly=False,
+ ack=False)
+ obstab.close()
+ image.putkeyword ("LOFAR_ORIGIN", subtab)
+ # Set the correct units of columns to update.
+ subtab.putcolkeyword ("CHANNEL_WIDTH", "QuantumUnits", ["Hz"])
+ subtab.putcolkeyword ("EXPOSURE", "QuantumUnits", ["s"])
+ # Update the columns not in OBSERVATION table.
+ # Get EXPOSURE from first row in main tables.
+ # Get NUM_CHAN from SPECTRAL_WINDOW subtables.
+ # Calculate CHANNEL_WIDTH (convert from MHz to Hz).
+ # Get SUBBAND from MS name.
+ for i in range(len(msNames)):
+ t = pt.table(msNames[i], ack=False)
+ subtab.putcell ("EXPOSURE", i, t.getcell("EXPOSURE", 0))
+ t1 = pt.table(t.getkeyword("SPECTRAL_WINDOW"), ack=False)
+ numchan = t1.getcell("NUM_CHAN", 0)
+ subtab.putcell ("NUM_CHAN", i, numchan)
+ w = subtab.getcell("FREQUENCY_MAX",i) - subtab.getcell("FREQUENCY_MIN",i)
+ subtab.putcell ("CHANNEL_WIDTH", i, w * 1e6 / numchan)
+ t1.close()
+ t.close()
+ subband = 0
+ inx = msNames[i].find ("SB")
+ if inx>= 0:
+ try:
+ subband = int(msNames[i][inx+2:inx+5])
+ except:
+ pass
+ subtab.putcell ("SUBBAND", i, subband)
+ # Ready
+ subtab.close()
+ sel.close()
+ print "Added subtable LOFAR_ORIGIN containing", len(msNames), "rows"
+
+""" Create the LOFAR_SOURCE subtable and fill from the SourceDB """
+def addSourceTable (image, sourcedbName, times):
+ # Create the table using TaQL.
+ tab = pt.taql ("create table '" + image.name() + "/LOFAR_SOURCE' " +
+ "SOURCE_ID int, \TIME double, INTERVAL double, " +
+ "NUM_LINES int, NAME string, " +
+ "DIRECTION double shape=[2], " +
+ "PROPER_MOTION double shape=[2], " +
+ "FLUX double shape=[4], " +
+ "SPINX double, " +
+ "SHAPE double shape=[3]")
+ tab.putcolkeyword ("TIME", "QuantumUnits", ["s"])
+ tab.putcolkeyword ("INTERVAL", "QuantumUnits", ["s"])
+ tab.putcolkeyword ("DIRECTION", "QuantumUnits", ["rad"])
+ tab.putcolkeyword ("PROPER_MOTION", "QuantumUnits", ["rad/s"])
+ tab.putcolkeyword ("FLUX", "QuantumUnits", ["Jy"])
+ tab.putcolkeyword ("SHAPE", "QuantumUnits", ["rad", "rad", "rad"])
+ tab.putcolkeyword ("TIME", "MEASINFO", {"Ref":"UTC", "type":"epoch"})
+ tab.putcolkeyword ("DIRECTION", "MEASINFO", {"Ref":"J2000", "type":"direction"})
+ image.putkeyword ("LOFAR_SOURCE", tab)
+ # Get all parameters from the source parmdb.
+ midtime = (times[0] + times[1]) / 2
+ inttime = times[1] - times[0]
+ sourcedb = pdb.parmdb(sourcedbName)
+ # Get all source names by getting the Ra parms from DEFAULTVALUES
+ names = [name[3:] for name in sourcedb.getDefNames ("Ra:*")]
+ values = sourcedb.getDefValues()
+ sourcedb = 0 # close
+ row = 0
+ tab.addrows (len(names))
+ # Add the info of all sources.
+ # The field names below are as used in SourceDB.
+ fldnames = ["Ra", "Dec", "I", "Q", "U", "V", "SpectralIndex:0",
+ "Orientation", "MajorAxis", "MinorAxis"]
+ vals = [0. for fld in fldnames]
+ for name in names:
+ for i in range(len(fldnames)):
+ key = fldnames[i] + ":" + name
+ if values.has_key (key):
+ vals[i] = values[key][0][0]
+ else:
+ vals[i] = 0.
+ tab.putcell ("SOURCE_ID", row, row)
+ tab.putcell ("TIME", row, midtime);
+ tab.putcell ("INTERVAL", row, inttime);
+ tab.putcell ("NUM_LINES", row, 0);
+ tab.putcell ("NAME", row, name);
+ tab.putcell ("DIRECTION", row, vals[:2]);
+ tab.putcell ("PROPER_MOTION", row, (0.,0.));
+ tab.putcell ("FLUX", row, vals[2:6]);
+ tab.putcell ("SPINX", row, vals[6]);
+ tab.putcell ("SHAPE", row, vals[7:10]);
+ print name, vals
+ row += 1
+ # Ready.
+ tab.close()
+ print "Added subtable LOFAR_SOURCE containing", row, "rows"
+
+""" Update times and frequencies in the LOFAR_OBSERVATION subtable """
+def updateObsTable (image):
+ obstab = pt.table (image.name() + "/LOFAR_OBSERVATION", readonly=False,
+ ack=False)
+ oritab = pt.table (image.name() + "/LOFAR_ORIGIN", ack=False)
+ minfreq = pt.taql ("calc min([select FREQUENCY_MIN from '" +
+ oritab.name() + "'])")
+ maxfreq = pt.taql ("calc max([select FREQUENCY_MAX from '" +
+ oritab.name() + "'])")
+ mintime = pt.taql ("calc min([select START from '" +
+ oritab.name() + "'])")
+ maxtime = pt.taql ("calc max([select END from '" +
+ oritab.name() + "'])")
+ obstab.putcell ("OBSERVATION_FREQUENCY_MIN", 0, minfreq[0]);
+ obstab.putcell ("OBSERVATION_FREQUENCY_MAX", 0, maxfreq[0]);
+ obstab.putcell ("OBSERVATION_FREQUENCY_CENTER", 0, (maxfreq[0]-minfreq[0])/2);
+ obstab.putcell ("OBSERVATION_START", 0, mintime[0]);
+ obstab.putcell ("OBSERVATION_END", 0, maxtime[0]);
+ obstab.putcell ("TIME_RANGE", 0, (mintime[0], maxtime[0]));
+ obstab.putcell ("FILETYPE", 0, "img")
+ pt.taql ("update '" + obstab.name() + "' set FILEDATE = mjd(date()), " +
+ "RELEASE_DATE = mjd(date()+365)")
+ obstab.close()
+ oritab.close()
+ print "Updated subtable LOFAR_OBSERVATION"
+ return (mintime[0], maxtime[0])
+
+""" Add all imaging info """
+def addImagingInfo (imageName, msNames, sourcedbName):
+ image = pt.table (imageName, readonly=False, ack=False)
+ # Add all subtables while removing obsolete columns
+ addSubTable (image, msNames[0], "POINTING")
+ addSubTable (image, msNames[0], "FIELD")
+ addSubTable (image, msNames[0], "ANTENNA")
+ addSubTable (image, msNames[0], "LOFAR_STATION")
+ addSubTable (image, msNames[0], "HISTORY")
+ addSubTable (image, msNames[0], "OBSERVATION",
+ ["LOG", "SCHEDULE_TYPE", "SCHEDULE"])
+ # Create the (empty) LOFAR_QUALITY subtable.
+ addQualityTable (image)
+ # Create the LOFAR_ORIGIN subtable from all MSs.
+ addOriginTable (image, msNames)
+ # Update times/frequencies in the LOFAR_OBSERVATION table.
+ times = updateObsTable (image)
+ # Add the LOFAR_SOURCE table.
+ addSourceTable (image, sourcedbName, times)
+ # Flush and close the image.
+ image.close()
+
+
+if __name__ == "__main__":
+ import sys
+ if len(sys.argv) < 4:
+ print "Insufficient arguments; run as:"
+ print " addImagingInfo image sourcedb ms1 ms2 ..."
+ print " image name of the image"
+ print " sourcedb name of SourceDB containing the sources found"
+ print " ms1 ms2 ... names of MSs image is made of"
+ print " the names can be individual arguments or a"
+ print " comma-separated list of names (or a mix)"
+ print " E.g.,"
+ print " addImagingInfo myimg mysdb ms1,ms2 ms3 ms4,ms5"
+ sys.exit(1)
+ msNames = []
+ for arg in sys.argv[3:]:
+ msNames.extend (arg.split(","))
+ addImagingInfo (sys.argv[1], msNames, sys.argv[2]);
diff --git a/CEP/Imager/ImageLofar/test/CMakeLists.txt b/CEP/Imager/ImageLofar/test/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fd22edb2aa39156a7e567f10fff2c0fdddb0330a
--- /dev/null
+++ b/CEP/Imager/ImageLofar/test/CMakeLists.txt
@@ -0,0 +1,5 @@
+# $Id$
+
+include(LofarCTest)
+
+#lofar_add_test(tLofarConvolutionFunction tLofarConvolutionFunction.cc)
diff --git a/LCS/ApplCommon/test/tObservation.stdout b/LCS/ApplCommon/test/tObservation.stdout
index dfd78a1103e353cc59f960eeb7cb71cf1f62d782..f187a53712491e8d2f908e0d47e3f5a3001b8687 100644
--- a/LCS/ApplCommon/test/tObservation.stdout
+++ b/LCS/ApplCommon/test/tObservation.stdout
@@ -1,5 +1,4 @@
>>>
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
Observation :
ObsID : 5025
@@ -36,7 +35,6 @@ Beam[0].TAB[1]: 5.990000, 0.490000, J2000, 20.000000, incoherent
beamlet2beams : [0,0,0,0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999]
nrAnaBeams : 0
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
Observation :
ObsID : 5025
@@ -74,7 +72,6 @@ beamlet2beams : [0,0,0,0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1
nrAnaBeams : 0
<<<
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
Observation :
ObsID : 5025
@@ -122,22 +119,15 @@ TABringsize : 0
beamlet2beams : [0,0,0,0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1,1,-1,1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999]
nrAnaBeams : 0
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
>>>
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
INFO LCS.ApplCommon - Clock of observation 5025 and 5026 conflict
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
INFO LCS.ApplCommon - Conflicting use of receivers between observation 5025 and 5027
DEBUG LCS.ApplCommon - receiverConflict: 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111100000000
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
INFO LCS.ApplCommon - Conflict in beamlets between observation 5025 and 5028
DEBUG LCS.ApplCommon - First conflicting beamlet: 2
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
INFO LCS.ApplCommon - Conflict in nrSlotsInFrame: 48<->54 for resp. observation 5025 and 5029
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
<<<
No conflict found in file 5 which is oke.
-DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
getRCUbitset(96,48,'') = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111
getRCUbitset(96,96,'') = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111
getRCUbitset(96,48,LBA_XXX) = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111
@@ -159,7 +149,6 @@ LBA_OUTER(true) = LBA
LBA_X(false) = LBA
LBA_X(true) = LBA
>>>
-DEBUG LCS.ApplCommon - Dataslots.CS001HBA0.DataslotList
OLD SYNTAX
Observation :
@@ -228,7 +217,6 @@ TABringsize : 0
beamlet2beams : [0,0,0,0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,999,999,999,999,999,999,999,999,999,999,999,999,999,1,1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,999,999,999,999,999,999,999,999,999,999,999,999,999]
nrAnaBeams : 0
-DEBUG LCS.ApplCommon - Dataslots.CS001HBA0.DataslotList
NEW SYNTAX
Observation :