Skip to content

Commit

Permalink
refs #5626. Uses function against wavelengths on input.
Browse files Browse the repository at this point in the history
  • Loading branch information
@OwenArnold
OwenArnold committed Jul 19, 2012
1 parent 8254afc commit 1954025
Show file tree
Hide file tree
Showing 3 changed files with 180 additions and 77 deletions.
3 changes: 3 additions & 0 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/NormaliseByDetector.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@

#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
#include "MantidGeometry/Instrument/FitParameter.h"
#include "MantidGeometry/Instrument/Detector.h"

namespace Mantid
{
Expand Down Expand Up @@ -44,6 +46,7 @@ namespace Algorithms
virtual const std::string category() const;

private:
const Mantid::Geometry::FitParameter tryParseFunctionParameter(Mantid::Geometry::Parameter_sptr parameter, Geometry::IDetector_const_sptr det);
virtual void initDocs();
void init();
void exec();
Expand Down
199 changes: 137 additions & 62 deletions Code/Mantid/Framework/Algorithms/src/NormaliseByDetector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,94 +4,169 @@ Normalise a workspace by the detector efficiency.

#include "MantidAlgorithms/NormaliseByDetector.h"
#include "MantidKernel/System.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/Component.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/FitParameter.h"
#include "MantidGeometry/muParser_Silent.h"

using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;

namespace Mantid
{
namespace Algorithms
{
namespace Algorithms
{

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(NormaliseByDetector)

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(NormaliseByDetector)


//----------------------------------------------------------------------------------------------
/** Constructor
*/
NormaliseByDetector::NormaliseByDetector()
{
}

//----------------------------------------------------------------------------------------------
/** Destructor
*/
NormaliseByDetector::~NormaliseByDetector()
{
}


//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string NormaliseByDetector::name() const { return "NormaliseByDetector";};

/// Algorithm's version for identification. @see Algorithm::version
int NormaliseByDetector::version() const { return 1;};

/// Algorithm's category for identification. @see Algorithm::category
const std::string NormaliseByDetector::category() const { return "General";}

//----------------------------------------------------------------------------------------------
/// Sets documentation strings for this algorithm
void NormaliseByDetector::initDocs()
{
this->setWikiSummary("Normalise the input workspace by the detector efficiency.");
this->setOptionalMessage("Normalise the input workspace by the detector efficiency.");
}

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void NormaliseByDetector::init()
{
declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace","",Direction::Input), "An input workspace.");
declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output), "An output workspace.");
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void NormaliseByDetector::exec()
{
MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");

const Geometry::ParameterMap& paramMap = inWS->instrumentParameters();
//----------------------------------------------------------------------------------------------
/** Constructor
*/
NormaliseByDetector::NormaliseByDetector()
{
}

//----------------------------------------------------------------------------------------------
/** Destructor
*/
NormaliseByDetector::~NormaliseByDetector()
{
}


//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string NormaliseByDetector::name() const { return "NormaliseByDetector";};

/// Algorithm's version for identification. @see Algorithm::version
int NormaliseByDetector::version() const { return 1;};

/// Algorithm's category for identification. @see Algorithm::category
const std::string NormaliseByDetector::category() const { return "General";}

//----------------------------------------------------------------------------------------------
/// Sets documentation strings for this algorithm
void NormaliseByDetector::initDocs()
{
this->setWikiSummary("Normalise the input workspace by the detector efficiency.");
this->setOptionalMessage("Normalise the input workspace by the detector efficiency.");
}

for(size_t wsIndex = 0; wsIndex < inWS->getNumberHistograms(); ++wsIndex)
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void NormaliseByDetector::init()
{
Geometry::IDetector_const_sptr det = inWS->getDetector( wsIndex );
const std::string type = "fitting";
Geometry::Parameter_sptr fittingParam = paramMap.getRecursiveByType(&(*det), type);
if(fittingParam == NULL)
auto compositeValidator = boost::make_shared<CompositeValidator>();
compositeValidator->add(boost::make_shared<API::WorkspaceUnitValidator>("Wavelength"));
compositeValidator->add(boost::make_shared<API::HistogramValidator>());

declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace","",Direction::Input, compositeValidator),
"An input workspace in wavelength");

declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output), "An output workspace.");
}

const Geometry::FitParameter NormaliseByDetector::tryParseFunctionParameter(Geometry::Parameter_sptr parameter, Geometry::IDetector_const_sptr det)
{
if(parameter == NULL)
{
std::stringstream stream;
stream << det->getName() << " and all of it's parent components, have no fitting type parameters. This algorithm cannot be run without fitting parameters.";
this->g_log.warning(stream.str());
throw std::invalid_argument(stream.str());
}
return parameter->value<Geometry::FitParameter>();
}

setProperty("OutputWorkspace", inWS); //HACK
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void NormaliseByDetector::exec()
{
MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");

const Geometry::ParameterMap& paramMap = inWS->instrumentParameters();

IAlgorithm_sptr cloneAlg = this->createSubAlgorithm("CloneWorkspace", 0.0, 1.0, true);
cloneAlg->setRethrows(true);
cloneAlg->setProperty("InputWorkspace", inWS);
cloneAlg->setPropertyValue("OutputWorkspace", "temp");
cloneAlg->executeAsSubAlg();
Workspace_sptr temp = cloneAlg->getProperty("OutputWorkspace");
MatrixWorkspace_sptr outputWS = boost::dynamic_pointer_cast<MatrixWorkspace>(temp);

for(size_t wsIndex = 0; wsIndex < inWS->getNumberHistograms(); ++wsIndex)
{
Geometry::IDetector_const_sptr det = inWS->getDetector( wsIndex );
const std::string type = "fitting";
Geometry::Parameter_sptr foundParam = paramMap.getRecursiveByType(&(*det), type);

const Geometry::FitParameter& foundFittingParam = tryParseFunctionParameter(foundParam, det);

std::string fitFunctionName = foundFittingParam.getFunction();
IFunction_sptr function = FunctionFactory::Instance().createFunction(fitFunctionName);
typedef std::vector<std::string> ParamNames;
ParamNames allParamNames = function->getParameterNames();

// Lookup each parameter name.
for(ParamNames::iterator it = allParamNames.begin(); it != allParamNames.end(); ++it)
{
Geometry::Parameter_sptr param = paramMap.getRecursive(&(*det), (*it), type);

const Geometry::FitParameter& fitParam = tryParseFunctionParameter(param, det);

if ( fitParam.getFormula().compare("") == 0 )
{
throw std::invalid_argument("A Forumla has not been provided for a fit function");
}
else
{
std::string resultUnitStr = fitParam.getResultUnit();
if ( !resultUnitStr.empty() && resultUnitStr.compare("Wavelength") != 0)
{
throw std::invalid_argument("Units for function parameters must be in Wavelength");
}
}
mu::Parser p;
p.SetExpr(fitParam.getFormula());
double paramValue = p.Eval();
//Set the function coeffiecents.
function->setParameter(fitParam.getName(), paramValue);
}

auto wavelengths = inWS->readX(wsIndex);
const size_t nInputBins = wavelengths.size() -1;
std::vector<double> centerPointWavelength(nInputBins);
std::vector<double> outIntensity(nInputBins);
for(size_t binIndex = 0; binIndex < nInputBins; ++binIndex)
{
centerPointWavelength[binIndex] = 0.5*(wavelengths[binIndex] + wavelengths[binIndex+1]);
}
FunctionDomain1DVector domain(centerPointWavelength);
FunctionValues values(domain);
function->function(domain,values);
for(size_t i = 0; i < domain.size(); ++i)
{
outIntensity[i] = values[i];
}
outputWS->dataY(wsIndex) = outIntensity;
}

setProperty("OutputWorkspace", outputWS);
}

} // namespace Mantid
} // namespace Mantid
} // namespace Algorithms
55 changes: 40 additions & 15 deletions Code/Mantid/Framework/Algorithms/test/NormaliseByDetectorTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@
#include <cxxtest/TestSuite.h>
#include "MantidKernel/Timer.h"
#include "MantidKernel/System.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidAPI/FrameworkManager.h"
#include "MantidDataHandling/LoadEmptyInstrument.h"
#include "MantidDataHandling/LoadParameterFile.h"
#include <iostream>
Expand All @@ -22,18 +24,36 @@ class NormaliseByDetectorTest : public CxxTest::TestSuite

private:

MatrixWorkspace_sptr create_empty_instrument_with_no_fitting_functions()
MatrixWorkspace_sptr create_workspace_with_no_fitting_functions()
{
const std::string instrumentWorkspace = "InstrumentWorkspace";
LoadEmptyInstrument loadAlg;
loadAlg.initialize();
loadAlg.setPropertyValue("Filename", "POLREF_Definition.xml");
loadAlg.setPropertyValue("OutputWorkspace", instrumentWorkspace);
loadAlg.execute();
return AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(instrumentWorkspace);
IAlgorithm* loadalg = FrameworkManager::Instance().createAlgorithm("Load");
loadalg->setRethrows(true);
loadalg->initialize();
loadalg->setPropertyValue("Filename", "POLREF00004699.nxs");
loadalg->setPropertyValue("OutputWorkspace", "testws");
loadalg->execute();

// Convert units to wavelength
IAlgorithm* unitsalg = FrameworkManager::Instance().createAlgorithm("ConvertUnits");
unitsalg->initialize();
unitsalg->setPropertyValue("InputWorkspace", "testws");
unitsalg->setPropertyValue("OutputWorkspace", "testws");
unitsalg->setPropertyValue("Target", "Wavelength");
unitsalg->execute();

// Convert the specturm axis ot signed_theta
IAlgorithm* specaxisalg = FrameworkManager::Instance().createAlgorithm("ConvertSpectrumAxis");
specaxisalg->initialize();
specaxisalg->setPropertyValue("InputWorkspace", "testws");
specaxisalg->setPropertyValue("OutputWorkspace", "testws");
specaxisalg->setPropertyValue("Target", "signed_theta");
specaxisalg->execute();

WorkspaceGroup_sptr group = API::AnalysisDataService::Instance().retrieveWS<WorkspaceGroup>("testws");
return boost::dynamic_pointer_cast<MatrixWorkspace>(group->getItem(0));
}

MatrixWorkspace_sptr create_empty_instrument_with_fitting_functions()
MatrixWorkspace_sptr create_workspace_with_fitting_functions()
{
class FileObject
{
Expand Down Expand Up @@ -73,18 +93,18 @@ class NormaliseByDetectorTest : public CxxTest::TestSuite
"<parameter-file instrument = \"POLREF\" date = \"2012-01-31T00:00:00\">\n" +
"<component-link name=\"POLREF\">\n" +
"<parameter name=\"LinearBackground:A0\" type=\"fitting\">\n" +
" <formula eq=\"1.0\" result-unit=\"TOF\"/>\n" +
" <formula eq=\"1.0\" result-unit=\"Wavelength\"/>\n" +
" <fixed />\n" +
"</parameter>\n" +
"<parameter name=\"LinearBackground:A1\" type=\"fitting\">\n" +
" <formula eq=\"2.0\" result-unit=\"TOF\"/>\n" +
" <formula eq=\"2.0\" result-unit=\"Wavelength\"/>\n" +
" <fixed />\n" +
"</parameter>\n" +
"</component-link>\n" +
"</parameter-file>\n";

FileObject file(parameterFileContents, "POLREF_Parameters.xml");
MatrixWorkspace_sptr ws = create_empty_instrument_with_no_fitting_functions();
MatrixWorkspace_sptr ws = create_workspace_with_no_fitting_functions();

LoadParameterFile loadParameterAlg;
loadParameterAlg.setRethrows(true);
Expand All @@ -111,19 +131,24 @@ class NormaliseByDetectorTest : public CxxTest::TestSuite

void test_throws_when_no_fit_function_on_detector_tree()
{
MatrixWorkspace_sptr inputWS = create_empty_instrument_with_no_fitting_functions();
int i;
std::cin >> i;
MatrixWorkspace_sptr inputWS = create_workspace_with_no_fitting_functions();

NormaliseByDetector alg;
alg.setRethrows(true);
alg.initialize();
alg.setPropertyValue("OutputWorkspace", "out");
alg.setProperty("InputWorkspace", inputWS);
TS_ASSERT_THROWS(alg.execute(), std::invalid_argument);
alg.execute();
}

void test_applies_instrument_function_to_child_detectors()
{
MatrixWorkspace_sptr inputWS = create_empty_instrument_with_fitting_functions();
int i;
std::cin >> i;
// Fitting function applied to whole instrument, so no detectors should complain.
MatrixWorkspace_sptr inputWS = create_workspace_with_fitting_functions();
NormaliseByDetector alg;
alg.setRethrows(true);
alg.initialize();
Expand Down

0 comments on commit 1954025

Please sign in to comment.