Commit 1f27ed92 authored by Alex Buts's avatar Alex Buts
Browse files

Re #16049 Unit test for GetEi v1 and minor changes to make it

more GetEi V2 compatible.
parent d33b7adb
......@@ -708,6 +708,7 @@ set ( TEST_FILES
GetDetOffsetsMultiPeaksTest.h
GetDetectorOffsetsTest.h
GetEiTest.h
GetEiV1Test.h
GetTimeSeriesLogInformationTest.h
GroupWorkspacesTest.h
HRPDSlabCanAbsorptionTest.h
......
......@@ -91,8 +91,8 @@ private:
specnum_t specNum2) const;
double timeToFly(double s, double E_KE) const;
double getPeakCentre(API::MatrixWorkspace_const_sptr WS,
const int64_t monitIn, const double peakTime);
void extractSpec(int64_t wsInd, double start, double end);
const size_t monitIn, const double peakTime);
void extractSpec(int wsInd, double start, double end);
void getPeakEstimates(double &height, int64_t &centreInd,
double &background) const;
double findHalfLoc(MantidVec::size_type startInd, const double height,
......
......@@ -7,6 +7,7 @@
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/Exception.h"
#include <boost/lexical_cast.hpp>
#include <cmath>
......@@ -81,7 +82,7 @@ void GetEi::init() {
* does not have common binning
*/
void GetEi::exec() {
MatrixWorkspace_const_sptr inWS = getProperty("InputWorkspace");
MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
const specnum_t mon1Spec = getProperty("Monitor1Spec");
const specnum_t mon2Spec = getProperty("Monitor2Spec");
double dist2moni0 = -1, dist2moni1 = -1;
......@@ -139,6 +140,11 @@ void GetEi::exec() {
<< " (your estimate was " << E_est << " meV)\n";
setProperty("IncidentEnergy", E_i);
// store property in input workspace
Property *incident_energy =
new PropertyWithValue<double>("Ei", E_i, Direction::Input);
inWS->mutableRun().addProperty(incident_energy, true);
}
/** Gets the distances between the source and detectors whose IDs you pass to it
* @param WS :: the input workspace
......@@ -273,14 +279,23 @@ double GetEi::timeToFly(double s, double E_KE) const {
* @throw runtime_error a Child Algorithm just falls over
*/
double GetEi::getPeakCentre(API::MatrixWorkspace_const_sptr WS,
const int64_t monitIn, const double peakTime) {
const size_t monitIn, const double peakTime) {
const MantidVec &timesArray = WS->readX(monitIn);
// we search for the peak only inside some window because there are often more
// peaks in the monitor histogram
double halfWin = (timesArray.back() - timesArray.front()) * HALF_WINDOW;
if (monitIn < std::numeric_limits<int>::max()) {
int ivsInd = static_cast<int>(monitIn);
// runs CropWorkspace as a Child Algorithm to and puts the result in a new
// temporary workspace that will be deleted when this algorithm has finished
extractSpec(monitIn, peakTime - halfWin, peakTime + halfWin);
extractSpec(ivsInd, peakTime - halfWin, peakTime + halfWin);
}
else {
throw Kernel::Exception::NotImplementedError("Spectra number exceeds maximal"
" integer number defined for this OS."
" This behaviour is not yet supported");
}
// converting the workspace to count rate is required by the fitting algorithm
// if the bin widths are not all the same
WorkspaceHelpers::makeDistribution(m_tempWS);
......@@ -318,7 +333,7 @@ double GetEi::getPeakCentre(API::MatrixWorkspace_const_sptr WS,
* @throw runtime_error if the algorithm just falls over
* @throw invalid_argument if the input workspace does not have common binning
*/
void GetEi::extractSpec(int64_t wsInd, double start, double end) {
void GetEi::extractSpec(int wsInd, double start, double end) {
IAlgorithm_sptr childAlg = createChildAlgorithm(
"CropWorkspace", 100 * m_fracCompl, 100 * (m_fracCompl + CROP));
m_fracCompl += CROP;
......@@ -328,14 +343,8 @@ void GetEi::extractSpec(int64_t wsInd, double start, double end) {
childAlg->setProperty("XMin", start);
childAlg->setProperty("XMax", end);
if (wsInd < std::numeric_limits<int>::max()) {
auto ivsInd = static_cast<int>(wsInd);
childAlg->setProperty("StartWorkspaceIndex", ivsInd);
childAlg->setProperty("EndWorkspaceIndex", ivsInd);
} else {
childAlg->setProperty("StartWorkspaceIndex", wsInd);
childAlg->setProperty("EndWorkspaceIndex", wsInd);
}
childAlg->executeAsChildAlg();
m_tempWS = childAlg->getProperty("OutputWorkspace");
......
#ifndef GETEIV1TEST_H_
#define GETEIV1TEST_H_
#include "MantidAPI/Axis.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/FrameworkManager.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
#include <cxxtest/TestSuite.h>
#include <cmath>
using namespace Mantid::Kernel;
using namespace Mantid::API;
using Mantid::MantidVecPtr;
class GetEiV1Test : 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 GetEiV1Test *createSuite() { return new GetEiV1Test(); }
static void destroySuite(GetEiV1Test *suite) { delete suite; }
GetEiV1Test() {
FrameworkManager::Instance(); // Load plugins
}
void test_Result_For_Good_Estimate() {
const double input_ei = 15.0;
const bool fixei = false;
do_test_on_result_values(input_ei, fixei);
}
void do_test_on_result_values(double input_ei, bool fixei) {
MatrixWorkspace_sptr testWS = createTestWorkspaceWithMonitors();
// This algorithm needs a name attached to the workspace
const std::string outputName("eitest");
AnalysisDataService::Instance().add(outputName, testWS);
IAlgorithm_sptr alg;
TS_ASSERT_THROWS_NOTHING(
alg = runGetEiUsingTestMonitors(outputName, input_ei, fixei));
// Test output answers
// The monitor peak should always be calculated from the data
const double expected_mon_peak = 6495.7499801169;
const int expected_mon_index = 0;
const double expected_ei = (fixei) ? input_ei : 15.001453367;
const double ei = alg->getProperty("IncidentEnergy");
const double first_mon_peak = alg->getProperty("FirstMonitorPeak");
TS_ASSERT_DELTA(ei, expected_ei, 1e-08);
TS_ASSERT_DELTA(first_mon_peak, expected_mon_peak, 1e-08);
// and verify it has been store on the run object
Property *ei_runprop = testWS->run().getProperty("Ei");
PropertyWithValue<double> *ei_propvalue =
dynamic_cast<PropertyWithValue<double> *>(ei_runprop);
TS_ASSERT_DELTA((*ei_propvalue)(), expected_ei, 1e-08);
AnalysisDataService::Instance().remove(outputName);
}
void testParametersOnWorkspace() {
MatrixWorkspace_sptr testWS = createTestWorkspaceWithMonitors();
testWS->instrumentParameters().addString(
testWS->getInstrument()->getChild(0).get(), "ei-mon1-spec", "1");
testWS->instrumentParameters().addString(
testWS->getInstrument()->getChild(0).get(), "ei-mon2-spec", "2");
Property *incident_energy_guess =
new PropertyWithValue<double>("EnergyRequest", 15.0, Direction::Input);
testWS->mutableRun().addProperty(incident_energy_guess, true);
// This algorithm needs a name attached to the workspace
const std::string outputName("eiNoParTest");
AnalysisDataService::Instance().add(outputName, testWS);
IAlgorithm_sptr alg;
TS_ASSERT_THROWS_NOTHING(
alg = runGetEiUsingTestMonitors(outputName, 15, false));
// Test output answers
const double expected_ei = 15.001453367;
const double ei = alg->getProperty("IncidentEnergy");
TS_ASSERT_DELTA(ei, expected_ei, 1e-08);
// and verify it has been store on the run object
Property *ei_runprop = testWS->run().getProperty("Ei");
PropertyWithValue<double> *ei_propvalue =
dynamic_cast<PropertyWithValue<double> *>(ei_runprop);
TS_ASSERT_DELTA((*ei_propvalue)(), expected_ei, 1e-08);
AnalysisDataService::Instance().remove(outputName);
}
void testThrowsMon1() {
MatrixWorkspace_sptr testWS = createTestWorkspaceWithMonitors();
// This algorithm needs a name attached to the workspace
const std::string outputName("eitest1");
AnalysisDataService::Instance().add(outputName, testWS);
IAlgorithm_sptr alg =
AlgorithmManager::Instance().createUnmanaged("GetEi", 1);
alg->initialize();
alg->setPropertyValue("InputWorkspace", outputName);
alg->setProperty("Monitor2Spec", 2);
alg->setProperty("EnergyEstimate", 15.0);
alg->setRethrows(true);
TS_ASSERT_THROWS_EQUALS(
alg->execute(), const std::runtime_error &e, std::string(e.what()),
"Some invalid Properties found");
AnalysisDataService::Instance().remove(outputName);
}
void testThrowsEi() {
MatrixWorkspace_sptr testWS = createTestWorkspaceWithMonitors();
// This algorithm needs a name attached to the workspace
const std::string outputName("eitest2");
AnalysisDataService::Instance().add(outputName, testWS);
IAlgorithm_sptr alg =
AlgorithmManager::Instance().createUnmanaged("GetEi", 1);
alg->initialize();
alg->setPropertyValue("InputWorkspace", outputName);
alg->setProperty("Monitor1Spec", 1);
alg->setProperty("Monitor2Spec", 2);
alg->setRethrows(true);
TS_ASSERT_THROWS_EQUALS(alg->execute(), const std::runtime_error &e,
std::string(e.what()),
"Some invalid Properties found");
AnalysisDataService::Instance().remove(outputName);
}
void test_throws_error_when_ei_not_fixed_and_no_peaks_found() {
const bool includePeaks(false);
MatrixWorkspace_sptr testWS = createTestWorkspaceWithMonitors(includePeaks);
// This algorithm needs a name attached to the workspace
const std::string outputName("eitest2");
AnalysisDataService::Instance().add(outputName, testWS);
IAlgorithm_sptr alg =
AlgorithmManager::Instance().createUnmanaged("GetEi", 1);
alg->initialize();
alg->setPropertyValue("InputWorkspace", outputName);
alg->setProperty("Monitor1Spec", 1);
alg->setProperty("Monitor2Spec", 2);
alg->setProperty("EnergyEstimate", 15.0);
alg->setRethrows(true);
TS_ASSERT_THROWS(alg->execute(), std::invalid_argument);
AnalysisDataService::Instance().remove(outputName);
}
private:
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");
MantidVecPtr xdata;
xdata.access().resize(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(150. * 150.),
peakTwoCentre(10625.), sigmaSqTwo(25. * 25.);
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->dataY(0)[i] =
peakOneHeight *
exp(-0.5 * pow(xValue - peakOneCentre, 2.) / sigmaSqOne);
testWS->dataY(1)[i] =
peakTwoHeight *
exp(-0.5 * pow(xValue - peakTwoCentre, 2.) / sigmaSqTwo);
}
xdata.access()[i] = xValue;
}
testWS->setX(0, xdata);
testWS->setX(1, xdata);
return testWS;
}
IAlgorithm_sptr runGetEiUsingTestMonitors(const std::string &inputWS,
const double energyGuess,
const bool fixei) {
IAlgorithm_sptr alg =
AlgorithmManager::Instance().createUnmanaged("GetEi", 1);
alg->initialize();
alg->setPropertyValue("InputWorkspace", inputWS);
alg->setProperty("Monitor1Spec", 1);
alg->setProperty("Monitor2Spec", 2);
alg->setProperty("EnergyEstimate", energyGuess);
alg->setRethrows(true);
alg->execute();
return alg;
}
};
#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