Commit 8ab67e8d authored by Martyn Gigg's avatar Martyn Gigg Committed by GitHub
Browse files

Merge pull request #17675 from mantidproject/17669_Histogram_rollout_algs_30-35

Histogram Rollout - Algorithms 30-35
parents 8b192fab 4b5f1ac0
......@@ -23,6 +23,7 @@
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
#include "MantidKernel/System.h"
#include "MantidHistogramData/Histogram.h"
#include <memory>
#include <utility>
#include <vector>
......@@ -61,15 +62,15 @@ public:
/// Constructor
MayersSampleCorrectionStrategy(
MayersSampleCorrectionStrategy::Parameters params,
const std::vector<double> &tof, const std::vector<double> &sigIn,
const std::vector<double> &errIn);
const Mantid::HistogramData::Histogram &inputHist);
/// Destructor - defined in cpp file to use forward declaration with
/// unique_ptr
~MayersSampleCorrectionStrategy();
/// Return the correction factors
void apply(std::vector<double> &sigOut, std::vector<double> &errOut);
/// Calculate the self-attentation factor for a single mu*r value
Mantid::HistogramData::Histogram getCorrectedHisto();
/// Calculate the self-attenuation factor for a single mu*r value
double calculateSelfAttenuation(const double muR);
/// Calculate the multiple scattering factor for a single mu*r value &
/// absorption value
......@@ -84,19 +85,15 @@ private:
double muR(const double flightPath, const double tof) const;
double muR(const double sigt) const;
double sigmaTotal(const double flightPath, const double tof) const;
double tof(const size_t i) const;
void seedRNG(const size_t seed);
/// A copy of the correction parameters
const Parameters m_pars;
/// A reference to the tof vluaes
const std::vector<double> &m_tof;
/// A reference to the input signal values
const std::vector<double> &m_sigin;
/// A reference to the input error values
const std::vector<double> &m_errin;
// True if we have binned TOF values
const bool m_histogram;
// Holds histogram to process
const HistogramData::Histogram &m_histogram;
const HistogramData::Points m_tofVals;
/// Holds the number of Y vals to process
const size_t m_histoYSize;
/// Limits for the range of mu*r values to cover
const std::pair<double, double> m_muRrange;
/// Random number generator
......
......@@ -58,19 +58,19 @@ void RRFMuon::exec() {
// Get phase
double phase = getProperty("Phase");
// Get number of histograms
size_t nHisto = inputWs->getNumberHistograms();
const size_t nHisto = inputWs->getNumberHistograms();
if (nHisto != 2) {
throw std::runtime_error("Invalid number of spectra in input workspace");
}
// Set number of data points
size_t nData = inputWs->blocksize();
const size_t nData = inputWs->blocksize();
// Compute the RRF polarization
const double twoPiFreq = 2. * M_PI * freq * factor;
MantidVec time = inputWs->readX(0); // X axis: time
MantidVec labRe = inputWs->readY(0); // Lab-frame polarization (real part)
MantidVec labIm =
inputWs->readY(1); // Lab-frame polarization (imaginary part)
const auto &time = inputWs->x(0); // X axis: time
const auto &labRe = inputWs->y(0); // Lab-frame polarization (real part)
const auto &labIm = inputWs->y(1); // Lab-frame polarization (imaginary part)
MantidVec rrfRe(nData),
rrfIm(nData); // Rotating Reference frame (RRF) polarizations
for (size_t t = 0; t < nData; t++) {
......@@ -90,12 +90,12 @@ void RRFMuon::exec() {
// Put results into output workspace
// Real RRF polarization
outputWs->getSpectrum(0).setSpectrumNo(1);
outputWs->dataX(0) = inputWs->readX(0);
outputWs->dataY(0) = rrfRe;
outputWs->setSharedX(0, inputWs->sharedX(0));
outputWs->mutableY(0) = rrfRe;
// Imaginary RRF polarization
outputWs->getSpectrum(1).setSpectrumNo(2);
outputWs->dataX(1) = inputWs->readX(1);
outputWs->dataY(1) = rrfIm;
outputWs->setSharedX(1, inputWs->sharedX(1));
outputWs->mutableY(1) = rrfIm;
// Set output workspace
setProperty("OutputWorkspace", outputWs);
......
......@@ -83,11 +83,11 @@ void SANSDirectBeamScaling::exec() {
Progress progress(this, 0.0, 1.0, numHists);
// Number of X bins
const int64_t xLength = inputWS->readY(0).size();
const int64_t xLength = inputWS->y(0).size();
// Monitor counts
double monitor = 0.0;
const MantidVec &MonIn = inputWS->readY(index[0]);
const auto &MonIn = inputWS->y(index[0]);
for (int64_t j = 0; j < xLength; j++)
monitor += MonIn[j];
......@@ -111,8 +111,8 @@ void SANSDirectBeamScaling::exec() {
if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
continue;
const MantidVec &YIn = inputWS->readY(i);
const MantidVec &EIn = inputWS->readE(i);
const auto &YIn = inputWS->y(i);
const auto &EIn = inputWS->e(i);
// Sum up all the counts
V3D pos = spectrumInfo.position(i) - V3D(sourcePos.X(), sourcePos.Y(), 0.0);
......
......@@ -105,9 +105,6 @@ void MayersSampleCorrection::exec() {
for (int64_t i = 0; i < static_cast<int64_t>(nhist); ++i) {
PARALLEL_START_INTERUPT_REGION
// Copy the X values over
const auto &inX = inputWS->readX(i);
outputWS->dataX(i) = inX;
IDetector_const_sptr det;
try {
det = inputWS->getDetector(i);
......@@ -128,9 +125,9 @@ void MayersSampleCorrection::exec() {
params.sigmaSc = sampleMaterial.totalScatterXSection();
params.cylRadius = radius;
params.cylHeight = height;
MayersSampleCorrectionStrategy correction(params, inX, inputWS->readY(i),
inputWS->readE(i));
correction.apply(outputWS->dataY(i), outputWS->dataE(i));
MayersSampleCorrectionStrategy correction(params, inputWS->histogram(i));
outputWS->setHistogram(i, correction.getCorrectedHisto());
prog.report();
PARALLEL_END_INTERUPT_REGION
......
......@@ -2,10 +2,10 @@
// Includes
//-----------------------------------------------------------------------------
#include "MantidAlgorithms/SampleCorrections/MayersSampleCorrectionStrategy.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/Statistics.h"
#include "MantidKernel/Math/ChebyshevPolyFit.h"
#include "MantidKernel/Math/Distributions/ChebyshevSeries.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/Statistics.h"
#include <cassert>
#include <cmath>
......@@ -73,23 +73,18 @@ namespace Algorithms {
/**
* Constructor
* @param params Defines the required parameters for the correction
* @param tof The TOF values corresponding to the signals to correct.
* Number of tof values must match number of signal/error or be 1 greater
* @param sigIn Values of the signal that will be corrected
* @param errIn Values of the errors that will be corrected
* @param inputHist A histogram containing the TOF values corresponding
* to the values to be corrected.
*/
MayersSampleCorrectionStrategy::MayersSampleCorrectionStrategy(
MayersSampleCorrectionStrategy::Parameters params,
const std::vector<double> &tof, const std::vector<double> &sigIn,
const std::vector<double> &errIn)
: m_pars(params), m_tof(tof), m_sigin(sigIn), m_errin(errIn),
m_histogram(tof.size() == sigIn.size() + 1),
m_muRrange(calculateMuRange()), m_rng(new MersenneTwister(1)) {
// Sanity check
assert(sigIn.size() == tof.size() || sigIn.size() == tof.size() - 1);
assert(errIn.size() == tof.size() || sigIn.size() == tof.size() - 1);
const Mantid::HistogramData::Histogram &inputHist)
: m_pars(params), m_histogram(inputHist), m_tofVals(inputHist.points()),
m_histoYSize(inputHist.y().size()), m_muRrange(calculateMuRange()),
m_rng(new MersenneTwister(1)) {
if (!(m_tof.front() < m_tof.back())) {
const auto &xVals = m_histogram.x();
if (!(xVals.front() < xVals.back())) {
throw std::invalid_argument(
"TOF values are expected to be monotonically increasing");
}
......@@ -104,15 +99,11 @@ MayersSampleCorrectionStrategy::~MayersSampleCorrectionStrategy() = default;
* Correct the data for absorption and multiple scattering effects. Allows
* both histogram or point data. For histogram the TOF is taken to be
* the mid point of a bin
* @param sigOut Signal values to correct [In/Out]
* @param errOut Error values to correct [In/Out]
*
* @return A histogram containing corrected values
*/
void MayersSampleCorrectionStrategy::apply(std::vector<double> &sigOut,
std::vector<double> &errOut) {
const size_t nsig(m_sigin.size());
// Sanity check
assert(sigOut.size() == m_sigin.size());
assert(errOut.size() == m_errin.size());
Mantid::HistogramData::Histogram
MayersSampleCorrectionStrategy::getCorrectedHisto() {
// Temporary storage
std::vector<double> xmur(N_MUR_PTS + 1, 0.0),
......@@ -162,8 +153,13 @@ void MayersSampleCorrectionStrategy::apply(std::vector<double> &sigOut,
const double rns = (vol * 1e6) * (m_pars.rho * 1e24) * 1e-22;
ChebyshevSeries chebyPoly(N_POLY_ORDER);
for (size_t i = 0; i < nsig; ++i) {
const double sigt = sigmaTotal(flightPath, tof(i));
auto outputHistogram = m_histogram;
auto &sigOut = outputHistogram.mutableY();
auto &errOut = outputHistogram.mutableE();
for (size_t i = 0; i < m_histoYSize; ++i) {
const double sigt = sigmaTotal(flightPath, m_tofVals[i]);
const double rmu = muR(sigt);
// Varies between [-1,+1]
const double xcap = ((rmu - muMin) - (muMax - rmu)) / (muMax - muMin);
......@@ -174,10 +170,11 @@ void MayersSampleCorrectionStrategy::apply(std::vector<double> &sigOut,
corrfact *= (1.0 - beta) / rns;
}
// apply correction
const double yin(m_sigin[i]), ein(m_errin[i]);
const double yin(m_histogram.y()[i]), ein(m_histogram.e()[i]);
sigOut[i] = yin * corrfact;
errOut[i] = sigOut[i] * ein / yin;
}
return outputHistogram;
}
/**
......@@ -290,7 +287,7 @@ MayersSampleCorrectionStrategy::calculateMS(const size_t irp, const double muR,
std::pair<double, double>
MayersSampleCorrectionStrategy::calculateMuRange() const {
const double flightPath(m_pars.l1 + m_pars.l2);
const double tmin(tof(0)), tmax(tof(m_sigin.size() - 1));
const double tmin(m_tofVals[0]), tmax(m_tofVals[m_histoYSize - 1]);
return std::make_pair(muR(flightPath, tmin), muR(flightPath, tmax));
}
......@@ -330,17 +327,6 @@ double MayersSampleCorrectionStrategy::sigmaTotal(const double flightPath,
return sigabs + m_pars.sigmaSc;
}
/**
* Return the TOF for the given index of the signal value, taking into account
* if we have a histogram. Histograms will use the mid point of the bin
* as the TOF value. Note that there is no range check for the index.
* @param i Index of the signal value
* @return The associated TOF value
*/
double MayersSampleCorrectionStrategy::tof(const size_t i) const {
return m_histogram ? 0.5 * (m_tof[i] + m_tof[i + 1]) : m_tof[i];
}
/**
* (Re-)seed the random number generator
* @param seed Seed value for the random number generator
......
......@@ -116,25 +116,24 @@ void SetUncertainties::exec() {
PARALLEL_START_INTERUPT_REGION
// copy the X/Y
outputWorkspace->setX(i, inputWorkspace->refX(i));
outputWorkspace->dataY(i) = inputWorkspace->readY(i);
outputWorkspace->setSharedX(i, inputWorkspace->sharedX(i));
outputWorkspace->setSharedY(i, inputWorkspace->sharedY(i));
// copy the E or set to zero depending on the mode
if (errorType.compare(ONE_IF_ZERO) == 0) {
outputWorkspace->dataE(i) = inputWorkspace->readE(i);
outputWorkspace->setSharedE(i, inputWorkspace->sharedE(i));
} else {
outputWorkspace->dataE(i) =
std::vector<double>(inputWorkspace->readE(i).size(), 0.);
outputWorkspace->mutableE(i) = 0.0;
}
// ZERO mode doesn't calculate anything further
if ((!zeroError) && (!isMasked(inputWorkspace, i))) {
MantidVec &E = outputWorkspace->dataE(i);
auto &E = outputWorkspace->mutableE(i);
if (takeSqrt) {
const MantidVec &Y = outputWorkspace->readY(i);
const auto &Y = outputWorkspace->y(i);
std::transform(Y.begin(), Y.end(), E.begin(),
sqrterror(resetOne ? 1. : 0.));
} else {
std::for_each(E.begin(), E.end(), resetzeroerror(1.));
std::transform(E.begin(), E.end(), E.begin(), resetzeroerror(1.));
}
}
......
......@@ -92,19 +92,19 @@ void SmoothData::exec() {
// Copy the X data over. Preserves data sharing if present in input
// workspace.
outputWorkspace->setX(i, inputWorkspace->refX(i));
outputWorkspace->setSharedX(i, inputWorkspace->sharedX(i));
// Now get references to the Y & E vectors in the input and output
// workspaces
const MantidVec &Y = inputWorkspace->readY(i);
const MantidVec &E = inputWorkspace->readE(i);
MantidVec &newY = outputWorkspace->dataY(i);
MantidVec &newE = outputWorkspace->dataE(i);
if (npts == 0) {
newY = Y;
newE = E;
outputWorkspace->setSharedY(i, inputWorkspace->sharedY(i));
outputWorkspace->setSharedE(i, inputWorkspace->sharedE(i));
continue;
}
const auto &Y = inputWorkspace->y(i);
const auto &E = inputWorkspace->e(i);
auto &newY = outputWorkspace->mutableY(i);
auto &newE = outputWorkspace->mutableE(i);
// Use total to help hold our moving average
double total = 0.0, totalE = 0.0;
// First push the values ahead of the current point onto total
......
......@@ -4,10 +4,12 @@
#include <cxxtest/TestSuite.h>
#include "MantidAlgorithms/SampleCorrections/MayersSampleCorrectionStrategy.h"
#include "MantidHistogramData/LinearGenerator.h"
#include <algorithm>
#include <cmath>
using Mantid::Algorithms::MayersSampleCorrectionStrategy;
using namespace Mantid::HistogramData;
class MayersSampleCorrectionStrategyTest : public CxxTest::TestSuite {
public:
......@@ -21,10 +23,8 @@ public:
}
void test_Attentuaton_Correction_For_Fixed_Mur() {
std::vector<double> dummy(2, 0.0);
dummy[1] = 1.0;
MayersSampleCorrectionStrategy mscat(createTestParameters(), dummy, dummy,
dummy);
Histogram histo(Points{0, 1}, Counts{0, 1});
MayersSampleCorrectionStrategy mscat(createTestParameters(), histo);
auto absFactor = mscat.calculateSelfAttenuation(0.01);
const double delta = 1e-8;
......@@ -33,10 +33,8 @@ public:
void
test_Multiple_Scattering_With_Fixed_Mur_And_Absorption_Correction_Factor() {
std::vector<double> dummy(2, 0.0);
dummy[1] = 1.0;
MayersSampleCorrectionStrategy mscat(createTestParameters(), dummy, dummy,
dummy);
Histogram histo(Points{0, 1}, Counts{0, 1});
MayersSampleCorrectionStrategy mscat(createTestParameters(), histo);
const size_t irp(1);
const double muR(0.01), abs(0.0003);
auto absFactor = mscat.calculateMS(irp, muR, abs);
......@@ -47,17 +45,16 @@ public:
}
void test_Corrects_Both_Absorption_And_Multiple_Scattering_For_Point_Data() {
using std::sqrt;
const size_t nypts(100);
std::vector<double> signal(nypts, 2.0), tof(nypts), error(nypts);
std::transform(signal.begin(), signal.end(), error.begin(),
(double (*)(double))sqrt);
std::generate(tof.begin(), tof.end(), Incrementer(100.0));
MayersSampleCorrectionStrategy mscat(createTestParameters(), tof, signal,
error);
Histogram histo(Points(nypts, LinearGenerator(100.0, 1.0)),
Counts(nypts, 2.0));
MayersSampleCorrectionStrategy mscat(createTestParameters(), histo);
// Correct it
mscat.apply(signal, error);
auto outHisto = mscat.getCorrectedHisto();
const auto &tof = outHisto.x();
const auto &signal = outHisto.y();
const auto &error = outHisto.e();
// Check some values
const double delta(1e-06);
......@@ -73,18 +70,16 @@ public:
void
test_Corrects_Both_Absorption_And_Multiple_Scattering_For_Histogram_Data() {
using std::sqrt;
const size_t nypts(100);
std::vector<double> signal(nypts, 2.0), tof(nypts + 1), error(nypts);
std::transform(signal.begin(), signal.end(), error.begin(),
(double (*)(double))sqrt);
// Generate a histogram with the same mid points as the point data example
std::generate(tof.begin(), tof.end(), Incrementer(99.5));
MayersSampleCorrectionStrategy mscat(createTestParameters(), tof, signal,
error);
Histogram histo(BinEdges(nypts + 1, LinearGenerator(99.5, 1.0)),
Counts(nypts, 2.0));
MayersSampleCorrectionStrategy mscat(createTestParameters(), histo);
auto outHisto = mscat.getCorrectedHisto();
// Correct it
mscat.apply(signal, error);
const auto &tof = outHisto.x();
const auto &signal = outHisto.y();
const auto &error = outHisto.e();
// Check some values
const double delta(1e-06);
......@@ -99,19 +94,17 @@ public:
}
void test_Corrects_For_Absorption_For_Histogram_Data() {
using std::sqrt;
const size_t nypts(100);
std::vector<double> signal(nypts, 2.0), tof(nypts + 1), error(nypts);
std::transform(signal.begin(), signal.end(), error.begin(),
(double (*)(double))sqrt);
// Generate a histogram with the same mid points as the point data example
std::generate(tof.begin(), tof.end(), Incrementer(99.5));
bool mscatOn(false);
MayersSampleCorrectionStrategy mscat(createTestParameters(mscatOn), tof,
signal, error);
Histogram histo(BinEdges(nypts + 1, LinearGenerator(99.5, 1.0)),
Counts(nypts, 2.0));
MayersSampleCorrectionStrategy mscat(createTestParameters(mscatOn), histo);
auto outHisto = mscat.getCorrectedHisto();
// Correct it
mscat.apply(signal, error);
auto tof = outHisto.x();
auto signal = outHisto.y();
auto error = outHisto.e();
// Check some values
const double delta(1e-06);
......@@ -127,29 +120,16 @@ public:
// ---------------------- Failure tests -----------------------------
void test_Tof_Not_Monotonically_Increasing_Throws_Invalid_Argument() {
using std::sqrt;
const size_t nypts(10);
std::vector<double> signal(nypts, 2.0), tof(nypts + 1), error(nypts);
std::transform(signal.begin(), signal.end(), error.begin(),
(double (*)(double))sqrt);
std::generate(tof.begin(), tof.end(), Decrementer(199.5));
TS_ASSERT_THROWS(MayersSampleCorrectionStrategy(createTestParameters(), tof,
signal, error),
std::invalid_argument);
Histogram histo(BinEdges(nypts + 1, LinearGenerator(199.5, -1.0)),
Counts(nypts, 2.0));
TS_ASSERT_THROWS(
MayersSampleCorrectionStrategy(createTestParameters(), histo),
std::invalid_argument);
}
private:
struct Incrementer {
Incrementer(double start) : current(start) {}
double operator()() { return current++; }
double current;
};
struct Decrementer {
Decrementer(double start) : current(start) {}
double operator()() { return current--; }
double current;
};
MayersSampleCorrectionStrategy::Parameters
createTestParameters(bool mscatOn = true) {
// A bit like a POLARIS spectrum
......
......@@ -34,9 +34,9 @@ public:
TS_ASSERT(alg->isExecuted());
MatrixWorkspace_sptr corrected = alg->getProperty("OutputWorkspace");
const auto &tof = corrected->readX(0);
const auto &signal = corrected->readY(0);
const auto &error = corrected->readE(0);
const auto &tof = corrected->x(0);
const auto &signal = corrected->y(0);
const auto &error = corrected->e(0);
const double delta(1e-06);
TS_ASSERT_DELTA(99.5, tof.front(), delta);
TS_ASSERT_DELTA(199.5, tof.back(), delta);
......@@ -57,9 +57,9 @@ public:
TS_ASSERT(alg->isExecuted());
MatrixWorkspace_sptr corrected = alg->getProperty("OutputWorkspace");
const auto &tof = corrected->readX(0);
const auto &signal = corrected->readY(0);
const auto &error = corrected->readE(0);
const auto &tof = corrected->x(0);
const auto &signal = corrected->y(0);
const auto &error = corrected->e(0);
const double delta(1e-06);
TS_ASSERT_DELTA(99.5, tof.front(), delta);
TS_ASSERT_DELTA(199.5, tof.back(), delta);
......
......@@ -7,6 +7,7 @@
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidTestHelpers/HistogramDataTestHelper.h"
using namespace Mantid::Algorithms;
using namespace Mantid::API;
......@@ -46,11 +47,11 @@ public:
// Checks
// X values
TS_ASSERT_EQUALS(ws->readX(0), ows->readX(0));
TS_ASSERT_EQUALS(ws->readX(1), ows->readX(1));
TS_ASSERT_EQUALS(ws->x(0), ows->x(0));
TS_ASSERT_EQUALS(ws->x(1), ows->x(1));
// Y values
TS_ASSERT_EQUALS(ws->readY(0), ows->readY(0));
TS_ASSERT_EQUALS(ws->readY(1), ows->readY(1));
TS_ASSERT_EQUALS(ws->y(0), ows->y(0));
TS_ASSERT_EQUALS(ws->y(1), ows->y(1));
}
void testRRFMuonNonZeroFrequency() {
......@@ -80,19 +81,19 @@ public:
// Checks
// X values
TS_ASSERT_EQUALS(ws->readX(0), ows->readX(0));
TS_ASSERT_EQUALS(ws->readX(1), ows->readX(1));
TS_ASSERT_EQUALS(ws->x(0), ows->x(0));
TS_ASSERT_EQUALS(ws->x(1), ows->x(1));
// Y values
// The input frequency is close to the precession frequency, so:
// The real part of the RRF polarization should be close to 1 for all X
// values
// The imaginary part should be close to 0 for all X values
TS_ASSERT_DELTA(ows->readY(0)[0], 1, 0.001);
TS_ASSERT_DELTA(ows->readY(0)[100], 1, 0.001);
TS_ASSERT_DELTA(ows->readY(0)[200], 1, 0.001);
TS_ASSERT_DELTA(ows->readY(1)[0], 0, 0.001);
TS_ASSERT_DELTA(ows->readY(1)[100], 0, 0.001);
TS_ASSERT_DELTA(ows->readY(1)[200], 0, 0.001);
TS_ASSERT_DELTA(ows->y(0)[0], 1, 0.001);
TS_ASSERT_DELTA(ows->y(0)[100], 1, 0.001);
TS_ASSERT_DELTA(ows->y(0)[200], 1, 0.001);
TS_ASSERT_DELTA(ows->y(1)[0], 0, 0.001);
TS_ASSERT_DELTA(ows->y(1)[100], 0, 0.001);
TS_ASSERT_DELTA(ows->y(1)[200], 0, 0.001);
}
void testRRFMuonUnits() {
......@@ -154,27 +155,27 @@ public:
// Check Y values
// ows1 vs ows2
// Results with different frequency units should be very similar
TS_ASSERT_DELTA(ows1->readY(0)[5], ows2->readY(0)[5], 0.000001);
TS_ASSERT_DELTA(ows1->readY(0)[98], ows2->readY(0)[98], 0.000001);
TS_ASSERT_DELTA(ows1->readY(0)[276], ows2->readY(0)[276], 0.000001);
TS_ASSERT_DELTA(ows1->y(0)[5], ows2->y(0)[5], 0.000001);
TS_ASSERT_DELTA(ows1->y(0)[98], ows2->y(0)[98], 0.000001);
TS_ASSERT_DELTA(ows1->y(0)[276], ows2->y(0)[276], 0.000001);
// But not exactly the same
// (They should only be the same if the input frequency in rrfMuon2 were
// exactly 1/2/M_PI)
TS_ASSERT_DIFFERS(ows1->readY(0)[5], ows2->readY(0)[5]);
TS_ASSERT_DIFFERS(ows1->readY(0)[98], ows2->readY(0)[98]);