Skip to content
Snippets Groups Projects
Commit e39609e8 authored by Jan David Mol's avatar Jan David Mol
Browse files

Task #3049: Merged trunk with branch of task #3064

parents 982b1b57 bd132766
No related branches found
No related tags found
No related merge requests found
Showing
with 393 additions and 53 deletions
...@@ -709,6 +709,7 @@ CEP/DP3/DPPP/test/CS1_IDPPP.parset -text ...@@ -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/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_cd -text
CEP/DP3/DPPP/test/tmwflagger.in_vd -text CEP/DP3/DPPP/test/tmwflagger.in_vd -text
CEP/Imager/ImageLofar/src/addImagingInfo -text
CEP/Imager/LofarFT/TODO.txt -text CEP/Imager/LofarFT/TODO.txt -text
CEP/Imager/LofarFT/include/LofarFT/LofarATerm.h -text CEP/Imager/LofarFT/include/LofarFT/LofarATerm.h -text
CEP/Imager/LofarFT/include/LofarFT/LofarCFDefs.h -text CEP/Imager/LofarFT/include/LofarFT/LofarCFDefs.h -text
......
...@@ -62,11 +62,11 @@ class HistogramPage : public Gtk::HBox { ...@@ -62,11 +62,11 @@ class HistogramPage : public Gtk::HBox {
void plotPolarization(const HistogramCollection &histogramCollection, unsigned polarization); void plotPolarization(const HistogramCollection &histogramCollection, unsigned polarization);
void plotPolarization(const class LogHistogram &totalHistogram, const class LogHistogram &rfiHistogram); void plotPolarization(const class LogHistogram &totalHistogram, const class LogHistogram &rfiHistogram);
void plotFit(const class LogHistogram &histogram, const std::string &title); 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 onPlotPropertiesClicked();
void onDataExportClicked(); void onDataExportClicked();
void readFromFile(); void readFromFile();
void updateSlopeFrame(); void updateSlopeFrame(const LogHistogram &histogram);
void addSlopeText(std::stringstream &str, const LogHistogram &histogram, bool updateRange); void addSlopeText(std::stringstream &str, const LogHistogram &histogram, bool updateRange);
void updateDataWindow(); void updateDataWindow();
...@@ -112,7 +112,7 @@ class HistogramPage : public Gtk::HBox { ...@@ -112,7 +112,7 @@ class HistogramPage : public Gtk::HBox {
Gtk::Frame _slopeFrame; Gtk::Frame _slopeFrame;
Gtk::VBox _slopeBox; Gtk::VBox _slopeBox;
Gtk::TextView _slopeTextView; Gtk::TextView _slopeTextView;
Gtk::CheckButton _drawSlopeButton; Gtk::CheckButton _drawSlopeButton, _drawSlope2Button;
Gtk::CheckButton _slopeAutoRangeButton; Gtk::CheckButton _slopeAutoRangeButton;
Gtk::Entry _slopeStartEntry, _slopeEndEntry, _slopeRFIRatio; Gtk::Entry _slopeStartEntry, _slopeEndEntry, _slopeRFIRatio;
......
...@@ -141,7 +141,9 @@ class HistogramTablesFormatter { ...@@ -141,7 +141,9 @@ class HistogramTablesFormatter {
{ {
Close(); Close();
openMainTable(true); openMainTable(true);
if(_measurementSet->keywordSet().isDefined(TableName(table)))
_measurementSet->rwKeywordSet().removeField(TableName(table)); _measurementSet->rwKeywordSet().removeField(TableName(table));
if(_measurementSet->isReadable(TableFilename(table)))
casa::Table::deleteTable(TableFilename(table)); casa::Table::deleteTable(TableFilename(table));
} }
} }
...@@ -165,6 +167,11 @@ class HistogramTablesFormatter { ...@@ -165,6 +167,11 @@ class HistogramTablesFormatter {
{ {
return TableExists(HistogramCountTable) && TableExists(HistogramTypeTable); return TableExists(HistogramCountTable) && TableExists(HistogramTypeTable);
} }
void RemoveAll()
{
RemoveTable(HistogramCountTable);
RemoveTable(HistogramTypeTable);
}
private: private:
HistogramTablesFormatter(const HistogramTablesFormatter &) { } // don't allow copies HistogramTablesFormatter(const HistogramTablesFormatter &) { } // don't allow copies
void operator=(const HistogramTablesFormatter &) { } // don't allow assignment void operator=(const HistogramTablesFormatter &) { } // don't allow assignment
......
...@@ -263,6 +263,34 @@ class LogHistogram : public Serializable ...@@ -263,6 +263,34 @@ class LogHistogram : public Serializable
return (double) (sumOffset/(long double) n); 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 double PowerLawUpperLimit(double constrainingAmplitude, double exponent, double factor) const
{ {
const double count = NormalizedCountAbove(constrainingAmplitude); const double count = NormalizedCountAbove(constrainingAmplitude);
...@@ -278,6 +306,14 @@ class LogHistogram : public Serializable ...@@ -278,6 +306,14 @@ class LogHistogram : public Serializable
return pow(term, 1.0/(exponent+1.0)); 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 void GetRFIRegion(double &start, double &end) const
{ {
double sigmaEstimate = AmplitudeWithMaxNormalizedCount(); double sigmaEstimate = AmplitudeWithMaxNormalizedCount();
......
...@@ -427,6 +427,23 @@ void actionRemove(const std::string &filename) ...@@ -427,6 +427,23 @@ void actionRemove(const std::string &filename)
formatter.RemoveAllQualityTables(); 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) void actionHistogram(const std::string &filename, const std::string &query)
{ {
HistogramTablesFormatter histogramFormatter(filename); HistogramTablesFormatter histogramFormatter(filename);
...@@ -454,19 +471,21 @@ void actionHistogram(const std::string &filename, const std::string &query) ...@@ -454,19 +471,21 @@ void actionHistogram(const std::string &filename, const std::string &query)
for(size_t a=0;a<antennaCount;++a) for(size_t a=0;a<antennaCount;++a)
antennae[a] = set.GetAntennaInfo(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) for(unsigned p=0;p<polarizationCount;++p)
{ {
const std::map<HistogramCollection::AntennaPair, LogHistogram*> &histogramMap = collection.GetRFIHistogram(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) printRFISlopeForHistogram(histogramMap, '0' + p, antennae);
{
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';
}
} }
} else if(query == "remove")
{
histogramFormatter.RemoveAll();
} else
{
std::cerr << "Unknown histogram command: " << query << "\n";
} }
} }
......
...@@ -52,6 +52,7 @@ HistogramPage::HistogramPage() : ...@@ -52,6 +52,7 @@ HistogramPage::HistogramPage() :
_dataExportButton("Data"), _dataExportButton("Data"),
_slopeFrame("Slope"), _slopeFrame("Slope"),
_drawSlopeButton("Draw"), _drawSlopeButton("Draw"),
_drawSlope2Button("Draw2"),
_slopeAutoRangeButton("Auto range"), _slopeAutoRangeButton("Auto range"),
_plotPropertiesWindow(0), _plotPropertiesWindow(0),
_histograms(0), _histograms(0),
...@@ -131,6 +132,8 @@ HistogramPage::HistogramPage() : ...@@ -131,6 +132,8 @@ HistogramPage::HistogramPage() :
_slopeBox.pack_start(_slopeTextView, Gtk::PACK_SHRINK); _slopeBox.pack_start(_slopeTextView, Gtk::PACK_SHRINK);
_drawSlopeButton.signal_clicked().connect(sigc::mem_fun(*this, &HistogramPage::updatePlot)); _drawSlopeButton.signal_clicked().connect(sigc::mem_fun(*this, &HistogramPage::updatePlot));
_slopeBox.pack_start(_drawSlopeButton, Gtk::PACK_SHRINK); _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); _slopeBox.pack_start(_slopeAutoRangeButton, Gtk::PACK_SHRINK);
_slopeAutoRangeButton.set_active(true); _slopeAutoRangeButton.set_active(true);
...@@ -224,7 +227,6 @@ void HistogramPage::updatePlot() ...@@ -224,7 +227,6 @@ void HistogramPage::updatePlot()
plotPolarization(*_summedPolarizationHistograms, 0); plotPolarization(*_summedPolarizationHistograms, 0);
_plotWidget.Update(); _plotWidget.Update();
updateSlopeFrame();
updateDataWindow(); updateDataWindow();
} }
} }
...@@ -248,6 +250,15 @@ void HistogramPage::plotPolarization(const LogHistogram &totalHistogram, const L ...@@ -248,6 +250,15 @@ void HistogramPage::plotPolarization(const LogHistogram &totalHistogram, const L
{ {
plotFit(totalHistogram, "Fit to total"); 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()) if(_rfiHistogramButton.get_active())
...@@ -259,10 +270,15 @@ void HistogramPage::plotPolarization(const LogHistogram &totalHistogram, const L ...@@ -259,10 +270,15 @@ void HistogramPage::plotPolarization(const LogHistogram &totalHistogram, const L
{ {
plotFit(rfiHistogram, "Fit to RFI"); plotFit(rfiHistogram, "Fit to RFI");
} }
} updateSlopeFrame(rfiHistogram);
if(_drawSlopeButton.get_active()) if(_drawSlopeButton.get_active())
{ {
plotSlope(rfiHistogram, "Fitted slope"); plotSlope(rfiHistogram, "Fitted slope", false);
}
if(_drawSlope2Button.get_active())
{
plotSlope(rfiHistogram, "Fitted slope", true);
}
} }
if(_notRFIHistogramButton.get_active()) if(_notRFIHistogramButton.get_active())
...@@ -375,7 +391,7 @@ void HistogramPage::addRayleighDifferenceToPlot(const LogHistogram &histogram, d ...@@ -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; double start, end;
if(_slopeAutoRangeButton.get_active()) if(_slopeAutoRangeButton.get_active())
...@@ -385,13 +401,30 @@ void HistogramPage::plotSlope(const LogHistogram &histogram, const std::string & ...@@ -385,13 +401,30 @@ void HistogramPage::plotSlope(const LogHistogram &histogram, const std::string &
start = atof(_slopeStartEntry.get_text().c_str()); start = atof(_slopeStartEntry.get_text().c_str());
end = atof(_slopeEndEntry.get_text().c_str()); end = atof(_slopeEndEntry.get_text().c_str());
} }
double slope = histogram.NormalizedSlope(start, end); double
double offset = histogram.NormalizedSlopeOffset(start, end, slope); 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)"); _plot.StartLine(title, "Amplitude in arbitrary units (log)", "Frequency (log)");
double xStart = log10(start / 10.0); if(useLowerLimit2 && std::isfinite(xMin))
double xEnd = log10(histogram.MaxAmplitude()); _plot.PushDataPoint(xMin, yStart);
double yStart = xStart*slope + offset;
double yEnd = xEnd*slope + offset;
_plot.PushDataPoint(xStart, yStart); _plot.PushDataPoint(xStart, yStart);
_plot.PushDataPoint(xEnd, yEnd); _plot.PushDataPoint(xEnd, yEnd);
} }
...@@ -415,21 +448,13 @@ void HistogramPage::onDataExportClicked() ...@@ -415,21 +448,13 @@ void HistogramPage::onDataExportClicked()
updateDataWindow(); updateDataWindow();
} }
void HistogramPage::updateSlopeFrame() void HistogramPage::updateSlopeFrame(const LogHistogram &histogram)
{ {
std::stringstream str; std::stringstream str;
str << "Slopes:"; str << "Slopes:";
LogHistogram summedHistogram; addSlopeText(str, histogram, true);
_summedPolarizationHistograms->GetRFIHistogramForCrossCorrelations(0, summedHistogram);
addSlopeText(str, summedHistogram, 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()); _slopeTextView.get_buffer()->set_text(str.str());
} }
...@@ -453,11 +478,17 @@ void HistogramPage::addSlopeText(std::stringstream &str, const LogHistogram &his ...@@ -453,11 +478,17 @@ void HistogramPage::addSlopeText(std::stringstream &str, const LogHistogram &his
} }
double rfiRatio = atof(_slopeRFIRatio.get_text().c_str()); double rfiRatio = atof(_slopeRFIRatio.get_text().c_str());
double slope = histogram.NormalizedSlope(minRange, maxRange); const double
double offset = histogram.NormalizedSlopeOffset(minRange, maxRange, slope); slope = histogram.NormalizedSlope(minRange, maxRange),
double upperLimit = histogram.PowerLawUpperLimit(minRange, slope, pow10(offset)); offset = histogram.NormalizedSlopeOffset(minRange, maxRange, slope),
double lowerLimit = histogram.PowerLawLowerLimit(minRange, slope, pow10(offset), rfiRatio); error = histogram.NormalizedSlopeStdDev(minRange, maxRange, slope, offset),
str << '\n' << slope << '[' << log10(lowerLimit) << ';' << log10(upperLimit) << ']'; 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() void HistogramPage::updateDataWindow()
......
# $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)
# $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})
# $Id$
include(LofarPackageVersion)
lofar_add_library(lofarimage
Package__Version.cc
)
install(PROGRAMS
addImagingInfo
DESTINATION bin)
#!/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]);
# $Id$
include(LofarCTest)
#lofar_add_test(tLofarConvolutionFunction tLofarConvolutionFunction.cc)
>>> >>>
DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
Observation : Observation :
ObsID : 5025 ObsID : 5025
...@@ -36,7 +35,6 @@ Beam[0].TAB[1]: 5.990000, 0.490000, J2000, 20.000000, incoherent ...@@ -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] 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 nrAnaBeams : 0
DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
Observation : Observation :
ObsID : 5025 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 ...@@ -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 nrAnaBeams : 0
<<< <<<
DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
Observation : Observation :
ObsID : 5025 ObsID : 5025
...@@ -122,22 +119,15 @@ TABringsize : 0 ...@@ -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] 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 nrAnaBeams : 0
DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
>>> >>>
DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
INFO LCS.ApplCommon - Clock of observation 5025 and 5026 conflict 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 INFO LCS.ApplCommon - Conflicting use of receivers between observation 5025 and 5027
DEBUG LCS.ApplCommon - receiverConflict: 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111100000000 DEBUG LCS.ApplCommon - receiverConflict: 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111100000000
DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
INFO LCS.ApplCommon - Conflict in beamlets between observation 5025 and 5028 INFO LCS.ApplCommon - Conflict in beamlets between observation 5025 and 5028
DEBUG LCS.ApplCommon - First conflicting beamlet: 2 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 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. No conflict found in file 5 which is oke.
DEBUG LCS.ApplCommon - Dataslots.CS016HBA.DataslotList
getRCUbitset(96,48,'') = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111 getRCUbitset(96,48,'') = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111
getRCUbitset(96,96,'') = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111 getRCUbitset(96,96,'') = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111
getRCUbitset(96,48,LBA_XXX) = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111 getRCUbitset(96,48,LBA_XXX) = 000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000111111111111
...@@ -159,7 +149,6 @@ LBA_OUTER(true) = LBA ...@@ -159,7 +149,6 @@ LBA_OUTER(true) = LBA
LBA_X(false) = LBA LBA_X(false) = LBA
LBA_X(true) = LBA LBA_X(true) = LBA
>>> >>>
DEBUG LCS.ApplCommon - Dataslots.CS001HBA0.DataslotList
OLD SYNTAX OLD SYNTAX
Observation : Observation :
...@@ -228,7 +217,6 @@ TABringsize : 0 ...@@ -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] 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 nrAnaBeams : 0
DEBUG LCS.ApplCommon - Dataslots.CS001HBA0.DataslotList
NEW SYNTAX NEW SYNTAX
Observation : Observation :
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment