diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h b/CEP/DP3/AOFlagger/include/AOFlagger/quality/loghistogram.h index 4f335323d7aeca79ed57eaafc953174ea186ef23..6b10279e8e3e9bc8095606a37d07ce1a8591f868 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); diff --git a/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp b/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp index 3053bae58f951ca74fc852ca080f7b8c65018342..0f7115f5f751339e62b32adfca3fa77577abba9e 100644 --- a/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp +++ b/CEP/DP3/AOFlagger/src/gui/quality/histogrampage.cpp @@ -453,11 +453,13 @@ 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); + str << '\n' << slope << "±" << error << '[' << log10(lowerLimit) << ';' << log10(upperLimit) << ']'; } void HistogramPage::updateDataWindow()