diff --git a/Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp b/Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp index b54db9de6e7ba28d3e7084d05b49d98d24f0f134..bfaa272e26cac04350802bbfdb9d207d433671e6 100644 --- a/Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/PawleyFunction.cpp @@ -267,9 +267,8 @@ DECLARE_FUNCTION(PawleyFunction) /// Constructor PawleyFunction::PawleyFunction() - : IPawleyFunction(), 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() {} void PawleyFunction::setMatrixWorkspace( boost::shared_ptr<const MatrixWorkspace> workspace, size_t wi, @@ -366,16 +365,24 @@ void PawleyFunction::function(const FunctionDomain &domain, FunctionValues &values) const { UnitCell cell = m_pawleyParameterFunction->getUnitCellFromParameters(); double zeroShift = m_pawleyParameterFunction->getParameter("ZeroShift"); + std::string centreName = + m_pawleyParameterFunction->getProfileFunctionCenterParameterName(); for (size_t i = 0; i < m_hkls.size(); ++i) { double centre = getTransformedCenter(cell.d(m_hkls[i])); - m_peakProfileComposite->getFunction(i)->setParameter( - m_pawleyParameterFunction->getProfileFunctionCenterParameterName(), - centre + zeroShift); + m_peakProfileComposite->getFunction(i) + ->setParameter(centreName, centre + zeroShift); } m_peakProfileComposite->function(domain, values); + + for (size_t i = 0; i < m_hkls.size(); ++i) { + m_peakProfileComposite->getFunction(i) + ->setParameter(centreName, m_peakProfileComposite->getFunction(i) + ->getParameter(centreName) - + zeroShift); + } } /// Removes all peaks from the function. diff --git a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h index dc2d381822e324b24fffd5f45d4f864a9865e578..343d2e7476482b480979cfff1660826e083029a0 100644 --- a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h +++ b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h @@ -5,6 +5,7 @@ #include "MantidSINQ/DllConfig.h" #include "MantidAPI/Algorithm.h" #include "MantidAPI/IFunction.h" +#include "MantidAPI/IPeakFunction.h" #include "MantidDataObjects/TableWorkspace.h" #include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h" #include "MantidSINQ/PoldiUtilities/PoldiTimeTransformer.h" @@ -58,6 +59,18 @@ public: virtual const std::string summary() const; + 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) const; + protected: PoldiPeakCollection_sptr getPeakCollection(const DataObjects::TableWorkspace_sptr &peakTable) const; diff --git a/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp b/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp index 735094834f0288d4e726c83a446d37d785719773..24023846edad6fd6e7e3c8b527919f1c4ae85119 100644 --- a/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp +++ b/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp @@ -12,6 +12,7 @@ use the Build/wiki_maker.py script to generate your full wiki page. #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 +20,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" @@ -50,9 +52,7 @@ const std::string PoldiFitPeaks2D::name() const { return "PoldiFitPeaks2D"; } int PoldiFitPeaks2D::version() const { return 1; } /// Algorithm's category for identification. @see Algorithm::category -const std::string PoldiFitPeaks2D::category() const { - return "SINQ\\Poldi"; -} +const std::string PoldiFitPeaks2D::category() const { return "SINQ\\Poldi"; } /// Very short algorithm summary. @see Algorith::summary const std::string PoldiFitPeaks2D::summary() const { @@ -71,6 +71,17 @@ void PoldiFitPeaks2D::init() { "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."); + declareProperty("CrystalSystem", "", "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, @@ -102,6 +113,47 @@ void PoldiFitPeaks2D::init() { "Table workspace with fitted peaks."); } +PoldiPeak_sptr +PoldiFitPeaks2D::getPeakFromPeakFunction(IPeakFunction_sptr profileFunction, + const V3D &hkl) const { + double centre = profileFunction->centre(); + double height = profileFunction->height(); + + size_t dIndex = 0; + size_t iIndex = 0; + size_t fIndex = 0; + + for (size_t j = 0; j < profileFunction->nParams(); ++j) { + if (profileFunction->getParameter(j) == centre) { + dIndex = j; + } else if (profileFunction->getParameter(j) == height) { + iIndex = j; + } else { + fIndex = j; + } + } + + // size_t dIndex = peakFunction->parameterIndex("Centre"); + UncertainValue d(profileFunction->getParameter(dIndex), + profileFunction->getError(dIndex)); + + // size_t iIndex = peakFunction->parameterIndex("Area"); + UncertainValue intensity(profileFunction->getParameter(iIndex), + profileFunction->getError(iIndex)); + + // size_t fIndex = peakFunction->parameterIndex("Sigma"); + double fwhmValue = profileFunction->fwhm(); + UncertainValue fwhm(fwhmValue, fwhmValue / + profileFunction->getParameter(fIndex) * + profileFunction->getError(fIndex)); + + PoldiPeak_sptr peak = + PoldiPeak::create(MillerIndices(hkl), d, intensity, UncertainValue(1.0)); + peak->setFwhm(fwhm, PoldiPeak::FwhmRelation::AbsoluteD); + + return peak; +} + /** * Construct a PoldiPeakCollection from a Poldi2DFunction * @@ -127,77 +179,72 @@ PoldiPeakCollection_sptr PoldiFitPeaks2D::getPeakCollectionFromFunction( boost::make_shared<PoldiPeakCollection>(PoldiPeakCollection::Integral); for (size_t i = 0; i < poldi2DFunction->nFunctions(); ++i) { - boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction = - boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>( + boost::shared_ptr<PoldiSpectrumPawleyFunction> poldiPawleyFunction = + boost::dynamic_pointer_cast<PoldiSpectrumPawleyFunction>( poldi2DFunction->getFunction(i)); - if (peakFunction) { - IPeakFunction_sptr profileFunction = - boost::dynamic_pointer_cast<IPeakFunction>( - peakFunction->getProfileFunction()); + if (poldiPawleyFunction) { + IPawleyFunction_sptr pawleyFunction = + poldiPawleyFunction->getPawleyFunction(); + if (pawleyFunction) { + CompositeFunction_sptr wrappedFn = + boost::dynamic_pointer_cast<CompositeFunction>( + pawleyFunction->getDecoratedFunction()); + IFunction_sptr pawleyParamFn = wrappedFn->getFunction(0); + for (size_t j = 0; j < pawleyParamFn->nParams(); ++j) { + std::cout << j << " " << pawleyParamFn->parameterName(j) << " " + << pawleyParamFn->getParameter(j) << " " + << pawleyParamFn->getError(j) << std::endl; + } + + for (size_t j = 0; j < pawleyFunction->getPeakCount(); ++j) { + IPeakFunction_sptr profileFunction = + pawleyFunction->getPeakFunction(j); + V3D peakHKL = pawleyFunction->getPeakHKL(j); - double centre = profileFunction->centre(); - double height = profileFunction->height(); + PoldiPeak_sptr peak = + getPeakFromPeakFunction(profileFunction, peakHKL); - size_t dIndex = 0; - size_t iIndex = 0; - size_t fIndex = 0; + normalizedPeaks->addPeak(peak); - for (size_t j = 0; j < profileFunction->nParams(); ++j) { - if (profileFunction->getParameter(j) == centre) { - dIndex = j; - } else if (profileFunction->getParameter(j) == height) { - iIndex = j; - } else { - fIndex = j; + for (size_t p = 0; p < profileFunction->nParams(); ++p) { + std::cout << j << " " << p << " " + << profileFunction->parameterName(p) << " " + << profileFunction->getParameter(p) << std::endl; + } } } + break; + } - // size_t dIndex = peakFunction->parameterIndex("Centre"); - UncertainValue d(peakFunction->getParameter(dIndex), - peakFunction->getError(dIndex)); - - // size_t iIndex = peakFunction->parameterIndex("Area"); - UncertainValue intensity(peakFunction->getParameter(iIndex), - peakFunction->getError(iIndex)); + boost::shared_ptr<PoldiSpectrumDomainFunction> peakFunction = + boost::dynamic_pointer_cast<PoldiSpectrumDomainFunction>( + poldi2DFunction->getFunction(i)); - // size_t fIndex = peakFunction->parameterIndex("Sigma"); - double fwhmValue = profileFunction->fwhm(); - UncertainValue fwhm(fwhmValue, fwhmValue / - peakFunction->getParameter(fIndex) * - peakFunction->getError(fIndex)); + if (peakFunction) { + IPeakFunction_sptr profileFunction = + boost::dynamic_pointer_cast<IPeakFunction>( + peakFunction->getProfileFunction()); 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); + continue; } } return normalizedPeaks; } -/** - * Constructs a proper function from a peak collection - * - * This method constructs a Poldi2DFunction and assigns one - *PoldiSpectrumDomainFunction to it for - * 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( +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( @@ -226,6 +273,86 @@ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionFromPeakCollection( return mdFunction; } +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()); + + IPeakFunction_sptr p = pawleyFunction->getPeakFunction(i); + V3D h = pawleyFunction->getPeakHKL(i); + std::cout << p->centre() << " " << p->fwhm() << " " << h << std::endl; + } + + pawleyFunction->fix(pawleyFunction->parameterIndex("f0.ZeroShift")); + + mdFunction->addFunction(poldiPawleyFunction); + + for (size_t i = 0; i < mdFunction->nParams(); ++i) { + std::cout << i << " " << mdFunction->parameterName(i) << " " + << mdFunction->getParameter(i) << std::endl; + } + + return mdFunction; +} + +/** + * Constructs a proper function from a peak collection + * + * This method constructs a Poldi2DFunction and assigns one + *PoldiSpectrumDomainFunction to it for + * 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"); diff --git a/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumPawleyFunction.cpp b/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumPawleyFunction.cpp index ef2d24267fff9ffc053cba2ee1c943348cadb9f5..32094fb95ae044c7e2fe711ed145c182729373f2 100644 --- a/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumPawleyFunction.cpp +++ b/Code/Mantid/Framework/SINQ/src/PoldiUtilities/PoldiSpectrumPawleyFunction.cpp @@ -34,13 +34,16 @@ void PoldiSpectrumPawleyFunction::function1DSpectrum( FunctionValues localValues(*helper->domain); for (size_t i = 0; i < helper->dOffsets.size(); ++i) { - double newDOffset = helper->dOffsets[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), localValues); for (size_t j = 0; j < localValues.size(); ++j) { - values.addToCalculated(j % domainSize, + values.addToCalculated((j + baseOffset) % domainSize, localValues[j] * helper->factors[j]); } } @@ -60,18 +63,23 @@ void PoldiSpectrumPawleyFunction::functionDeriv1DSpectrum( if (helper) { for (size_t i = 0; i < helper->dOffsets.size(); ++i) { - double newDOffset = helper->dOffsets[i]; + double newDOffset = + helper->dOffsets[i] * helper->deltaD + helper->dFractionalOffsets[i]; m_pawleyFunction->setParameter("f0.ZeroShift", newDOffset); + size_t baseOffset = helper->minTOFN; + LocalJacobian localJacobian(ny, nParams); m_pawleyFunction->functionDeriv(*(helper->domain), localJacobian); for (size_t j = 0; j < ny; ++j) { - size_t wrapped = j % domainSize; + size_t wrapped = (j + baseOffset) % domainSize; for (size_t p = 0; p < nParams; ++p) { - jacobian.set(wrapped, p, - jacobian.get(wrapped, p) + - localJacobian.getRaw(j, p) * helper->factors[j]); + if (m_pawleyFunction->isActive(p)) { + jacobian.set(wrapped, p, + jacobian.get(wrapped, p) + + localJacobian.getRaw(j, p) * helper->factors[j]); + } } } } @@ -86,19 +94,20 @@ PoldiSpectrumPawleyFunction::poldiFunction1D(const std::vector<int> &indices, FunctionValues &values) const {} IPawleyFunction_sptr PoldiSpectrumPawleyFunction::getPawleyFunction() const { - return m_pawleyFunction; + return m_pawleyFunction; } -void PoldiSpectrumPawleyFunction::beforeDecoratedFunctionSet(const IFunction_sptr &fn) -{ - IPawleyFunction_sptr pawleyFunction = - boost::dynamic_pointer_cast<IPawleyFunction>(fn); +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."); - } + if (!pawleyFunction) { + throw std::invalid_argument("Function is not a pawley function."); + } - m_pawleyFunction = pawleyFunction; + m_pawleyFunction = pawleyFunction; + m_pawleyFunction->fix(m_pawleyFunction->parameterIndex("f0.ZeroShift")); } DECLARE_FUNCTION(PoldiSpectrumPawleyFunction) diff --git a/Code/Mantid/Framework/SINQ/test/PoldiSpectrumPawleyFunctionTest.h b/Code/Mantid/Framework/SINQ/test/PoldiSpectrumPawleyFunctionTest.h index 5c61783146da7d191be00e7d24a4d63395f5eb57..aebc73942545167053f24ad199ee8069523f089e 100644 --- a/Code/Mantid/Framework/SINQ/test/PoldiSpectrumPawleyFunctionTest.h +++ b/Code/Mantid/Framework/SINQ/test/PoldiSpectrumPawleyFunctionTest.h @@ -5,7 +5,7 @@ #include "MantidSINQ/PoldiUtilities/PoldiSpectrumPawleyFunction.h" -using Mantid::SINQ::PoldiSpectrumPawleyFunction; +using namespace Mantid::Poldi; using namespace Mantid::API; class PoldiSpectrumPawleyFunctionTest : public CxxTest::TestSuite @@ -26,4 +26,4 @@ public: }; -#endif /* MANTID_SINQ_POLDISPECTRUMPAWLEYFUNCTIONTEST_H_ */ \ No newline at end of file +#endif /* MANTID_SINQ_POLDISPECTRUMPAWLEYFUNCTIONTEST_H_ */ diff --git a/Code/Mantid/instrument/POLDI_Parameters_2014b.xml b/Code/Mantid/instrument/POLDI_Parameters_2014b.xml index e7c318fa2fec41be006d7b9b9e5e1365a4f41bdf..4867cd326c1a49022812bc68dac0e4697ca39e43 100644 --- a/Code/Mantid/instrument/POLDI_Parameters_2014b.xml +++ b/Code/Mantid/instrument/POLDI_Parameters_2014b.xml @@ -3,22 +3,22 @@ <component-link name="chopper"> <parameter name="t0" type="number"> - <value val="-0.00801" /> + <value val="-0.0069" /> </parameter> <parameter name="t0_const" type="number"> - <value val="-18.3" /> + <value val="-6.9" /> </parameter> </component-link> <component-link name="detector"> <parameter name="two_theta" type="number"> - <value val="90.95" /> + <value val="91.12" /> </parameter> <parameter name="x"> - <value val="-1.000" /> + <value val="-0.8322" /> </parameter> <parameter name="y"> - <value val="-0.900" /> + <value val="-0.910" /> </parameter> </component-link>