Skip to content

Commit

Permalink
GetEi2 throws if ei not fixed & peaks are not found.
Browse files Browse the repository at this point in the history 
Also added extra tests around these scenarios. Refs #6782
  • Loading branch information
@martyngigg
martyngigg committed Mar 27, 2013
1 parent d864a25 commit 864ebaa
Show file tree
Hide file tree
Showing 2 changed files with 134 additions and 98 deletions.
31 changes: 11 additions & 20 deletions Code/Mantid/Framework/Algorithms/src/GetEi2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,10 +129,6 @@ void GetEi2::exec()
}
}
double incident_energy = calculateEi(initial_guess);
if( !m_fixedei )
{
g_log.notice() << "Incident energy = " << incident_energy << " meV from initial guess = " << initial_guess << " meV\n";
}

storeEi(incident_energy);

Expand Down Expand Up @@ -240,39 +236,34 @@ double GetEi2::calculateEi(const double initial_guess)
}
catch(std::invalid_argument &)
{
if(!m_fixedei) throw;
peak_times[i] = 0.0;
g_log.information() << "No peak found for monitor " << (i+1) << " (at " << det_distances[i] << " metres). Setting peak time to zero\n";
}

if( i == 0 )
{
// Store for later adjustment of bins
if(i == 0)
{
m_peak1_pos = std::make_pair(static_cast<int>(ws_index), peak_times[i]);
if(m_fixedei)
{
g_log.information() << "Using fixed Ei=" << initial_guess << " meV\n";
break;
}
if(m_fixedei) break;
}
}

if( m_fixedei )
{
g_log.notice() << "Incident energy fixed at Ei=" << initial_guess << " meV\n";
return initial_guess;
}
else
{
double mean_speed = (det_distances[1] - det_distances[0])/(peak_times[1] - peak_times[0]);

double tzero = peak_times[1] - ((1.0/mean_speed)*det_distances[1]);
setProperty("Tzero", tzero);

g_log.debug() << "T0 = " << tzero << std::endl;

g_log.debug() << "Mean Speed = " << mean_speed << std::endl;
g_log.debug() << "Energy (meV) = " << mean_speed*mean_speed*m_t_to_mev << std::endl;

return mean_speed*mean_speed*m_t_to_mev;
setProperty("Tzero", tzero);

const double energy = mean_speed*mean_speed*m_t_to_mev;
g_log.notice() << "Incident energy calculated at Ei= " << energy
<< " meV from initial guess = " << initial_guess << " meV\n";
return energy;
}
}

Expand Down
201 changes: 123 additions & 78 deletions Code/Mantid/Framework/Algorithms/test/GetEiTest.h
View file
Original file line number Diff line number Diff line change
@@ -1,24 +1,30 @@
#ifndef GETE_ITEST_H_
#define GETE_ITEST_H_
#ifndef GETEITEST_H_
#define GETEITEST_H_

#include <cxxtest/TestSuite.h>
#include "MantidAlgorithms/GetEi2.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/FrameworkManager.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
#include <cxxtest/TestSuite.h>
#include <cmath>
#include "MantidDataHandling/LoadNexusMonitors.h"
#include "MantidDataHandling/SaveNexus.h"


using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Algorithms;
using Mantid::DataObjects::Workspace2D_sptr;
using Mantid::MantidVecPtr;

class GetEiTest : 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 GetEiTest *createSuite() { return new GetEiTest(); }
static void destroySuite( GetEiTest *suite ) { delete suite; }

GetEiTest()
{
FrameworkManager::Instance(); // Load plugins
}

void test_Result_For_Good_Estimate()
{
Expand All @@ -37,29 +43,22 @@ class GetEiTest : public CxxTest::TestSuite

void do_test_on_result_values(double input_ei, bool fixei)
{
Workspace2D_sptr testWS = createTestWorkspaceWithMonitors();
MatrixWorkspace_sptr testWS = createTestWorkspaceWithMonitors();
// This algorithm needs a name attached to the workspace
const std::string outputName("eitest");
AnalysisDataService::Instance().add(outputName, testWS);

GetEi2 alg;
alg.initialize();
alg.setPropertyValue("InputWorkspace", outputName);
alg.setProperty("Monitor1Spec", 1);
alg.setProperty("Monitor2Spec", 2);
alg.setProperty("EnergyEstimate", input_ei);
alg.setProperty("FixEi", fixei);
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg.execute());
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 = 6496.00571578;
const int expected_mon_index = 0;
const double expected_ei = (fixei) ? input_ei : 15.00322845;
const double ei = alg.getProperty("IncidentEnergy");
const double first_mon_peak = alg.getProperty("FirstMonitorPeak");
const int mon_index = alg.getProperty("FirstMonitorIndex");
const double ei = alg->getProperty("IncidentEnergy");
const double first_mon_peak = alg->getProperty("FirstMonitorPeak");
const int mon_index = alg->getProperty("FirstMonitorIndex");

TS_ASSERT_DELTA(ei, expected_ei, 1e-08);
TS_ASSERT_DELTA(first_mon_peak, expected_mon_peak, 1e-08);
Expand All @@ -69,49 +68,40 @@ class GetEiTest : public CxxTest::TestSuite
PropertyWithValue<double> *ei_propvalue = dynamic_cast<PropertyWithValue<double> *>(ei_runprop);
TS_ASSERT_DELTA((*ei_propvalue)(), expected_ei, 1e-08);

const Mantid::Kernel::Property *tzeroProp = alg.getProperty("Tzero");
const Mantid::Kernel::Property *tzeroProp = alg->getProperty("Tzero");
if(fixei)
{
TS_ASSERT(tzeroProp->isDefault());
}
else
{
// T0 value
const double tzero = alg.getProperty("Tzero");
const double tzero = alg->getProperty("Tzero");
const double expected_tzero = 3.2641273;
TS_ASSERT_DELTA(tzero, expected_tzero, 1e-08);
}

AnalysisDataService::Instance().remove(outputName);
}


void testParametersOnWorkspace()
{
Mantid::DataObjects::Workspace2D_sptr testWS = createTestWorkspaceWithMonitors();
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);

GetEi2 alg1;
alg1.initialize();
alg1.setPropertyValue("InputWorkspace", outputName);
//alg1.setProperty("Monitor1Spec", 1);
alg1.setProperty("Monitor2Spec", 2);
alg1.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg1.execute());
IAlgorithm_sptr alg;
TS_ASSERT_THROWS_NOTHING(alg = runGetEiUsingTestMonitors(outputName, 15, false));

// Test output answers
const double expected_ei = 15.00322845;
const double ei = alg1.getProperty("IncidentEnergy");
const double ei = alg->getProperty("IncidentEnergy");

TS_ASSERT_DELTA(ei, expected_ei, 1e-08);
// and verify it has been store on the run object
Expand All @@ -120,7 +110,7 @@ class GetEiTest : public CxxTest::TestSuite
TS_ASSERT_DELTA((*ei_propvalue)(), expected_ei, 1e-08);

// T0 value
const double tzero = alg1.getProperty("Tzero");
const double tzero = alg->getProperty("Tzero");
const double expected_tzero = 3.2641273;
TS_ASSERT_DELTA(tzero, expected_tzero, 1e-08);

Expand All @@ -129,69 +119,110 @@ class GetEiTest : public CxxTest::TestSuite

void testThrowsMon1()
{
Workspace2D_sptr testWS = createTestWorkspaceWithMonitors();
MatrixWorkspace_sptr testWS = createTestWorkspaceWithMonitors();
// This algorithm needs a name attached to the workspace
const std::string outputName("eitest1");
AnalysisDataService::Instance().add(outputName, testWS);

GetEi2 alg;
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::invalid_argument &e, std::string(e.what()), "Could not determine spectrum number to use. Try to set it explicitly" );
IAlgorithm_sptr alg = AlgorithmManager::Instance().createUnmanaged("GetEi", 2);
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::invalid_argument &e, std::string(e.what()), "Could not determine spectrum number to use. Try to set it explicitly" );
AnalysisDataService::Instance().remove(outputName);
}
void testThrowsEi()
{
Workspace2D_sptr testWS = createTestWorkspaceWithMonitors();
MatrixWorkspace_sptr testWS = createTestWorkspaceWithMonitors();
// This algorithm needs a name attached to the workspace
const std::string outputName("eitest2");
AnalysisDataService::Instance().add(outputName, testWS);

GetEi2 alg;
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::invalid_argument &e, std::string(e.what()), "Could not find an energy guess" );
IAlgorithm_sptr alg = AlgorithmManager::Instance().createUnmanaged("GetEi", 2);
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::invalid_argument &e,
std::string(e.what()), "Could not find an energy guess" );
AnalysisDataService::Instance().remove(outputName);
}

void testCNCS()
void test_throws_error_when_ei_not_fixed_and_no_peaks_found()
{
Mantid::DataHandling::LoadNexusMonitors ld;
std::string outws_name = "cncs";
ld.initialize();
ld.setPropertyValue("Filename","CNCS_7860_event.nxs");
ld.setPropertyValue("OutputWorkspace", outws_name);

ld.execute();
TS_ASSERT( ld.isExecuted() );

GetEi2 alg;
alg.initialize();
alg.setPropertyValue("InputWorkspace", outws_name);
TS_ASSERT_THROWS_NOTHING(alg.execute());

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", 2);
alg->initialize();
alg->setPropertyValue("InputWorkspace", outputName);
alg->setProperty("Monitor1Spec", 1);
alg->setProperty("Monitor2Spec", 2);
alg->setProperty("FixEi", false);
alg->setProperty("EnergyEstimate", 15.0);
alg->setRethrows(true);
TS_ASSERT_THROWS(alg->execute(), std::invalid_argument);

AnalysisDataService::Instance().remove(outputName);
}

// T0 value
const double tzero = alg.getProperty("Tzero");
const double expected_tzero = 61.7708;
TS_ASSERT_DELTA(tzero, expected_tzero, 1e-04);
void test_peak_time_is_zero_when_ei_fixed_and_no_peak_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;
TS_ASSERT_THROWS_NOTHING(alg = runGetEiUsingTestMonitors(outputName, 15.0, true));
if(alg)
{
const double firstMonPeak=alg->getProperty("FirstMonitorPeak");
TS_ASSERT_DELTA(firstMonPeak, 0.0, 1e-12);
}

AnalysisDataService::Instance().remove(outputName);
}

AnalysisDataService::Instance().remove(outws_name);
void testCNCS()
{
IAlgorithm_sptr ld = AlgorithmManager::Instance().createUnmanaged("LoadNexusMonitors");
std::string outws_name = "cncs";
ld->initialize();
ld->setPropertyValue("Filename","CNCS_7860_event.nxs");
ld->setPropertyValue("OutputWorkspace", outws_name);

ld->execute();
TS_ASSERT( ld->isExecuted() );

IAlgorithm_sptr alg = AlgorithmManager::Instance().createUnmanaged("GetEi", 2);
alg->initialize();
alg->setPropertyValue("InputWorkspace", outws_name);
TS_ASSERT_THROWS_NOTHING(alg->execute());


// T0 value
const double tzero = alg->getProperty("Tzero");
const double expected_tzero = 61.7708;
TS_ASSERT_DELTA(tzero, expected_tzero, 1e-04);

AnalysisDataService::Instance().remove(outws_name);
}

private:

Workspace2D_sptr createTestWorkspaceWithMonitors()
MatrixWorkspace_sptr createTestWorkspaceWithMonitors(const bool includePeaks=true)
{
const int numHists(2);
const int numBins(2000);
Workspace2D_sptr testWS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(numHists, numBins, true);
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);
Expand All @@ -206,7 +237,7 @@ class GetEiTest : public CxxTest::TestSuite
for( int i = 0; i <= numBins; ++i)
{
const double xValue = 5.0 + 5.5*i;
if( i < numBins )
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);
Expand All @@ -219,6 +250,20 @@ class GetEiTest : public CxxTest::TestSuite
return testWS;
}

IAlgorithm_sptr runGetEiUsingTestMonitors(const std::string & inputWS, const double energyGuess, const bool fixei)
{
IAlgorithm_sptr alg = AlgorithmManager::Instance().createUnmanaged("GetEi", 2);
alg->initialize();
alg->setPropertyValue("InputWorkspace", inputWS);
alg->setProperty("Monitor1Spec", 1);
alg->setProperty("Monitor2Spec", 2);
alg->setProperty("FixEi", fixei);
alg->setProperty("EnergyEstimate", energyGuess);
alg->setRethrows(true);
alg->execute();
return alg;
}

};

#endif /*GETE_ITEST_H_*/
#endif /*GETEITEST_H_*/

0 comments on commit 864ebaa

Please sign in to comment.