diff --git a/Code/Mantid/Framework/API/CMakeLists.txt b/Code/Mantid/Framework/API/CMakeLists.txt
index 1a8ece42503b648102159f60ac84b77ebfd39880..8c24cfb88fecc34355c8fad07a0d006103cfc774 100644
--- a/Code/Mantid/Framework/API/CMakeLists.txt
+++ b/Code/Mantid/Framework/API/CMakeLists.txt
@@ -40,7 +40,7 @@ set ( SRC_FILES
src/FunctionDomain1D.cpp
src/FunctionDomainMD.cpp
src/FunctionFactory.cpp
- src/FunctionParameterDecorator.cpp
+ src/FunctionParameterDecorator.cpp
src/FunctionProperty.cpp
src/FunctionValues.cpp
src/GridDomain.cpp
@@ -61,6 +61,7 @@ set ( SRC_FILES
src/IMDHistoWorkspace.cpp
src/IMDIterator.cpp
src/IMDWorkspace.cpp
+ src/IPawleyFunction.cpp
src/IPeak.cpp
src/IPeakFunction.cpp
src/IPeaksWorkspace.cpp
@@ -185,7 +186,7 @@ set ( INC_FILES
inc/MantidAPI/FunctionDomain1D.h
inc/MantidAPI/FunctionDomainMD.h
inc/MantidAPI/FunctionFactory.h
- inc/MantidAPI/FunctionParameterDecorator.h
+ inc/MantidAPI/FunctionParameterDecorator.h
inc/MantidAPI/FunctionProperty.h
inc/MantidAPI/FunctionValues.h
inc/MantidAPI/GridDomain.h
@@ -218,6 +219,7 @@ set ( INC_FILES
inc/MantidAPI/IMDNode.h
inc/MantidAPI/IMDWorkspace.h
inc/MantidAPI/IMaskWorkspace.h
+ inc/MantidAPI/IPawleyFunction.h
inc/MantidAPI/IPeak.h
inc/MantidAPI/IPeakFunction.h
inc/MantidAPI/IPeaksWorkspace.h
@@ -327,8 +329,8 @@ set ( TEST_FILES
FuncMinimizerFactoryTest.h
FunctionAttributeTest.h
FunctionDomainTest.h
- FunctionParameterDecoratorTest.h
FunctionFactoryTest.h
+ FunctionParameterDecoratorTest.h
FunctionPropertyTest.h
FunctionTest.h
FunctionValuesTest.h
diff --git a/Code/Mantid/Framework/API/inc/MantidAPI/IPawleyFunction.h b/Code/Mantid/Framework/API/inc/MantidAPI/IPawleyFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..f560667902f0fcfd7acb7e3269363875c7e66791
--- /dev/null
+++ b/Code/Mantid/Framework/API/inc/MantidAPI/IPawleyFunction.h
@@ -0,0 +1,80 @@
+#ifndef MANTID_API_IPAWLEYFUNCTION_H_
+#define MANTID_API_IPAWLEYFUNCTION_H_
+
+#include "MantidAPI/DllConfig.h"
+#include "MantidAPI/FunctionParameterDecorator.h"
+#include "MantidAPI/IPeakFunction.h"
+
+namespace Mantid {
+namespace API {
+
+/** IPawleyFunction
+
+ This abstract class defines the interface of a PawleyFunction. An
+ implementation can be found in CurveFitting/PawleyFunction. This interface
+ exists so that the function can be used in modules outside CurveFitting.
+
+ @author Michael Wedel, Paul Scherrer Institut - SINQ
+ @date 11/03/2015
+
+ Copyright © 2015 PSI-NXMM
+
+ This file is part of Mantid.
+
+ 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.
+
+ 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.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ File change history is stored at: <https://github.com/mantidproject/mantid>
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+*/
+class MANTID_API_DLL IPawleyFunction : public FunctionParameterDecorator {
+public:
+ IPawleyFunction();
+ /// Virtual destructor.
+ virtual ~IPawleyFunction() {}
+
+ /// A string that names the crystal system.
+ virtual void setCrystalSystem(const std::string &crystalSystem) = 0;
+
+ /// Sets the name of the profile function used for the reflections.
+ virtual void setProfileFunction(const std::string &profileFunction) = 0;
+
+ /// Set the function parameters according to the supplied unit cell.
+ virtual void setUnitCell(const std::string &unitCellString) = 0;
+
+ /// Assign several peaks with the same fwhm/height parameters.
+ virtual void setPeaks(const std::vector<Kernel::V3D> &hkls, double fwhm,
+ double height) = 0;
+
+ /// Removes all peaks from the function.
+ virtual void clearPeaks() = 0;
+
+ /// Add a peak with the given parameters.
+ virtual void addPeak(const Kernel::V3D &hkl, double fwhm, double height) = 0;
+
+ /// Returns the number of peaks in the function
+ virtual size_t getPeakCount() const = 0;
+
+ /// Returns the profile function stored for the i-th peak.
+ virtual IPeakFunction_sptr getPeakFunction(size_t i) const = 0;
+
+ /// Returns the Miller indices stored for the i-th peak.
+ virtual Kernel::V3D getPeakHKL(size_t i) const = 0;
+};
+
+typedef boost::shared_ptr<IPawleyFunction> IPawleyFunction_sptr;
+
+} // namespace API
+} // namespace Mantid
+
+#endif /* MANTID_API_IPAWLEYFUNCTION_H_ */
diff --git a/Code/Mantid/Framework/API/src/IPawleyFunction.cpp b/Code/Mantid/Framework/API/src/IPawleyFunction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..caeae3cf971ea985f2ac294af5c778c749e345e5
--- /dev/null
+++ b/Code/Mantid/Framework/API/src/IPawleyFunction.cpp
@@ -0,0 +1,9 @@
+#include "MantidAPI/IPawleyFunction.h"
+
+namespace Mantid {
+namespace API {
+/// Default constructor
+IPawleyFunction::IPawleyFunction() : FunctionParameterDecorator() {}
+
+} // namespace API
+} // namespace Mantid
diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PawleyFunction.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PawleyFunction.h
index cb0e9fb7acae149213d0e76acb2d45c81e135a4a..a015b74a52fc59377592ae5cae5a15b5f3edd5c9 100644
--- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PawleyFunction.h
+++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/PawleyFunction.h
@@ -3,7 +3,7 @@
#include "MantidKernel/System.h"
#include "MantidAPI/CompositeFunction.h"
-#include "MantidAPI/FunctionParameterDecorator.h"
+#include "MantidAPI/IPawleyFunction.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/ParamFunction.h"
@@ -106,7 +106,7 @@ typedef boost::shared_ptr<PawleyParameterFunction> PawleyParameterFunction_sptr;
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
-class DLLExport PawleyFunction : public API::FunctionParameterDecorator {
+class DLLExport PawleyFunction : public API::IPawleyFunction {
public:
PawleyFunction();
virtual ~PawleyFunction() {}
@@ -143,6 +143,13 @@ public:
PawleyParameterFunction_sptr getPawleyParameterFunction() const;
protected:
+ void setPeakPositions(std::string centreName, double zeroShift,
+ const Geometry::UnitCell &cell) const;
+
+ size_t calculateFunctionValues(const API::IPeakFunction_sptr &peak,
+ const API::FunctionDomain1D &domain,
+ API::FunctionValues &localValues) const;
+
double getTransformedCenter(double d) const;
void init();
@@ -156,6 +163,8 @@ protected:
Kernel::Unit_sptr m_dUnit;
Kernel::Unit_sptr m_wsUnit;
+
+ int m_peakRadius;
};
typedef boost::shared_ptr<PawleyFunction> PawleyFunction_sptr;
diff --git a/Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp b/Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp
index 58f597657779b5a42010fba68d7f23a0acfd4a13..781e67b9c6325c37a1f27bebd12a40e468edcca5 100644
--- a/Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp
+++ b/Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp
@@ -267,9 +267,15 @@ DECLARE_FUNCTION(PawleyFunction)
/// Constructor
PawleyFunction::PawleyFunction()
- : FunctionParameterDecorator(), m_compositeFunction(),
- m_pawleyParameterFunction(), m_peakProfileComposite(), m_hkls(),
- m_dUnit(), m_wsUnit() {}
+ : IPawleyFunction(), m_compositeFunction(), m_pawleyParameterFunction(),
+ m_peakProfileComposite(), m_hkls(), m_dUnit(), m_wsUnit(),
+ m_peakRadius(5) {
+ int peakRadius;
+ if (Kernel::ConfigService::Instance().getValue("curvefitting.peakRadius",
+ peakRadius)) {
+ m_peakRadius = peakRadius;
+ }
+}
void PawleyFunction::setMatrixWorkspace(
boost::shared_ptr<const MatrixWorkspace> workspace, size_t wi,
@@ -350,6 +356,44 @@ double PawleyFunction::getTransformedCenter(double d) const {
return d;
}
+void PawleyFunction::setPeakPositions(std::string centreName, double zeroShift,
+ const UnitCell &cell) const {
+ for (size_t i = 0; i < m_hkls.size(); ++i) {
+ double centre = getTransformedCenter(cell.d(m_hkls[i]));
+
+ m_peakProfileComposite->getFunction(i)
+ ->setParameter(centreName, centre + zeroShift);
+ }
+}
+
+size_t PawleyFunction::calculateFunctionValues(
+ const API::IPeakFunction_sptr &peak, const API::FunctionDomain1D &domain,
+ API::FunctionValues &localValues) const {
+ size_t domainSize = domain.size();
+ const double *domainBegin = domain.getPointerAt(0);
+ const double *domainEnd = domain.getPointerAt(domainSize);
+
+ double centre = peak->centre();
+ double dx = m_peakRadius * peak->fwhm();
+
+ auto lb = std::lower_bound(domainBegin, domainEnd, centre - dx);
+ auto ub = std::upper_bound(lb, domainEnd, centre + dx);
+
+ size_t n = std::distance(lb, ub);
+
+ if (n == 0) {
+ throw std::invalid_argument("Null-domain");
+ }
+
+ FunctionDomain1DView localDomain(lb, n);
+ localValues.reset(localDomain);
+
+ peak->functionLocal(localValues.getPointerToCalculated(0),
+ localDomain.getPointerAt(0), n);
+
+ return std::distance(domainBegin, lb);
+}
+
/**
* Calculates the function values on the supplied domain
*
@@ -364,18 +408,38 @@ double PawleyFunction::getTransformedCenter(double d) const {
*/
void PawleyFunction::function(const FunctionDomain &domain,
FunctionValues &values) const {
- UnitCell cell = m_pawleyParameterFunction->getUnitCellFromParameters();
- double zeroShift = m_pawleyParameterFunction->getParameter("ZeroShift");
+ values.zeroCalculated();
+ try {
+ const FunctionDomain1D &domain1D =
+ dynamic_cast<const FunctionDomain1D &>(domain);
- for (size_t i = 0; i < m_hkls.size(); ++i) {
- double centre = getTransformedCenter(cell.d(m_hkls[i]));
+ UnitCell cell = m_pawleyParameterFunction->getUnitCellFromParameters();
+ double zeroShift = m_pawleyParameterFunction->getParameter("ZeroShift");
+ std::string centreName =
+ m_pawleyParameterFunction->getProfileFunctionCenterParameterName();
- m_peakProfileComposite->getFunction(i)->setParameter(
- m_pawleyParameterFunction->getProfileFunctionCenterParameterName(),
- centre + zeroShift);
- }
+ setPeakPositions(centreName, zeroShift, cell);
+
+ FunctionValues localValues;
+
+ for (size_t i = 0; i < m_peakProfileComposite->nFunctions(); ++i) {
+ IPeakFunction_sptr peak = boost::dynamic_pointer_cast<IPeakFunction>(
+ m_peakProfileComposite->getFunction(i));
- m_peakProfileComposite->function(domain, values);
+ try {
+ size_t offset = calculateFunctionValues(peak, domain1D, localValues);
+ values.addToCalculated(offset, localValues);
+ }
+ catch (std::invalid_argument) {
+ // do nothing
+ }
+ }
+
+ setPeakPositions(centreName, 0.0, cell);
+ }
+ catch (std::bad_cast) {
+ // do nothing
+ }
}
/// Removes all peaks from the function.
diff --git a/Code/Mantid/Framework/SINQ/CMakeLists.txt b/Code/Mantid/Framework/SINQ/CMakeLists.txt
index eca5e9872469ab76fc6665f55709fd70f607c36d..cbcac43f077f6b5a7d1da08316646ff66d9a3d5c 100644
--- a/Code/Mantid/Framework/SINQ/CMakeLists.txt
+++ b/Code/Mantid/Framework/SINQ/CMakeLists.txt
@@ -32,6 +32,7 @@ set ( SRC_FILES
src/PoldiUtilities/PoldiSourceSpectrum.cpp
src/PoldiUtilities/PoldiSpectrumConstantBackground.cpp
src/PoldiUtilities/PoldiSpectrumLinearBackground.cpp
+ src/PoldiUtilities/PoldiSpectrumPawleyFunction.cpp
src/PoldiUtilities/PoldiTimeTransformer.cpp
src/PoldiUtilities/UncertainValue.cpp
src/ProjectMD.cpp
@@ -78,6 +79,7 @@ set ( INC_FILES
inc/MantidSINQ/PoldiUtilities/PoldiSpectrumConstantBackground.h
inc/MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h
inc/MantidSINQ/PoldiUtilities/PoldiSpectrumLinearBackground.h
+ inc/MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h
inc/MantidSINQ/PoldiUtilities/PoldiTimeTransformer.h
inc/MantidSINQ/PoldiUtilities/UncertainValue.h
inc/MantidSINQ/PoldiUtilities/UncertainValueIO.h
@@ -119,6 +121,7 @@ set ( TEST_FILES
PoldiSpectrumDomainFunctionTest.h
PoldiSpectrumConstantBackgroundTest.h
PoldiSpectrumLinearBackgroundTest.h
+ PoldiSpectrumPawleyFunctionTest.h
PoldiTimeTransformerTest.h
PoldiTruncateDataTest.h
ProjectMDTest.h
diff --git a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h
index 165df2ca7c0bb74e5922db6f39b58832342b1ed5..a7217d53e063a97becbcd694f891652508413fdc 100644
--- a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h
+++ b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h
@@ -62,11 +62,28 @@ public:
virtual const std::string summary() const;
+ std::map<std::string, std::string> validateInputs();
+
boost::shared_ptr<Kernel::DblMatrix> getLocalCovarianceMatrix(
const boost::shared_ptr<const Kernel::DblMatrix> &covarianceMatrix,
size_t parameterOffset, size_t nParams) const;
protected:
+ Poldi2DFunction_sptr getFunctionIndividualPeaks(
+ std::string profileFunctionName,
+ const PoldiPeakCollection_sptr &peakCollection) const;
+
+ Poldi2DFunction_sptr
+ getFunctionPawley(std::string profileFunctionName,
+ const PoldiPeakCollection_sptr &peakCollection) const;
+
+ PoldiPeak_sptr
+ getPeakFromPeakFunction(API::IPeakFunction_sptr profileFunction,
+ const Kernel::V3D &hkl);
+
+ API::ITableWorkspace_sptr
+ getRefinedCellParameters(const API::IFunction_sptr &fitFunction) const;
+
PoldiPeakCollection_sptr
getPeakCollection(const DataObjects::TableWorkspace_sptr &peakTable) const;
PoldiPeakCollection_sptr getIntegratedPeakCollection(
diff --git a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h
index 707023899791b15a974caff5511213aba51df0da..83022d4639410796dabefdbcde90d80f0c41d974 100644
--- a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h
+++ b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h
@@ -41,7 +41,7 @@ namespace Poldi {
struct MANTID_SINQ_DLL Poldi2DHelper {
/// Default constructor
Poldi2DHelper()
- : dFractionalOffsets(), dOffsets(), domain(), factors(), deltaD(),
+ : dFractionalOffsets(), dOffsets(), domain(), values(), factors(), deltaD(),
minTOFN() {}
/// Transforms the chopper slit offsets for a given 2theta/distance pair.
@@ -71,6 +71,7 @@ struct MANTID_SINQ_DLL Poldi2DHelper {
}
domain = boost::make_shared<API::FunctionDomain1DVector>(current);
+ values.reset(*domain);
}
/// Calculates intensity factors for each point in the spectrum domain.
@@ -90,6 +91,7 @@ struct MANTID_SINQ_DLL Poldi2DHelper {
std::vector<int> dOffsets;
API::FunctionDomain1D_sptr domain;
+ API::FunctionValues values;
std::vector<double> factors;
double deltaD;
@@ -98,6 +100,30 @@ struct MANTID_SINQ_DLL Poldi2DHelper {
typedef boost::shared_ptr<Poldi2DHelper> Poldi2DHelper_sptr;
+class WrapAroundJacobian : public API::Jacobian {
+public:
+ WrapAroundJacobian(API::Jacobian &jacobian, size_t offset,
+ const std::vector<double> &factors, size_t factorOffset,
+ size_t domainSize)
+ : m_jacobian(jacobian), m_offset(offset), m_factors(factors),
+ m_factorOffset(factorOffset), m_domainSize(domainSize) {}
+
+ double get(size_t iY, size_t iP) { return m_jacobian.get(iY, iP); }
+
+ void set(size_t iY, size_t iP, double value) {
+ size_t realY = (m_offset + iY) % m_domainSize;
+ m_jacobian.set(realY, iP, m_jacobian.get(realY, iP) +
+ value * m_factors[m_factorOffset + iY]);
+ }
+
+protected:
+ API::Jacobian &m_jacobian;
+ size_t m_offset;
+ const std::vector<double> &m_factors;
+ size_t m_factorOffset;
+ size_t m_domainSize;
+};
+
/**
* @brief The LocalJacobian class
*
@@ -239,7 +265,6 @@ protected:
const PoldiInstrumentAdapter_sptr &poldiInstrument);
void beforeDecoratedFunctionSet(const API::IFunction_sptr &fn);
- void setProfileFunction(const std::string &profileFunctionName);
std::vector<double>
getChopperSlitOffsets(const PoldiAbstractChopper_sptr &chopper);
diff --git a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..97d973c68b3f7b03bcf341923cdaf5803e83aad9
--- /dev/null
+++ b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h
@@ -0,0 +1,65 @@
+#ifndef MANTID_SINQ_POLDISPECTRUMPAWLEYFUNCTION_H_
+#define MANTID_SINQ_POLDISPECTRUMPAWLEYFUNCTION_H_
+
+#include "MantidSINQ/DllConfig.h"
+#include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
+#include "MantidAPI/IPawleyFunction.h"
+
+namespace Mantid {
+namespace Poldi {
+
+/** PoldiSpectrumPawleyFunction : TODO: DESCRIPTION
+
+ Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
+ National Laboratory & European Spallation Source
+
+ This file is part of Mantid.
+
+ 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.
+
+ 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.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ File change history is stored at: <https://github.com/mantidproject/mantid>
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+*/
+class MANTID_SINQ_DLL PoldiSpectrumPawleyFunction
+ : public PoldiSpectrumDomainFunction {
+public:
+ PoldiSpectrumPawleyFunction();
+ virtual ~PoldiSpectrumPawleyFunction() {}
+
+ std::string name() const { return "PoldiSpectrumPawleyFunction"; }
+
+ void
+ setMatrixWorkspace(boost::shared_ptr<const API::MatrixWorkspace> workspace,
+ size_t wi, double startX, double endX);
+
+ void function1DSpectrum(const API::FunctionDomain1DSpectrum &domain,
+ API::FunctionValues &values) const;
+ void functionDeriv1DSpectrum(const API::FunctionDomain1DSpectrum &domain,
+ API::Jacobian &jacobian);
+ void poldiFunction1D(const std::vector<int> &indices,
+ const API::FunctionDomain1D &domain,
+ API::FunctionValues &values) const;
+
+ API::IPawleyFunction_sptr getPawleyFunction() const;
+
+protected:
+ void beforeDecoratedFunctionSet(const API::IFunction_sptr &fn);
+
+ API::IPawleyFunction_sptr m_pawleyFunction;
+};
+
+} // namespace Poldi
+} // namespace Mantid
+
+#endif /* MANTID_SINQ_POLDISPECTRUMPAWLEYFUNCTION_H_ */
diff --git a/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp b/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp
index ce552fe18c6fa3c9f00dc5b70fd85ea99ce4f156..d1bd62548f5d7a926540a957fb8d66ecfe9df3ea 100644
--- a/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp
+++ b/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp
@@ -1,17 +1,14 @@
-/*WIKI*
-TODO: Enter a full wiki-markup description of your algorithm here. You can then
-use the Build/wiki_maker.py script to generate your full wiki page.
-*WIKI*/
-
#include "MantidSINQ/PoldiFitPeaks2D.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/TableWorkspace.h"
+#include "MantidKernel/ListValidator.h"
#include "MantidAPI/TableRow.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/MultiDomainFunction.h"
#include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
#include "MantidSINQ/PoldiUtilities/PoldiSpectrumLinearBackground.h"
+#include "MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidSINQ/PoldiUtilities/IPoldiFunction1D.h"
@@ -19,6 +16,7 @@ use the Build/wiki_maker.py script to generate your full wiki page.
#include "MantidSINQ/PoldiUtilities/PoldiInstrumentAdapter.h"
#include "MantidSINQ/PoldiUtilities/PoldiDeadWireDecorator.h"
#include "MantidAPI/IPeakFunction.h"
+#include "MantidAPI/IPawleyFunction.h"
#include "MantidSINQ/PoldiUtilities/Poldi2DFunction.h"
@@ -57,6 +55,30 @@ const std::string PoldiFitPeaks2D::summary() const {
return "Calculates a POLDI 2D-spectrum.";
}
+/// Validate inputs for algorithm in case PawleyFit is used.
+std::map<std::string, std::string> PoldiFitPeaks2D::validateInputs() {
+ std::map<std::string, std::string> errorMap;
+
+ bool isPawleyFit = getProperty("PawleyFit");
+
+ if (isPawleyFit) {
+ Property *cellProperty = getPointerToProperty("InitialCell");
+ if (cellProperty->isDefault()) {
+ errorMap["InitialCell"] = "Initial cell must be given for PawleyFit.";
+ }
+
+ Property *refinedCellParameters =
+ getPointerToProperty("RefinedCellParameters");
+ if (refinedCellParameters->isDefault()) {
+ errorMap["RefinedCellParameters"] = "Workspace name for refined cell "
+ "parameters must be supplied for "
+ "PawleyFit.";
+ }
+ }
+
+ return errorMap;
+}
+
/// Initialization of algorithm properties.
void PoldiFitPeaks2D::init() {
declareProperty(new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "",
@@ -65,10 +87,38 @@ void PoldiFitPeaks2D::init() {
declareProperty(new WorkspaceProperty<TableWorkspace>("PoldiPeakWorkspace",
"", Direction::Input),
"Table workspace with peak information.");
- declareProperty("PeakProfileFunction", "Gaussian",
+
+ auto peakFunctionValidator = boost::make_shared<StringListValidator>(
+ FunctionFactory::Instance().getFunctionNames<IPeakFunction>());
+ declareProperty("PeakProfileFunction", "Gaussian", peakFunctionValidator,
"Profile function to use for integrating the peak profiles "
"before calculating the spectrum.");
+ declareProperty("PawleyFit", false, "Instead of refining individual peaks, "
+ "refine a unit cell. Peaks must be "
+ "indexed, an initial cell must be "
+ "provided, as well as a crystal system.");
+ declareProperty("InitialCell", "", "Initial unit cell parameters as 6 "
+ "space-separated numbers. Only used when "
+ "PawleyFit is checked.");
+
+ std::vector<std::string> crystalSystems;
+ crystalSystems.push_back("Cubic");
+ crystalSystems.push_back("Tetragonal");
+ crystalSystems.push_back("Hexagonal");
+ crystalSystems.push_back("Trigonal");
+ crystalSystems.push_back("Orthorhombic");
+ crystalSystems.push_back("Monoclinic");
+ crystalSystems.push_back("Triclinic");
+
+ auto crystalSystemValidator =
+ boost::make_shared<StringListValidator>(crystalSystems);
+
+ declareProperty("CrystalSystem", "Cubic", crystalSystemValidator,
+ "Crystal system to use for constraining "
+ "unit cell parameters. Only used when "
+ "PawleyFit is checked.");
+
declareProperty("FitConstantBackground", true,
"Add a constant background term to the fit.");
declareProperty("ConstantBackgroundParameter", 0.0,
@@ -98,6 +148,91 @@ void PoldiFitPeaks2D::init() {
declareProperty(new WorkspaceProperty<TableWorkspace>(
"RefinedPoldiPeakWorkspace", "", Direction::Output),
"Table workspace with fitted peaks.");
+
+ declareProperty(new WorkspaceProperty<ITableWorkspace>(
+ "RefinedCellParameters", "", Direction::Output, PropertyMode::Optional));
+}
+
+/// Creates a PoldiPeak from the given profile function/hkl pair.
+PoldiPeak_sptr
+PoldiFitPeaks2D::getPeakFromPeakFunction(IPeakFunction_sptr profileFunction,
+ const V3D &hkl) {
+
+ IAlgorithm_sptr errorAlg = createChildAlgorithm("EstimatePeakErrors");
+ errorAlg->setProperty(
+ "Function", boost::dynamic_pointer_cast<IFunction>(profileFunction));
+ errorAlg->setPropertyValue("OutputWorkspace", "Errors");
+ errorAlg->execute();
+
+ double centre = profileFunction->centre();
+ double height = profileFunction->height();
+ double fwhmValue = profileFunction->fwhm();
+
+ ITableWorkspace_sptr errorTable = errorAlg->getProperty("OutputWorkspace");
+
+ double centreError = errorTable->cell<double>(0, 2);
+ double heightError = errorTable->cell<double>(1, 2);
+ double fwhmError = errorTable->cell<double>(2, 2);
+
+ UncertainValue d(centre, centreError);
+ UncertainValue intensity(height, heightError);
+ UncertainValue fwhm(fwhmValue, fwhmError);
+
+ PoldiPeak_sptr peak =
+ PoldiPeak::create(MillerIndices(hkl), d, intensity, UncertainValue(1.0));
+ peak->setFwhm(fwhm, PoldiPeak::FwhmRelation::AbsoluteD);
+
+ return peak;
+}
+
+/// Returns a TableWorkspace with refined cell parameters and error.
+ITableWorkspace_sptr PoldiFitPeaks2D::getRefinedCellParameters(
+ const IFunction_sptr &fitFunction) const {
+ Poldi2DFunction_sptr poldi2DFunction =
+ boost::dynamic_pointer_cast<Poldi2DFunction>(fitFunction);
+
+ if (!poldi2DFunction || poldi2DFunction->nFunctions() < 1) {
+ throw std::invalid_argument(
+ "Cannot process function that is not a Poldi2DFunction.");
+ }
+
+ ITableWorkspace_sptr latticeParameterTable =
+ WorkspaceFactory::Instance().createTable();
+
+ latticeParameterTable->addColumn("str", "Parameter");
+ latticeParameterTable->addColumn("double", "Value");
+ latticeParameterTable->addColumn("double", "Error");
+
+ // The first function should be PoldiSpectrumPawleyFunction
+ boost::shared_ptr<PoldiSpectrumPawleyFunction> poldiPawleyFunction =
+ boost::dynamic_pointer_cast<PoldiSpectrumPawleyFunction>(
+ poldi2DFunction->getFunction(0));
+
+ if (!poldiPawleyFunction) {
+ throw std::invalid_argument("First function in Poldi2DFunction is not "
+ "PoldiSpectrumPawleyFunction.");
+ }
+
+ IPawleyFunction_sptr pawleyFunction =
+ boost::dynamic_pointer_cast<IPawleyFunction>(
+ poldiPawleyFunction->getDecoratedFunction());
+
+ if (pawleyFunction) {
+ CompositeFunction_sptr pawleyParts =
+ boost::dynamic_pointer_cast<CompositeFunction>(
+ pawleyFunction->getDecoratedFunction());
+
+ IFunction_sptr pawleyParameters = pawleyParts->getFunction(0);
+
+ for (size_t i = 0; i < pawleyParameters->nParams(); ++i) {
+ TableRow newRow = latticeParameterTable->appendRow();
+ newRow << pawleyParameters->parameterName(i)
+ << pawleyParameters->getParameter(i)
+ << pawleyParameters->getError(i);
+ }
+ }
+
+ return latticeParameterTable;
}
/**
@@ -134,7 +269,8 @@ boost::shared_ptr<DblMatrix> PoldiFitPeaks2D::getLocalCovarianceMatrix(
/**
* Construct a PoldiPeakCollection from a Poldi2DFunction
*
- * This method performs the opposite operation of getFunctionFromPeakCollection.
+ * This method performs the opposite operation of
+ *getFunctionFromPeakCollection.
* It takes a function, checks if it's of the proper type and turns the
*information
* into a PoldiPeakCollection.
@@ -161,6 +297,44 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
size_t offset = 0;
for (size_t i = 0; i < poldi2DFunction->nFunctions(); ++i) {
+ boost::shared_ptr<PoldiSpectrumPawleyFunction> poldiPawleyFunction =
+ boost::dynamic_pointer_cast<PoldiSpectrumPawleyFunction>(
+ poldi2DFunction->getFunction(i));
+
+ if (poldiPawleyFunction) {
+ IPawleyFunction_sptr pawleyFunction =
+ poldiPawleyFunction->getPawleyFunction();
+
+ if (pawleyFunction) {
+ CompositeFunction_sptr decoratedFunction =
+ boost::dynamic_pointer_cast<CompositeFunction>(
+ pawleyFunction->getDecoratedFunction());
+
+ offset = decoratedFunction->getFunction(0)->nParams();
+
+ for (size_t j = 0; j < pawleyFunction->getPeakCount(); ++j) {
+ IPeakFunction_sptr profileFunction =
+ pawleyFunction->getPeakFunction(j);
+
+ size_t nLocalParams = profileFunction->nParams();
+ boost::shared_ptr<Kernel::DblMatrix> localCov =
+ getLocalCovarianceMatrix(covarianceMatrix, offset, nLocalParams);
+ profileFunction->setCovarianceMatrix(localCov);
+
+ // Increment offset for next function
+ offset += nLocalParams;
+
+ V3D peakHKL = pawleyFunction->getPeakHKL(j);
+
+ PoldiPeak_sptr peak =
+ getPeakFromPeakFunction(profileFunction, peakHKL);
+
+ normalizedPeaks->addPeak(peak);
+ }
+ }
+ break;
+ }
+
boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction =
boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>(
poldi2DFunction->getFunction(i));
@@ -179,30 +353,8 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
// Increment offset for next function
offset += nLocalParams;
- IAlgorithm_sptr errorAlg = createChildAlgorithm("EstimatePeakErrors");
- errorAlg->setProperty(
- "Function", boost::dynamic_pointer_cast<IFunction>(profileFunction));
- errorAlg->setPropertyValue("OutputWorkspace", "Errors");
- errorAlg->execute();
-
- double centre = profileFunction->centre();
- double height = profileFunction->height();
- double fwhmValue = profileFunction->fwhm();
-
- ITableWorkspace_sptr errorTable =
- errorAlg->getProperty("OutputWorkspace");
-
- double centreError = errorTable->cell<double>(0, 2);
- double heightError = errorTable->cell<double>(1, 2);
- double fwhmError = errorTable->cell<double>(2, 2);
-
- UncertainValue d(centre, centreError);
- UncertainValue intensity(height, heightError);
- UncertainValue fwhm(fwhmValue, fwhmError);
-
PoldiPeak_sptr peak =
- PoldiPeak::create(MillerIndices(), d, intensity, UncertainValue(1.0));
- peak->setFwhm(fwhm, PoldiPeak::FwhmRelation::AbsoluteD);
+ getPeakFromPeakFunction(profileFunction, V3D(0, 0, 0));
normalizedPeaks->addPeak(peak);
}
@@ -212,25 +364,25 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction(
}
/**
- * Constructs a proper function from a peak collection
+ * Returns a Poldi2DFunction that encapsulates individual peaks
*
- * This method constructs a Poldi2DFunction and assigns one
- *PoldiSpectrumDomainFunction to it for
- * each peak contained in the peak collection.
+ * This function takes all peaks from the supplied peak collection and
+ *generates
+ * an IPeakFunction of the type given in the name parameter, wraps them
+ * in a Poldi2DFunction and returns it.
*
- * @param peakCollection :: PoldiPeakCollection containing peaks with integral
- *intensities
- * @return Poldi2DFunction with one PoldiSpectrumDomainFunction per peak
+ * @param profileFunctionName :: Profile function name.
+ * @param peakCollection :: Peak collection with peaks to be used in the fit.
+ * @return :: A Poldi2DFunction with peak profile functions.
*/
-Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionFromPeakCollection(
+Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionIndividualPeaks(
+ std::string profileFunctionName,
const PoldiPeakCollection_sptr &peakCollection) const {
Poldi2DFunction_sptr mdFunction(new Poldi2DFunction);
for (size_t i = 0; i < peakCollection->peakCount(); ++i) {
PoldiPeak_sptr peak = peakCollection->peak(i);
- std::string profileFunctionName = getProperty("PeakProfileFunction");
-
boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction =
boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>(
FunctionFactory::Instance().createFunction(
@@ -259,6 +411,91 @@ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionFromPeakCollection(
return mdFunction;
}
+/**
+ * Returns a Poldi2DFunction that encapsulates a PawleyFunction
+ *
+ * This function creates a PawleyFunction using the supplied profile function
+ * name and the crystal system as well as initial cell from the input
+ *properties
+ * of the algorithm and wraps it in a Poldi2DFunction.
+ *
+ * Because the peak intensities are integral at this step but PawleyFunction
+ * expects peak heights, a profile function is created and
+ *setIntensity/height-
+ * methods are used to convert.
+ *
+ * @param profileFunctionName :: Profile function name for PawleyFunction.
+ * @param peakCollection :: Peak collection with peaks to be used in the fit.
+ * @return :: A Poldi2DFunction with a PawleyFunction.
+ */
+Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionPawley(
+ std::string profileFunctionName,
+ const PoldiPeakCollection_sptr &peakCollection) const {
+ Poldi2DFunction_sptr mdFunction(new Poldi2DFunction);
+
+ boost::shared_ptr<PoldiSpectrumPawleyFunction> poldiPawleyFunction =
+ boost::dynamic_pointer_cast<PoldiSpectrumPawleyFunction>(
+ FunctionFactory::Instance().createFunction(
+ "PoldiSpectrumPawleyFunction"));
+
+ if (!poldiPawleyFunction) {
+ throw std::invalid_argument("Could not create pawley function.");
+ }
+
+ poldiPawleyFunction->setDecoratedFunction("PawleyFunction");
+
+ IPawleyFunction_sptr pawleyFunction =
+ poldiPawleyFunction->getPawleyFunction();
+ pawleyFunction->setProfileFunction(profileFunctionName);
+
+ std::string crystalSystem = getProperty("CrystalSystem");
+ pawleyFunction->setCrystalSystem(crystalSystem);
+
+ std::string initialCell = getProperty("InitialCell");
+ pawleyFunction->setUnitCell(initialCell);
+
+ IPeakFunction_sptr pFun = boost::dynamic_pointer_cast<IPeakFunction>(
+ FunctionFactory::Instance().createFunction(profileFunctionName));
+
+ for (size_t i = 0; i < peakCollection->peakCount(); ++i) {
+ PoldiPeak_sptr peak = peakCollection->peak(i);
+
+ pFun->setCentre(peak->d());
+ pFun->setFwhm(peak->fwhm(PoldiPeak::AbsoluteD));
+ pFun->setIntensity(peak->intensity());
+
+ pawleyFunction->addPeak(peak->hkl().asV3D(),
+ peak->fwhm(PoldiPeak::AbsoluteD), pFun->height());
+ }
+
+ pawleyFunction->fix(pawleyFunction->parameterIndex("f0.ZeroShift"));
+ mdFunction->addFunction(poldiPawleyFunction);
+
+ return mdFunction;
+}
+
+/**
+ * Constructs a proper function from a peak collection
+ *
+ * This method constructs a Poldi2DFunction depending on whether or not a
+ * Pawley fit is performed each peak contained in the peak collection.
+ *
+ * @param peakCollection :: PoldiPeakCollection containing peaks with integral
+ *intensities
+ * @return Poldi2DFunction with one PoldiSpectrumDomainFunction per peak
+ */
+Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionFromPeakCollection(
+ const PoldiPeakCollection_sptr &peakCollection) const {
+ std::string profileFunctionName = getProperty("PeakProfileFunction");
+
+ bool pawleyFit = getProperty("PawleyFit");
+ if (pawleyFit) {
+ return getFunctionPawley(profileFunctionName, peakCollection);
+ }
+
+ return getFunctionIndividualPeaks(profileFunctionName, peakCollection);
+}
+
/// Executes the algorithm
void PoldiFitPeaks2D::exec() {
TableWorkspace_sptr peakTable = getProperty("PoldiPeakWorkspace");
@@ -296,6 +533,18 @@ void PoldiFitPeaks2D::exec() {
setProperty("OutputWorkspace", getWorkspace(fitAlgorithm));
setProperty("RefinedPoldiPeakWorkspace", integralPeaks->asTableWorkspace());
setProperty("Calculated1DSpectrum", outWs1D);
+
+ bool isPawleyFit = getProperty("PawleyFit");
+ if (isPawleyFit) {
+ ITableWorkspace_sptr cellParameters = getRefinedCellParameters(fitFunction);
+
+ if (cellParameters->rowCount() > 0) {
+ setProperty("RefinedCellParameters",
+ getRefinedCellParameters(fitFunction));
+ } else {
+ g_log.warning() << "Warning: Cell parameter table is empty.";
+ }
+ }
}
/**
@@ -370,7 +619,11 @@ IAlgorithm_sptr PoldiFitPeaks2D::calculateSpectrum(
fit->setProperty("Minimizer", "Levenberg-MarquardtMD");
+ // Setting the level to Notice to avoid problems with Levenberg-MarquardtMD.
+ int oldLogLevel = g_log.getLevel();
+ g_log.setLevelForAll(Poco::Message::PRIO_NOTICE);
fit->execute();
+ g_log.setLevelForAll(oldLogLevel);
return fit;
}
@@ -497,14 +750,17 @@ void PoldiFitPeaks2D::setTimeTransformer(
/**
* Extracts time bin width from workspace parameter
*
- * The method uses the difference between first and second x-value of the first
+ * The method uses the difference between first and second x-value of the
+ *first
*spectrum as
* time bin width. If the workspace does not contain proper data (0 spectra or
*less than
- * 2 x-values), the method throws an std::invalid_argument-exception. Otherwise
+ * 2 x-values), the method throws an std::invalid_argument-exception.
+ *Otherwise
*it calls setDeltaT.
*
- * @param matrixWorkspace :: MatrixWorkspace with at least one spectrum with at
+ * @param matrixWorkspace :: MatrixWorkspace with at least one spectrum with
+ *at
*least two x-values.
*/
void PoldiFitPeaks2D::setDeltaTFromWorkspace(
@@ -520,12 +776,14 @@ void PoldiFitPeaks2D::setDeltaTFromWorkspace(
"Cannot process MatrixWorkspace with less than 2 x-values.");
}
- // difference between first and second x-value is assumed to be the bin width.
+ // difference between first and second x-value is assumed to be the bin
+ // width.
setDeltaT(matrixWorkspace->readX(0)[1] - matrixWorkspace->readX(0)[0]);
}
/**
- * Assigns delta t, throws std::invalid_argument on invalid value (determined by
+ * Assigns delta t, throws std::invalid_argument on invalid value (determined
+ *by
*isValidDeltaT).
*
* @param newDeltaT :: Value to be used as delta t for calculations.
@@ -567,7 +825,8 @@ PoldiFitPeaks2D::getPeakCollection(const TableWorkspace_sptr &peakTable) const {
/**
* Return peak collection with integrated peaks
*
- * This method takes a PoldiPeakCollection where the intensity is represented by
+ * This method takes a PoldiPeakCollection where the intensity is represented
+ *by
*the maximum. Then
* it takes the profile function stored in the peak collection, which must be
*the name of a registered
@@ -607,7 +866,8 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getIntegratedPeakCollection(
}
/* If no profile function is specified, it's not possible to get integrated
- * intensities at all and we try to use the one specified by the user instead.
+ * intensities at all and we try to use the one specified by the user
+ * instead.
*/
std::string profileFunctionName = rawPeakCollection->getProfileFunctionName();
@@ -688,7 +948,8 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getNormalizedPeakCollection(
/**
* Converts normalized peak intensities to count based integral intensities
*
- * This operation is the opposite of getNormalizedPeakCollection and is used to
+ * This operation is the opposite of getNormalizedPeakCollection and is used
+ *to
*convert
* the intensities back to integral intensities.
*
diff --git a/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumDomainFunction.cpp b/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumDomainFunction.cpp
index 46d186b5309063320907c7ca3f967035459c3302..1eb8ca7a8a2ade535938ea5e17afe2ea7d3ef953 100644
--- a/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumDomainFunction.cpp
+++ b/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumDomainFunction.cpp
@@ -116,7 +116,7 @@ void PoldiSpectrumDomainFunction::functionDeriv1DSpectrum(
Poldi2DHelper_sptr helper = m_2dHelpers[index];
if (helper) {
- int domainSize = static_cast<int>(domain.size());
+ size_t domainSize = domain.size();
double fwhm = m_profileFunction->fwhm();
double centre = m_profileFunction->centre();
@@ -135,29 +135,18 @@ void PoldiSpectrumDomainFunction::functionDeriv1DSpectrum(
}
}
- size_t np = m_profileFunction->nParams();
-
size_t baseOffset = static_cast<size_t>(pos + helper->minTOFN);
for (size_t i = 0; i < helper->dOffsets.size(); ++i) {
- LocalJacobian smallJ(dWidthN, np);
-
double newD = centre + helper->dFractionalOffsets[i];
size_t offset = static_cast<size_t>(helper->dOffsets[i]) + baseOffset;
+ WrapAroundJacobian smallJ(jacobian, offset, helper->factors, pos,
+ domainSize);
m_profileFunction->setCentre(newD);
m_profileFunction->functionDerivLocal(
&smallJ, helper->domain->getPointerAt(pos), dWidthN);
-
- for (size_t j = 0; j < dWidthN; ++j) {
- size_t off = (offset + j) % domainSize;
- for (size_t p = 0; p < np; ++p) {
- jacobian.set(off, p,
- jacobian.get(off, p) +
- smallJ.getRaw(j, p) * helper->factors[pos + j]);
- }
- }
}
m_profileFunction->setCentre(centre);
diff --git a/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumPawleyFunction.cpp b/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumPawleyFunction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ddab49664864c7fafd2a70072690847caeb4e121
--- /dev/null
+++ b/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumPawleyFunction.cpp
@@ -0,0 +1,104 @@
+#include "MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h"
+#include "MantidAPI/FunctionFactory.h"
+#include "MantidAPI/FunctionValues.h"
+
+namespace Mantid {
+namespace Poldi {
+
+using namespace API;
+
+/// Default constructor
+PoldiSpectrumPawleyFunction::PoldiSpectrumPawleyFunction()
+ : PoldiSpectrumDomainFunction(), m_pawleyFunction() {}
+
+/// This function does nothing to prevent setting the workspace on the wrapped
+/// function (unit conversion will not work and is not needed).
+void PoldiSpectrumPawleyFunction::setMatrixWorkspace(
+ boost::shared_ptr<const MatrixWorkspace> workspace, size_t wi,
+ double startX, double endX) {
+ UNUSED_ARG(workspace);
+ UNUSED_ARG(wi);
+ UNUSED_ARG(startX);
+ UNUSED_ARG(endX);
+}
+
+/// Implementation of function1DSpectrum that transforms PawleyFunction output.
+void PoldiSpectrumPawleyFunction::function1DSpectrum(
+ const FunctionDomain1DSpectrum &domain, FunctionValues &values) const {
+ values.zeroCalculated();
+
+ size_t domainSize = domain.size();
+ size_t index = domain.getWorkspaceIndex();
+ Poldi2DHelper_sptr helper = m_2dHelpers[index];
+
+ if (helper) {
+ for (size_t i = 0; i < helper->dOffsets.size(); ++i) {
+ double newDOffset =
+ helper->dOffsets[i] * helper->deltaD + helper->dFractionalOffsets[i];
+ m_pawleyFunction->setParameter("f0.ZeroShift", newDOffset);
+
+ size_t baseOffset = helper->minTOFN;
+
+ m_pawleyFunction->function(*(helper->domain), helper->values);
+
+ for (size_t j = 0; j < helper->values.size(); ++j) {
+ values.addToCalculated((j + baseOffset) % domainSize,
+ helper->values[j] * helper->factors[j]);
+ }
+ }
+
+ m_pawleyFunction->setParameter("f0.ZeroShift", 0.0);
+ }
+}
+
+/// Using numerical derivatives turned out to be faster for this case.
+void PoldiSpectrumPawleyFunction::functionDeriv1DSpectrum(
+ const FunctionDomain1DSpectrum &domain, Jacobian &jacobian) {
+ calNumericalDeriv(domain, jacobian);
+}
+
+/// Calculate 1D function by adding the functions calculated for each index.
+void
+PoldiSpectrumPawleyFunction::poldiFunction1D(const std::vector<int> &indices,
+ const FunctionDomain1D &domain,
+ FunctionValues &values) const {
+ FunctionValues localValues(domain);
+
+ m_pawleyFunction->function(domain, localValues);
+
+ double chopperSlitCount = static_cast<double>(m_chopperSlitOffsets.size());
+
+ for (auto index = indices.begin(); index != indices.end(); ++index) {
+ std::vector<double> factors(domain.size());
+
+ for (size_t i = 0; i < factors.size(); ++i) {
+ values.addToCalculated(i,
+ chopperSlitCount * localValues[i] *
+ m_timeTransformer->detectorElementIntensity(
+ domain[i], static_cast<size_t>(*index)));
+ }
+ }
+}
+
+/// Returns the internally stored Pawley function.
+IPawleyFunction_sptr PoldiSpectrumPawleyFunction::getPawleyFunction() const {
+ return m_pawleyFunction;
+}
+
+/// Makes sure that the decorated function is of the right type.
+void PoldiSpectrumPawleyFunction::beforeDecoratedFunctionSet(
+ const IFunction_sptr &fn) {
+ IPawleyFunction_sptr pawleyFunction =
+ boost::dynamic_pointer_cast<IPawleyFunction>(fn);
+
+ if (!pawleyFunction) {
+ throw std::invalid_argument("Function is not a pawley function.");
+ }
+
+ m_pawleyFunction = pawleyFunction;
+}
+
+DECLARE_FUNCTION(PoldiSpectrumPawleyFunction)
+
+} // namespace Poldi
+} // namespace Mantid
diff --git a/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h b/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h
index 1ac7a69edaaf3814ff1e7001da18d7d75adc1249..88f8e3f7f6a6d8091fa3ca2809cd59489089b625 100644
--- a/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h
+++ b/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h
@@ -147,10 +147,6 @@ public:
PoldiPeakCollection_sptr noProfilePeaks(new PoldiPeakCollection);
TS_ASSERT_THROWS_NOTHING(spectrumCalculator.getIntegratedPeakCollection(noProfilePeaks));
- // While setting an invalid function name throws.
- spectrumCalculator.setProperty("PeakProfileFunction", "InvalidFunctionName");
- TS_ASSERT_THROWS(spectrumCalculator.getIntegratedPeakCollection(noProfilePeaks), std::runtime_error);
-
// When there is no valid PoldiPeakCollection, the method also throws
PoldiPeakCollection_sptr invalidPeakCollection;
TS_ASSERT_THROWS(spectrumCalculator.getIntegratedPeakCollection(invalidPeakCollection), std::invalid_argument);
diff --git a/Code/Mantid/Framework/SINQ/test/PoldiSpectrumPawleyFunctionTest.h b/Code/Mantid/Framework/SINQ/test/PoldiSpectrumPawleyFunctionTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..6adb145afa8521a2c8b2a71776b6165b696e8a4d
--- /dev/null
+++ b/Code/Mantid/Framework/SINQ/test/PoldiSpectrumPawleyFunctionTest.h
@@ -0,0 +1,191 @@
+#ifndef MANTID_SINQ_POLDISPECTRUMPAWLEYFUNCTIONTEST_H_
+#define MANTID_SINQ_POLDISPECTRUMPAWLEYFUNCTIONTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidSINQ/PoldiUtilities/PoldiSpectrumDomainFunction.h"
+#include "MantidSINQ/PoldiUtilities/PoldiMockInstrumentHelpers.h"
+#include "MantidSINQ/PoldiUtilities/PoldiInstrumentAdapter.h"
+#include "MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h"
+#include "MantidKernel/V3D.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include "MantidAPI/FunctionFactory.h"
+
+#include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
+using namespace Mantid::Poldi;
+using namespace Mantid::API;
+using namespace Mantid::Kernel;
+
+using ::testing::_;
+using ::testing::Mock;
+
+class MockPawleyFunction : public IPawleyFunction {
+public:
+ MockPawleyFunction() {}
+ ~MockPawleyFunction() {}
+ // IFunction interface
+ MOCK_CONST_METHOD0(name, std::string());
+ MOCK_CONST_METHOD2(function, void(const FunctionDomain &, FunctionValues &));
+
+ // IPawleyFunction interface
+ MOCK_METHOD1(setCrystalSystem, void(const std::string &));
+ MOCK_METHOD1(setProfileFunction, void(const std::string &));
+ MOCK_METHOD1(setUnitCell, void(const std::string &));
+
+ MOCK_METHOD3(setPeaks, void(const std::vector<V3D> &, double, double));
+ MOCK_METHOD0(clearPeaks, void());
+ MOCK_METHOD3(addPeak, void(const V3D &, double, double));
+ MOCK_CONST_METHOD0(getPeakCount, size_t(void));
+ MOCK_CONST_METHOD1(getPeakFunction, IPeakFunction_sptr(size_t));
+ MOCK_CONST_METHOD1(getPeakHKL, V3D(size_t));
+
+ MOCK_METHOD4(setMatrixWorkspace,
+ void(MatrixWorkspace_const_sptr, size_t, double, double));
+
+protected:
+ void init() { setDecoratedFunction("Gaussian"); }
+};
+
+DECLARE_FUNCTION(MockPawleyFunction)
+
+class PoldiSpectrumPawleyFunctionTest : 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 PoldiSpectrumPawleyFunctionTest *createSuite() {
+ return new PoldiSpectrumPawleyFunctionTest();
+ }
+ static void destroySuite(PoldiSpectrumPawleyFunctionTest *suite) {
+ delete suite;
+ }
+
+ PoldiSpectrumPawleyFunctionTest() {
+ m_detector = boost::shared_ptr<ConfiguredHeliumDetector>(
+ new ConfiguredHeliumDetector);
+ m_chopper = boost::shared_ptr<MockChopper>(new MockChopper);
+
+ m_spectrum = PoldiSourceSpectrum_sptr(new ConfiguredSpectrum);
+
+ EXPECT_CALL(*m_chopper, distanceFromSample())
+ .WillRepeatedly(Return(11800.0));
+
+ EXPECT_CALL(*m_chopper, zeroOffset()).WillRepeatedly(Return(0.15));
+
+ m_instrument = PoldiInstrumentAdapter_sptr(new FakePoldiInstrumentAdapter);
+ }
+
+ void test_setDecoratedFunction() {
+ PoldiSpectrumPawleyFunction fn;
+ fn.initialize();
+
+ TS_ASSERT_THROWS_NOTHING(fn.setDecoratedFunction("PawleyFunction"));
+ TS_ASSERT_THROWS(fn.setDecoratedFunction("Gaussian"),
+ std::invalid_argument);
+ TS_ASSERT_THROWS(fn.setDecoratedFunction("CompositeFunction"),
+ std::invalid_argument);
+ }
+
+ void test_getPawleyFunction() {
+ PoldiSpectrumPawleyFunction fn;
+ fn.initialize();
+
+ TS_ASSERT(!fn.getPawleyFunction());
+ fn.setDecoratedFunction("PawleyFunction");
+
+ IFunction_sptr pawleyFn = fn.getDecoratedFunction();
+ TS_ASSERT(pawleyFn);
+ TS_ASSERT(boost::dynamic_pointer_cast<IPawleyFunction>(pawleyFn));
+
+ TS_ASSERT(fn.getPawleyFunction());
+ }
+
+ void test_setMatrixWorkspace() {
+ PoldiSpectrumPawleyFunction fn;
+ fn.initialize();
+ fn.setDecoratedFunction("MockPawleyFunction");
+
+ MatrixWorkspace_const_sptr ws =
+ WorkspaceCreationHelper::Create2DWorkspace123(4, 10);
+
+ // Make sure the setMatrixWorkspace method can be called directly.
+ IPawleyFunction_sptr pFn = fn.getPawleyFunction();
+ boost::shared_ptr<MockPawleyFunction> mpFn =
+ boost::dynamic_pointer_cast<MockPawleyFunction>(pFn);
+ EXPECT_CALL(*mpFn, setMatrixWorkspace(_, _, _, _)).Times(1);
+
+ mpFn->setMatrixWorkspace(ws, 0, 0.0, 0.0);
+
+ /* Make sure the decorated function does not get the matrix workspace
+ * so that there are no unit problems (poldi workspaces are in time,
+ * calculation needs to be done in d-spacing).
+ */
+ fn.setMatrixWorkspace(ws, 0, 0.0, 0.0);
+
+ TS_ASSERT(Mock::VerifyAndClearExpectations(&(*mpFn)));
+ }
+
+ void test_function1DSpectrum() {
+ TestablePoldiSpectrumPawleyFunction fn;
+ fn.setDecoratedFunction("PawleyFunction");
+
+ IPawleyFunction_sptr pFn = fn.getPawleyFunction();
+ pFn->setProfileFunction("Gaussian");
+ pFn->setCrystalSystem("Cubic");
+ // Only the first figure matters, because of cubic
+ pFn->setUnitCell("5.43122617238802162554 5.431 5.431 90 90 90");
+ pFn->addPeak(V3D(4, 2, 2), 0.0027446316797104233, 679.59369981039407842726);
+
+ fn.m_deltaT = 3.0;
+ fn.initializeInstrumentParameters(m_instrument);
+
+ std::vector<double> xvalues(500, 1.0);
+
+ FunctionDomain1DSpectrum domain(342, xvalues);
+ TS_ASSERT_EQUALS(domain.getWorkspaceIndex(), 342);
+ FunctionValues values(domain);
+ values.setCalculated(0.0);
+
+ fn.function(domain, values);
+
+ std::vector<double> reference;
+ reference.push_back(0.214381692355321);
+ reference.push_back(1.4396533098854);
+ reference.push_back(7.69011673999647);
+ reference.push_back(32.6747845396612);
+ reference.push_back(110.432605589092);
+ reference.push_back(296.883931458002);
+ reference.push_back(634.864220660384);
+ reference.push_back(1079.89069118744);
+ reference.push_back(1461.11207069126);
+ reference.push_back(1572.50503614829);
+ reference.push_back(1346.18685763306);
+ reference.push_back(916.691981263516);
+ reference.push_back(496.502218342172);
+ reference.push_back(213.861997764049);
+ reference.push_back(73.2741206547921);
+ reference.push_back(19.9697293956518);
+ reference.push_back(4.32910692237627);
+ reference.push_back(0.746498624291666);
+ reference.push_back(0.102391587633906);
+
+ for(size_t i = 0; i < reference.size(); ++i) {
+ TS_ASSERT_DELTA(values[479 + i] / reference[i], 1.0, 1e-12);
+ }
+ }
+
+private:
+ class TestablePoldiSpectrumPawleyFunction
+ : public PoldiSpectrumPawleyFunction {
+ friend class PoldiSpectrumPawleyFunctionTest;
+ };
+
+ boost::shared_ptr<ConfiguredHeliumDetector> m_detector;
+ boost::shared_ptr<MockChopper> m_chopper;
+ PoldiSourceSpectrum_sptr m_spectrum;
+
+ PoldiInstrumentAdapter_sptr m_instrument;
+};
+
+#endif /* MANTID_SINQ_POLDISPECTRUMPAWLEYFUNCTIONTEST_H_ */
diff --git a/Code/Mantid/Testing/SystemTests/tests/analysis/POLDIFitPeaks2DTest.py b/Code/Mantid/Testing/SystemTests/tests/analysis/POLDIFitPeaks2DTest.py
index 99a0154e0b54f4b1548827b7d4cee21d68e53ca5..2accbd059d2f94895a423017eab7692483bec224 100644
--- a/Code/Mantid/Testing/SystemTests/tests/analysis/POLDIFitPeaks2DTest.py
+++ b/Code/Mantid/Testing/SystemTests/tests/analysis/POLDIFitPeaks2DTest.py
@@ -1,11 +1,11 @@
-#pylint: disable=no-init,invalid-name
+#pylint: disable=no-init,invalid-name,too-many-locals
import stresstesting
from mantid.simpleapi import *
import numpy as np
-'''The system test currently checks that the calculation of 2D spectra
-works correctly.'''
class POLDIFitPeaks2DTest(stresstesting.MantidStressTest):
+ """The system test currently checks that the calculation of 2D spectra
+works correctly."""
def runTest(self):
dataFiles = ["poldi2013n006904"]
@@ -81,3 +81,27 @@ class POLDIFitPeaks2DTest(stresstesting.MantidStressTest):
self.assertTrue(np.all(xDataCalc == xDataRef))
self.assertLessThan(maxDifference, 0.07)
+
+class POLDIFitPeaks2DPawleyTest(stresstesting.MantidStressTest):
+ def runTest(self):
+ si = PoldiLoadRuns(2013, 6903, 6904, 2)
+ corr = PoldiAutoCorrelation('si_data_6904')
+ peaks = PoldiPeakSearch(corr)
+ peaks_ref, fit_plots = PoldiFitPeaks1D(corr, PoldiPeakTable='peaks')
+ si_refs = PoldiCreatePeaksFromCell("F d -3 m", "Si 0 0 0", a=5.431, LatticeSpacingMin=0.7)
+ indexed = PoldiIndexKnownCompounds(peaks_ref, "si_refs")
+
+ DeleteTableRows("indexed_si_refs", "10-30")
+
+ fit2d, fit1d, peaks_ref_2d, cell = PoldiFitPeaks2D('si_data_6904', 'indexed_si_refs',
+ PawleyFit=True,
+ InitialCell="5.431 5.431 5.431 90 90 90",
+ CrystalSystem="Cubic",
+ MaximumIterations=100)
+
+ cell_a = cell.cell(0, 1)
+ cell_a_err = cell.cell(0, 2)
+
+ self.assertLessThan(np.abs(cell_a_err), 5.0e-5)
+ self.assertLessThan(np.abs(cell_a - 5.4311946) / cell_a_err, 1.5)
+
diff --git a/Code/Mantid/docs/source/algorithms/PoldiFitPeaks2D-v1.rst b/Code/Mantid/docs/source/algorithms/PoldiFitPeaks2D-v1.rst
index 2e2b532ad0eedd6db92125293d3123dce0412773..374e32c6fd36c56d482851862222cf5cbe5a237d 100644
--- a/Code/Mantid/docs/source/algorithms/PoldiFitPeaks2D-v1.rst
+++ b/Code/Mantid/docs/source/algorithms/PoldiFitPeaks2D-v1.rst
@@ -13,7 +13,7 @@ PoldiFitPeaks2D is an algorithm that can be used to fit a set of individual peak
The 1D-peak intensities need to be integral intensities, so the peaks are integrated if necessary. If there is no profile information supplied in the peak table (:ref:`algm-PoldiFitPeaks1D` adds this automatically), it's possible to supply a profile function as parameter to this algorithm. If a profile function name is present in the peak table, the one supplied in the parameters has priority.
-At the moment all profiles are calculated independently, using Gaussian functions. In future versions of the algorithm this will be much more flexible.
+There are two modes for performing the fit. In the default mode, all peak profiles are fitted independently using the same function that is used for integration. The other possibility is to perform a Pawley-type fit, where peak positions are calculated using lattice parameters and Miller indices (see :ref:`algm-PawleyFit` for a more general explanation of the method). This mode is controlled by the PawleyFit parameter, if it is activated, InitialCell and CrystalSystem need to be specified as well. For these fits, an additional table will be created which contains the refined lattice parameters. Please note that the peaks need to be indexed to use this mode (:ref:`algm-PoldiCreatePeaksFromCell`, :ref:`algm-PoldiIndexKnownCompounds`).
PoldiFitPeaks2D can also be used to calculate a theoretical 2D pattern from a set of peaks by limiting the iterations to 0.
@@ -24,19 +24,17 @@ Usage
.. include:: ../usagedata-note.txt
+**Individual peak profiles**
+
PoldiFitPeaks2D operates on a MatrixWorkspace with a valid POLDI instrument definition. The following short example demonstrates how to use the algorithm, processing data obtained from recording the spectrum of a Silicon standard material (powder) and calculating a theoretical 2D-spectrum.
.. testcode:: ExSilicon2D
# Load data file with Si spectrum and instrument definition
- raw_6904 = LoadSINQFile(Filename = "poldi2013n006904.hdf", Instrument = "POLDI")
- LoadInstrument(raw_6904, InstrumentName = "POLDI")
-
- # Data needs to be truncated to correct dimensions.
- truncated_6904 = PoldiTruncateData(raw_6904)
-
+ truncated = PoldiLoadRuns(2013, 6904)
+
# Perform correlation, peak search and fit
- correlated_6904 = PoldiAutoCorrelation(truncated_6904)
+ correlated_6904 = PoldiAutoCorrelation("truncated_data_6904")
peaks_6904 = PoldiPeakSearch(correlated_6904)
PoldiFitPeaks1D(InputWorkspace = correlated_6904, FwhmMultiples = 4.0,
@@ -45,7 +43,7 @@ PoldiFitPeaks2D operates on a MatrixWorkspace with a valid POLDI instrument defi
FitPlotsWorkspace = "fit_plots_6904")
# Calculate a 2D spectrum using the refined peaks
- PoldiFitPeaks2D(InputWorkspace=truncated_6904,
+ PoldiFitPeaks2D(InputWorkspace="truncated_data_6904",
PoldiPeakWorkspace="peaks_refined_6904",
RefinedPoldiPeakWorkspace="peaks_fit_2d_6904",
Calculated1DSpectrum="simulated_1d_6904",
@@ -65,14 +63,10 @@ In general, there is a background in POLDI data that depends on :math:`2\theta`.
.. testcode:: ExSilicon2DBackground
# Load data file with Si spectrum and instrument definition
- raw_6904 = LoadSINQFile(Filename = "poldi2013n006904.hdf", Instrument = "POLDI")
- LoadInstrument(raw_6904, InstrumentName = "POLDI")
-
- # Data needs to be truncated to correct dimensions.
- truncated_6904 = PoldiTruncateData(raw_6904)
+ truncated = PoldiLoadRuns(2013, 6904)
# Perform correlation, peak search and fit
- correlated_6904 = PoldiAutoCorrelation(truncated_6904)
+ correlated_6904 = PoldiAutoCorrelation("truncated_data_6904")
peaks_6904 = PoldiPeakSearch(correlated_6904)
PoldiFitPeaks1D(InputWorkspace = correlated_6904, FwhmMultiples = 4.0,
@@ -81,7 +75,7 @@ In general, there is a background in POLDI data that depends on :math:`2\theta`.
FitPlotsWorkspace = "fit_plots_6904")
# Calculate a 2D spectrum using the refined peaks - with background linear in 2theta
- PoldiFitPeaks2D(InputWorkspace=truncated_6904,
+ PoldiFitPeaks2D(InputWorkspace="truncated_data_6904",
PoldiPeakWorkspace="peaks_refined_6904",
OutputWorkspace="simulated_6904",
RefinedPoldiPeakWorkspace="peaks_fit_2d_6904",
@@ -106,4 +100,62 @@ Furthermore, a 1D diffractogram is also calculated, which shows all peaks that w
Calculated diffractogram for Silicon powder standard.
+**Pawley-type fit**
+
+The following example shows an example for refinement of lattice parameters using the PawleyFit-option.
+
+.. testcode:: ExSilicon2DPawley
+
+ import numpy as np
+
+ # Load and merge 2 data files for better statistics.
+ truncated = PoldiLoadRuns(2013, 6903, 6904, 2)
+
+ # Perform correlation, peak search and fit
+ correlated_6904 = PoldiAutoCorrelation("truncated_data_6904")
+ peaks_6904 = PoldiPeakSearch(correlated_6904)
+
+ PoldiFitPeaks1D(InputWorkspace = correlated_6904, FwhmMultiples = 4.0,
+ PeakFunction = "Gaussian", PoldiPeakTable = peaks_6904,
+ OutputWorkspace = "peaks_refined_6904",
+ FitPlotsWorkspace = "fit_plots_6904")
+
+ # Generate reflections for Silicon
+ si_peaks = PoldiCreatePeaksFromCell(SpaceGroup = "F d -3 m",
+ Atoms = "Si 0.0 0.0 0.0",
+ a = 5.431,
+ LatticeSpacingMin = 0.7)
+ # Index the refined peaks
+ indexed = PoldiIndexKnownCompounds("peaks_refined_6904",
+ CompoundWorkspaces = "si_peaks")
+
+ # Only consider the first 8 peaks
+ DeleteTableRows("indexed_si_peaks", "8-30")
+
+ # Fit a unit cell.
+ PoldiFitPeaks2D(InputWorkspace="truncated_data_6904",
+ PoldiPeakWorkspace="indexed_si_peaks",
+ OutputWorkspace="fitted_6904",
+ PawleyFit = True,
+ InitialCell = "5.431 5.431 5.431 90 90 90",
+ CrystalSystem = "Cubic",
+ MaximumIterations=100,
+ RefinedPoldiPeakWorkspace="peaks_fit_2d_6904",
+ Calculated1DSpectrum="simulated_1d_6904",
+ RefinedCellParameters="refined_cell_6904")
+
+
+ lattice_parameters = AnalysisDataService.retrieve("refined_cell_6904")
+
+ cell_a = np.round(lattice_parameters.cell(0, 1), 5)
+ cell_a_error = np.round(lattice_parameters.cell(0, 2), 5)
+
+ print "Refined lattice parameter a =", cell_a, "+/-", cell_a_error
+
+The refined lattice parameter is printed at the end:
+
+.. testoutput:: ExSilicon2DPawley
+
+ Refined lattice parameter a = 5.43126 +/- 5e-05
+
.. categories::