Skip to content

Commit

Permalink
refs #8948. Unverified changes. Need proper tests for this.
Browse files Browse the repository at this point in the history
  • Loading branch information
@OwenArnold
OwenArnold committed Mar 3, 2014
1 parent d9b7d28 commit 7d30703
Show file tree
Hide file tree
Showing 3 changed files with 222 additions and 45 deletions.
81 changes: 46 additions & 35 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/PolarisationCorrection.h
View file
Expand Up @@ -3,54 +3,65 @@

#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
#include <boost/shared_ptr.hpp>
#include <boost/optional.hpp>

namespace Mantid
{
namespace Algorithms
{
namespace API
{
class WorkspaceGroup;
class MatrixWorkspace;
}
namespace Algorithms
{

/** PolarisationCorrection : TODO: DESCRIPTION
Copyright &copy; 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
/** PolarisationCorrection : Algorithm to perform polarisation corrections on multi-period group workspaces.
This file is part of Mantid.
Copyright &copy; 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
Mantid is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
This file is part of Mantid.
Mantid is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
Mantid is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
Mantid is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class DLLExport PolarisationCorrection : public API::Algorithm
{
public:
PolarisationCorrection();
virtual ~PolarisationCorrection();

virtual const std::string name() const;
virtual int version() const;
virtual const std::string category() const;
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
private:
virtual void initDocs();
void init();
void exec();
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class DLLExport PolarisationCorrection: public API::Algorithm
{
public:
PolarisationCorrection();
virtual ~PolarisationCorrection();

virtual const std::string name() const;
virtual int version() const;
virtual const std::string category() const;

};
private:
virtual void initDocs();
void init();
void exec();
boost::shared_ptr<Mantid::API::MatrixWorkspace> execPolynomialCorrection(boost::shared_ptr<Mantid::API::MatrixWorkspace>& input, const std::vector<double>& coefficients);
boost::shared_ptr<Mantid::API::WorkspaceGroup> execPA(boost::shared_ptr<Mantid::API::WorkspaceGroup> inWS);
boost::shared_ptr<Mantid::API::MatrixWorkspace> multiply(boost::shared_ptr<Mantid::API::MatrixWorkspace> & lhs, boost::shared_ptr<Mantid::API::MatrixWorkspace> & rhs);
boost::shared_ptr<Mantid::API::MatrixWorkspace> multiply(boost::shared_ptr<Mantid::API::MatrixWorkspace> & lhs, const double & rhs);
boost::shared_ptr<Mantid::API::MatrixWorkspace> add(boost::shared_ptr<Mantid::API::MatrixWorkspace> & lhs, boost::shared_ptr<Mantid::API::MatrixWorkspace> & rhs);
boost::shared_ptr<Mantid::API::MatrixWorkspace> add(boost::shared_ptr<Mantid::API::MatrixWorkspace>& lhs, const double& toAdd);

};

} // namespace Algorithms
} // namespace Algorithms
} // namespace Mantid

#endif /* MANTID_ALGORITHMS_POLARISATIONCORRECTION_H_ */
#endif /* MANTID_ALGORITHMS_POLARISATIONCORRECTION_H_ */
177 changes: 167 additions & 10 deletions Code/Mantid/Framework/Algorithms/src/PolarisationCorrection.cpp
View file
Expand Up @@ -6,11 +6,16 @@

#include "MantidAlgorithms/PolarisationCorrection.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/ListValidator.h"
#include <algorithm>

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

namespace
{
Expand All @@ -25,13 +30,46 @@ namespace
return "PA";
}

const std::string crhoLabel()
{
return "crho";
}

const std::string cppLabel()
{
return "cPp";
}

const std::string cAlphaLabel()
{
return "cAlpha";
}

const std::string cApLabel()
{
return "cAp";
}

std::vector<std::string> modes()
{
std::vector<std::string> modes;
modes.push_back(pALabel());
modes.push_back(pNRLabel());
return modes;
}

Instrument_const_sptr fetchInstrument(WorkspaceGroup const * const groupWS)
{
if (groupWS->size() == 0)
{
throw std::invalid_argument("Input group workspace has no children.");
}
Workspace_sptr firstWS = groupWS->getItem(0);
MatrixWorkspace_sptr matrixWS = boost::dynamic_pointer_cast<MatrixWorkspace>(firstWS);
return matrixWS->getInstrument();
}

typedef std::vector<double> VecDouble;
}

namespace Mantid
Expand Down Expand Up @@ -90,19 +128,132 @@ namespace Mantid
*/
void PolarisationCorrection::init()
{
declareProperty(new WorkspaceProperty<WorkspaceGroup>("InputWorkspace", "", Direction::Input),
"An input workspace.");
declareProperty(
new WorkspaceProperty<Mantid::API::WorkspaceGroup>("InputWorkspace", "", Direction::Input),
"An input workspace to process.");

auto propOptions = modes();
declareProperty("PolarisationAnalysis", "PA",
boost::make_shared<StringListValidator>(propOptions),
declareProperty("PolarisationAnalysis", "PA", boost::make_shared<StringListValidator>(propOptions),
"What Polarisation mode will be used?\n"
"PNR: Polarised Neutron Reflectivity mode\n"
"PA: Full Polarisation Analysis PNR-PA");
declareProperty(new WorkspaceProperty<WorkspaceGroup>("OutputWorkspace", "", Direction::Output),
"PNR: Polarised Neutron Reflectivity mode\n"
"PA: Full Polarisation Analysis PNR-PA");

VecDouble emptyVec;
auto mandatoryArray = boost::make_shared<MandatoryValidator<VecDouble> >();

declareProperty(new ArrayProperty<double>(crhoLabel(), mandatoryArray, Direction::Input),
"TODO-Description");
declareProperty(new ArrayProperty<double>(cppLabel(), mandatoryArray, Direction::Input),
"TODO - Description");
declareProperty(new ArrayProperty<double>(cAlphaLabel(), mandatoryArray, Direction::Input),
"TODO - Description");
declareProperty(new ArrayProperty<double>(cApLabel(), mandatoryArray, Direction::Input),
"TODO - Description");

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

MatrixWorkspace_sptr PolarisationCorrection::execPolynomialCorrection(MatrixWorkspace_sptr& input,
const VecDouble& coefficients)
{
auto polyCorr = this->createChildAlgorithm("PolynomialCorrection");
polyCorr->initialize();
polyCorr->setProperty("InputWorkspace", input);
polyCorr->setProperty("Coefficients", coefficients);
polyCorr->execute();
MatrixWorkspace_sptr corrected = polyCorr->getProperty("OutputWorkspace");
return corrected;
}

MatrixWorkspace_sptr PolarisationCorrection::multiply(MatrixWorkspace_sptr& lhs,
MatrixWorkspace_sptr& rhs)
{
auto multiply = this->createChildAlgorithm("Multiply");
multiply->initialize();
multiply->setProperty("LHSWorkspace", lhs);
multiply->setProperty("RHSWorkspace", rhs);
multiply->execute();
MatrixWorkspace_sptr outWS = multiply->getProperty("OutputWorskapce");
return outWS;
}

MatrixWorkspace_sptr PolarisationCorrection::multiply(MatrixWorkspace_sptr& lhs, const double& rhs)
{
auto multiply = this->createChildAlgorithm("Multiply");
multiply->initialize();
multiply->setProperty("LHSWorkspace", lhs);
multiply->setProperty("RHSWorkspace", rhs);
multiply->execute();
MatrixWorkspace_sptr outWS = multiply->getProperty("OutputWorskapce");
return outWS;
}

MatrixWorkspace_sptr PolarisationCorrection::add(MatrixWorkspace_sptr& lhs, const double& rhs)
{
auto plus = this->createChildAlgorithm("Plus");
plus->initialize();
plus->setProperty("LHSWorkspace", lhs);
plus->setProperty("RHSWorkspace", rhs);
plus->execute();
MatrixWorkspace_sptr outWS = plus->getProperty("OutputWorskapce");
return outWS;
}

WorkspaceGroup_sptr PolarisationCorrection::execPA(WorkspaceGroup_sptr inWS)
{

size_t itemIndex = 0;
MatrixWorkspace_sptr Ipp = boost::dynamic_pointer_cast<MatrixWorkspace>(
inWS->getItem(itemIndex++));
MatrixWorkspace_sptr Iaa = boost::dynamic_pointer_cast<MatrixWorkspace>(
inWS->getItem(itemIndex++));
MatrixWorkspace_sptr Ipa = boost::dynamic_pointer_cast<MatrixWorkspace>(
inWS->getItem(itemIndex++));
MatrixWorkspace_sptr Iap = boost::dynamic_pointer_cast<MatrixWorkspace>(
inWS->getItem(itemIndex++));

MatrixWorkspace_sptr ones = WorkspaceFactory::Instance().create(Iaa);
ones = this->add(ones, 1.0);

const VecDouble c_rho = getProperty(crhoLabel());
const VecDouble c_alpha = getProperty(cAlphaLabel());
const VecDouble c_pp = getProperty(cppLabel());
const VecDouble c_ap = getProperty(cApLabel());

const auto rho = this->execPolynomialCorrection(ones, c_rho);
const auto pp = this->execPolynomialCorrection(ones, c_pp);
const auto alpha = this->execPolynomialCorrection(ones, c_alpha);
const auto ap = this->execPolynomialCorrection(ones, c_ap);

const auto A0 = Iaa * pp + ap * Ipa * rho * pp + ap * Iap * pp * alpha
+ Ipp * ap * alpha * rho * pp;
const auto A1 = pp * Iaa;
const auto A2 = pp * Iap;
const auto A3 = ap * Iaa;
const auto A4 = ap * Ipa;
const auto A5 = ap * alpha * Ipp;
const auto A6 = ap * alpha * Iap;
const auto A7 = pp * rho * Ipp;
const auto A8 = pp * rho * Ipa;

const auto D = pp * ap * (rho + alpha + 1.0 + rho * alpha);

const auto nIpp = (A0 - A1 + A2 - A3 + A4 + A5 - A6 + A7 - A8 + Ipp + Iaa - Ipa - Iap) / D;
const auto nIaa = (A0 + A1 - A2 + A3 - A4 - A5 + A6 - A7 + A8 + Ipp + Iaa - Ipa - Iap) / D;
const auto nIpa = (A0 - A1 + A2 + A3 - A4 - A5 + A6 + A7 - A8 - Ipp - Iaa + Ipa + Iap) / D;
const auto nIap = (A0 + A1 - A2 - A3 + A4 + A5 - A6 - A7 + A8 - Ipp - Iaa + Ipa + Iap) / D;

WorkspaceGroup_sptr dataOut = boost::make_shared<WorkspaceGroup>();
dataOut->addWorkspace(nIpp);
dataOut->addWorkspace(nIaa);
dataOut->addWorkspace(nIpa);
dataOut->addWorkspace(nIap);

return DataOut;
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
Expand All @@ -111,19 +262,25 @@ namespace Mantid
WorkspaceGroup_sptr inWS = getProperty("InputWorkspace");
const std::string analysisMode = getProperty("PolarisationAnalysis");
const size_t nWorkspaces = inWS->size();
if(analysisMode == pALabel())

Instrument_const_sptr instrument = fetchInstrument(inWS.get());

WorkspaceGroup_sptr outWS;
if (analysisMode == pALabel())
{
if( nWorkspaces != 4)
if (nWorkspaces != 4)
{
throw std::invalid_argument("For PA analysis, input group must have 4 periods.");
}
outWS = execPA(inWS);
}
else if(analysisMode == pNRLabel())
else if (analysisMode == pNRLabel())
{
if (nWorkspaces != 2)
{
throw std::invalid_argument("For PNR analysis, input group must have 2 periods.");
}

}
}

Expand Down
9 changes: 9 additions & 0 deletions Code/Mantid/Framework/Algorithms/test/PolarisationCorrectionTest.h
View file
Expand Up @@ -79,6 +79,11 @@ class PolarisationCorrectionTest: public CxxTest::TestSuite
alg.initialize();
alg.setProperty("InputWorkspace", inWS);
alg.setProperty("PolarisationAnalysis", "PA");
alg.setPropertyValue("crho", "1,1,1,1");
alg.setPropertyValue("calpha", "1,1,1,1");
alg.setPropertyValue("cAp", "1,1,1,1");
alg.setPropertyValue("cPp", "1,1,1,1");

alg.setPropertyValue("OutputWorkspace", outWSName);
TSM_ASSERT_THROWS("Wrong number of grouped workspaces, should throw", alg.execute(),
std::invalid_argument&);
Expand All @@ -98,6 +103,10 @@ class PolarisationCorrectionTest: public CxxTest::TestSuite
alg.setProperty("InputWorkspace", inWS);
alg.setProperty("PolarisationAnalysis", "PNR");
alg.setPropertyValue("OutputWorkspace", outWSName);
alg.setPropertyValue("crho", "1,1,1,1");
alg.setPropertyValue("calpha", "1,1,1,1");
alg.setPropertyValue("cAp", "1,1,1,1");
alg.setPropertyValue("cPp", "1,1,1,1");
TSM_ASSERT_THROWS("Wrong number of grouped workspaces, should throw", alg.execute(),
std::invalid_argument&);
}
Expand Down

0 comments on commit 7d30703

Please sign in to comment.