Skip to content
GitLab
Commit adebbce4 authored by Dimitar Tasev's avatar Dimitar Tasev
Browse files

Re #16959 refactore GetEi and test

parent 9ce87603
Hide whitespace changes
Inline Side-by-side
Framework/Algorithms/src/GetEi.cpp
View file @ adebbce4
......@@ -7,7 +7,6 @@
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/Exception.h"
#include <boost/lexical_cast.hpp>
#include <cmath>
......@@ -367,22 +366,26 @@ void GetEi::extractSpec(int wsInd, double start, double end) {
*/
void GetEi::getPeakEstimates(double &height, int64_t &centreInd,
double &background) const {
const auto &X = m_tempWS->x(0);
const auto &Y = m_tempWS->y(0);
// take note of the number of background counts as error checking, do we have
// a peak or just a bump in the background
background = 0;
// start at the first Y value
height = m_tempWS->y(0)[0];
height = Y[0];
centreInd = 0;
// then loop through all the Y values and find the tallest peak
for (size_t i = 1; i < m_tempWS->y(0).size() - 1; ++i) {
background += m_tempWS->y(0)[i];
if (m_tempWS->y(0)[i] > height) {
for (size_t i = 1; i < Y.size() - 1; ++i) {
background += Y[i];
if (Y[i] > height) {
centreInd = i;
height = m_tempWS->y(0)[centreInd];
height = Y[centreInd];
}
}
background = background / static_cast<double>(m_tempWS->y(0).size());
background = background / static_cast<double>(Y.size());
if (height < PEAK_THRESH_H * background) {
throw std::invalid_argument(
"No peak was found or its height is less than the threshold of " +
......@@ -393,10 +396,8 @@ void GetEi::getPeakEstimates(double &height, int64_t &centreInd,
g_log.debug() << "Peak position is the bin that has the maximum Y value in "
"the monitor spectrum, which is at TOF "
<< (m_tempWS->x(0)[centreInd] +
m_tempWS->x(0)[centreInd + 1]) /
2 << " (peak height " << height
<< " counts/microsecond)\n";
<< (X[centreInd] + X[centreInd + 1]) / 2 << " (peak height "
<< height << " counts/microsecond)\n";
}
/** Estimates the closest time, looking either or back, when the number of
* counts is
......@@ -416,8 +417,9 @@ double GetEi::findHalfLoc(size_t startInd, const double height,
const double noise, const direction go) const {
auto endInd = startInd;
// todo move y(0) to const auto
while (m_tempWS->y(0)[endInd] > (height + noise) / 2.0) {
const auto &X = m_tempWS->x(0);
const auto &Y = m_tempWS->y(0);
while (Y[endInd] > (height + noise) / 2.0) {
endInd += go;
if (endInd < 1) {
throw std::out_of_range(
......@@ -426,7 +428,7 @@ double GetEi::findHalfLoc(size_t startInd, const double height,
"% window, at TOF values that are too low. Was the energy estimate "
"close enough?");
}
if (endInd > m_tempWS->y(0).size() - 2) {
if (endInd > Y.size() - 2) {
throw std::out_of_range(
"Can't analyse peak, some of the peak is outside the " +
boost::lexical_cast<std::string>(HALF_WINDOW * 100) +
......@@ -458,8 +460,7 @@ double GetEi::findHalfLoc(size_t startInd, const double height,
}
// get the TOF value in the middle of the bin with the first value below the
// half height
double halfTime =
(m_tempWS->x(0)[endInd] + m_tempWS->x(0)[endInd + 1]) / 2;
double halfTime = (X[endInd] + X[endInd + 1]) / 2;
// interpolate back between the first bin with less than half the counts to
// the bin before it
if (endInd != startInd) { // let the bin that we found have coordinates (x_1,
......@@ -468,11 +469,10 @@ double GetEi::findHalfLoc(size_t startInd, const double height,
// (y_3-y_1)/(x_3-x_1) where (x_3, y_3) are the
// coordinates of the other bin we are using
double gradient =
(m_tempWS->y(0)[endInd] - m_tempWS->y(0)[endInd - go]) /
(m_tempWS->x(0)[endInd] - m_tempWS->x(0)[endInd - go]);
(Y[endInd] - Y[endInd - go]) / (X[endInd] - X[endInd - go]);
// we don't need to check for a zero or negative gradient if we assume the
// endInd bin was found when the Y-value dropped below the threshold
double deltaY = m_tempWS->y(0)[endInd] - (height + noise) / 2.0;
double deltaY = Y[endInd] - (height + noise) / 2.0;
// correct for the interpolation back in the direction towards the peak
// centre
halfTime -= deltaY / gradient;
......
Framework/Algorithms/test/GetEiTest.h
View file @ adebbce4
......@@ -12,7 +12,6 @@
using namespace Mantid::Kernel;
using namespace Mantid::API;
using Mantid::MantidVecPtr;
using Mantid::HistogramData::BinEdges;
class GetEiTest : public CxxTest::TestSuite {
......@@ -241,10 +240,10 @@ private:
for (int i = 0; i <= numBins; ++i) {
const double xValue = 5.0 + 5.5 * i;
if (includePeaks && i < numBins) {
testWS->dataY(0)[i] =
testWS->mutableY(0)[i] =
peakOneHeight *
exp(-0.5 * pow(xValue - peakOneCentre, 2.) / sigmaSqOne);
testWS->dataY(1)[i] =
testWS->mutableY(1)[i] =
peakTwoHeight *
exp(-0.5 * pow(xValue - peakTwoCentre, 2.) / sigmaSqTwo);
}
......@@ -271,5 +270,113 @@ private:
return alg;
}
};
class GetEiTestPerformance : public CxxTest::TestSuite {
public:
// This pair of boilerplate methods prevent the suite being created statically
// This means the constructor isn't called when running other tests
static GetEiTestPerformance *createSuite() {
return new GetEiTestPerformance();
}
static void destroySuite(GetEiTestPerformance *suite) { delete suite; }
void setUp() {
inputWS1 = createTestWorkspaceWithMonitors();
outputName1 = "eitest1";
AnalysisDataService::Instance().add(outputName1, inputWS1);
inputWS2 = createTestWorkspaceWithMonitors();
outputName2 = "eitest2";
AnalysisDataService::Instance().add(outputName2, inputWS1);
// This algorithm needs a name attached to the workspace
}
void tearDown() {
AnalysisDataService::Instance().remove(outputName1);
AnalysisDataService::Instance().remove(outputName2);
}
void test_Result_For_Good_Estimate() {
const double input_ei = 15.0;
const bool fixei = false;
IAlgorithm_sptr alg;
TS_ASSERT_THROWS_NOTHING(
alg = runGetEiUsingTestMonitors(outputName1, input_ei, fixei));
}
void test_Result_When_Fixing_Ei() {
const double input_ei = 15.0;
const bool fixei = true;
IAlgorithm_sptr alg;
TS_ASSERT_THROWS_NOTHING(
alg = runGetEiUsingTestMonitors(outputName2, input_ei, fixei));
}
private:
/// nearly the same method as the unit test above
/// changed from std::string to MatrixWorkspace_sptr
IAlgorithm_sptr runGetEiUsingTestMonitors(const std::string inWSName,
const double energyGuess,
const bool fixei) {
IAlgorithm_sptr alg =
AlgorithmManager::Instance().createUnmanaged("GetEi", 2);
alg->initialize();
alg->setPropertyValue("InputWorkspace", inWSName);
alg->setProperty("Monitor1Spec", 1);
alg->setProperty("Monitor2Spec", 2);
alg->setProperty("FixEi", fixei);
alg->setProperty("EnergyEstimate", energyGuess);
alg->setRethrows(true);
alg->execute();
return alg;
}
/// Same method as the unit test above
MatrixWorkspace_sptr
createTestWorkspaceWithMonitors(const bool includePeaks = true) {
const int numHists(2);
const int numBins(2000);
MatrixWorkspace_sptr testWS =
WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(
numHists, numBins, true);
testWS->getAxis(0)->unit() =
Mantid::Kernel::UnitFactory::Instance().create("TOF");
BinEdges xdata(numBins + 1);
// Update X data to a sensible values. Looks roughly like the MARI binning
// Update the Y values. We don't care about errors here
// Instrument geometry + incident energy of ~15 mev (purely made up) gives
// these neceesary peak values.
// We'll simply use a gaussian as a test
const double peakOneCentre(6493.0), sigmaSqOne(250 * 250.),
peakTwoCentre(10625.), sigmaSqTwo(50 * 50);
const double peakOneHeight(3000.), peakTwoHeight(1000.);
for (int i = 0; i <= numBins; ++i) {
const double xValue = 5.0 + 5.5 * i;
if (includePeaks && i < numBins) {
testWS->mutableY(0)[i] =
peakOneHeight *
exp(-0.5 * pow(xValue - peakOneCentre, 2.) / sigmaSqOne);
testWS->mutableY(1)[i] =
peakTwoHeight *
exp(-0.5 * pow(xValue - peakTwoCentre, 2.) / sigmaSqTwo);
}
xdata.mutableData()[i] = xValue;
}
testWS->setBinEdges(0, xdata);
testWS->setBinEdges(1, xdata);
return testWS;
}
MatrixWorkspace_sptr inputWS1;
MatrixWorkspace_sptr inputWS2;
std::string outputName1;
std::string outputName2;
};
#endif /*GETEITEST_H_*/
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment