diff --git a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h index a7217d53e063..b408a1fc3d38 100644 --- a/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h +++ b/Code/Mantid/Framework/SINQ/inc/MantidSINQ/PoldiFitPeaks2D.h @@ -75,7 +75,11 @@ class MANTID_SINQ_DLL PoldiFitPeaks2D : public API::Algorithm { Poldi2DFunction_sptr getFunctionPawley(std::string profileFunctionName, - const PoldiPeakCollection_sptr &peakCollection) const; + const PoldiPeakCollection_sptr &peakCollection); + + std::string getRefinedStartingCell(const std::string &initialCell, + const std::string &crystalSystem, + const PoldiPeakCollection_sptr &peakCollection); PoldiPeak_sptr getPeakFromPeakFunction(API::IPeakFunction_sptr profileFunction, @@ -99,7 +103,7 @@ class MANTID_SINQ_DLL PoldiFitPeaks2D : public API::Algorithm { PoldiPeakCollection_sptr getPeakCollectionFromFunction(const API::IFunction_sptr &fitFunction); Poldi2DFunction_sptr getFunctionFromPeakCollection( - const PoldiPeakCollection_sptr &peakCollection) const; + const PoldiPeakCollection_sptr &peakCollection); void addBackgroundTerms(Poldi2DFunction_sptr poldi2DFunction) const; API::IAlgorithm_sptr diff --git a/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp b/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp index d1bd62548f5d..8a6fca0481d1 100644 --- a/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp +++ b/Code/Mantid/Framework/SINQ/src/PoldiFitPeaks2D.cpp @@ -15,8 +15,11 @@ #include "MantidSINQ/PoldiUtilities/PoldiPeakCollection.h" #include "MantidSINQ/PoldiUtilities/PoldiInstrumentAdapter.h" #include "MantidSINQ/PoldiUtilities/PoldiDeadWireDecorator.h" + +#include "MantidAPI/ILatticeFunction.h" #include "MantidAPI/IPeakFunction.h" #include "MantidAPI/IPawleyFunction.h" +#include "MantidGeometry/Crystal/UnitCell.h" #include "MantidSINQ/PoldiUtilities/Poldi2DFunction.h" @@ -416,13 +419,13 @@ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionIndividualPeaks( * * 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. + * properties of the algorithm and wraps it in a Poldi2DFunction. + * + * The cell is refined using LatticeFunction to get better starting values. * * 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. + * 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. @@ -430,7 +433,7 @@ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionIndividualPeaks( */ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionPawley( std::string profileFunctionName, - const PoldiPeakCollection_sptr &peakCollection) const { + const PoldiPeakCollection_sptr &peakCollection) { Poldi2DFunction_sptr mdFunction(new Poldi2DFunction); boost::shared_ptr poldiPawleyFunction = @@ -452,7 +455,8 @@ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionPawley( pawleyFunction->setCrystalSystem(crystalSystem); std::string initialCell = getProperty("InitialCell"); - pawleyFunction->setUnitCell(initialCell); + pawleyFunction->setUnitCell( + getRefinedStartingCell(initialCell, crystalSystem, peakCollection)); IPeakFunction_sptr pFun = boost::dynamic_pointer_cast( FunctionFactory::Instance().createFunction(profileFunctionName)); @@ -474,6 +478,61 @@ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionPawley( return mdFunction; } +/** + * Tries to refine the initial cell using the supplied peaks + * + * This method tries to refine the initial unit cell using the indexed peaks + * that are supplied in the PoldiPeakCollection. If there are unindexed peaks, + * the cell will not be refined at all, instead the unmodified initial cell + * is returned. + * + * @param initialCell :: String with the initial unit cell + * @param crystalSystem :: Crystal system name + * @param peakCollection :: Collection of bragg peaks, must be indexed + * + * @return String for refined unit cell + */ +std::string PoldiFitPeaks2D::getRefinedStartingCell( + const std::string &initialCell, const std::string &crystalSystem, + const PoldiPeakCollection_sptr &peakCollection) { + + Geometry::UnitCell cell = Geometry::strToUnitCell(initialCell); + + ILatticeFunction_sptr latticeFunction = + boost::dynamic_pointer_cast( + FunctionFactory::Instance().createFunction("LatticeFunction")); + + latticeFunction->setCrystalSystem(crystalSystem); + latticeFunction->fix(latticeFunction->parameterIndex("ZeroShift")); + latticeFunction->setUnitCell(cell); + + // Remove errors from d-values + PoldiPeakCollection_sptr clone = peakCollection->clone(); + for (size_t i = 0; i < clone->peakCount(); ++i) { + PoldiPeak_sptr peak = clone->peak(i); + + // If there are unindexed peaks, don't refine, just return the initial cell + if(peak->hkl() == MillerIndices()) { + return initialCell; + } + + peak->setD(UncertainValue(peak->d().value())); + } + + TableWorkspace_sptr peakTable = clone->asTableWorkspace(); + + IAlgorithm_sptr fit = createChildAlgorithm("Fit"); + fit->setProperty("Function", + boost::static_pointer_cast(latticeFunction)); + fit->setProperty("InputWorkspace", peakTable); + fit->setProperty("CostFunction", "Unweighted least squares"); + fit->execute(); + + Geometry::UnitCell refinedCell = latticeFunction->getUnitCell(); + + return Geometry::unitCellToStr(refinedCell); +} + /** * Constructs a proper function from a peak collection * @@ -485,7 +544,7 @@ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionPawley( * @return Poldi2DFunction with one PoldiSpectrumDomainFunction per peak */ Poldi2DFunction_sptr PoldiFitPeaks2D::getFunctionFromPeakCollection( - const PoldiPeakCollection_sptr &peakCollection) const { + const PoldiPeakCollection_sptr &peakCollection) { std::string profileFunctionName = getProperty("PeakProfileFunction"); bool pawleyFit = getProperty("PawleyFit"); diff --git a/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h b/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h index 88f8e3f7f6a6..4eb4a00e1339 100644 --- a/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h +++ b/Code/Mantid/Framework/SINQ/test/PoldiFitPeaks2DTest.h @@ -322,6 +322,39 @@ class PoldiFitPeaks2DTest : public CxxTest::TestSuite TS_ASSERT_EQUALS(funConstant->nFunctions(), 1); } + void testGetRefinedStartingCell() { + // Get Silicon peaks for testing + PoldiPeakCollection_sptr peaks = PoldiPeakCollectionHelpers::createPoldiPeakCollectionNormalized(); + + TestablePoldiFitPeaks2D alg; + + std::string refinedCell; + TS_ASSERT_THROWS_NOTHING( + refinedCell = alg.getRefinedStartingCell( + "5.4 5.4 5.4 90 90 90", + "Cubic", peaks)); + + UnitCell cell = strToUnitCell(refinedCell); + TS_ASSERT_DELTA(cell.a(), 5.43111972, 1e-8); + + // With less accurate starting parameters the result should not change + TS_ASSERT_THROWS_NOTHING( + refinedCell = alg.getRefinedStartingCell( + "5 5 5 90 90 90", + "Cubic", peaks)); + + cell = strToUnitCell(refinedCell); + TS_ASSERT_DELTA(cell.a(), 5.43111972, 1e-8); + + // Adding an unindexed peak should make the function return the initial + peaks->addPeak(PoldiPeak::create(UncertainValue(1.0))); + TS_ASSERT_THROWS_NOTHING( + refinedCell = alg.getRefinedStartingCell( + "5 5 5 90 90 90", + "Cubic", peaks)); + TS_ASSERT_EQUALS(refinedCell, "5 5 5 90 90 90"); + } + private: PoldiInstrumentAdapter_sptr m_instrument; PoldiTimeTransformer_sptr m_timeTransformer;