diff --git a/Framework/Crystal/test/OptimizeCrystalPlacementTest.h b/Framework/Crystal/test/OptimizeCrystalPlacementTest.h index 582a25b2039fcc4ad448a83249b777f13d95b06e..89b510b2b32d40d1ed2e8e9548dcb033fb4cbf1b 100644 --- a/Framework/Crystal/test/OptimizeCrystalPlacementTest.h +++ b/Framework/Crystal/test/OptimizeCrystalPlacementTest.h @@ -1,39 +1,65 @@ -/* - * OptimizeCrystalPlacementTest.h - * - * Created on: Mar 22, 2013 - * Author: Ruth Mikkelson - */ - #ifndef OPTIMIZECRYSTALPLACEMENTTEST_H_ #define OPTIMIZECRYSTALPLACEMENTTEST_H_ + #include <cxxtest/TestSuite.h> +#include "MantidAPI/AnalysisDataService.h" #include "MantidCrystal/LoadIsawPeaks.h" +#include "MantidCrystal/LoadIsawUB.h" #include "MantidCrystal/OptimizeCrystalPlacement.h" #include "MantidCrystal/PeakHKLErrors.h" #include "MantidDataObjects/PeaksWorkspace.h" -#include "MantidAPI/IFunction.h" -#include "MantidCrystal/LoadIsawUB.h" -#include "MantidCrystal/IndexPeaks.h" -#include "MantidKernel/Matrix.h" -#include "MantidDataObjects/TableWorkspace.h" -#include "MantidAPI/AnalysisDataService.h" #include "MantidGeometry/Instrument/Goniometer.h" -#include "MantidKernel/V3D.h" -#include "MantidGeometry/Instrument/ParameterMap.h" -#include "MantidGeometry/IObjComponent.h" -#include "MantidAPI/ITableWorkspace.h" - -using namespace Mantid; -using namespace Crystal; +using Mantid::API::AnalysisDataService; +using Mantid::API::ITableWorkspace; +using Mantid::API::ITableWorkspace_sptr; +using Mantid::API::Workspace_sptr; +using Mantid::Crystal::PeakHKLErrors; +using Mantid::Crystal::OptimizeCrystalPlacement; +using Mantid::DataObjects::PeaksWorkspace; using Mantid::DataObjects::PeaksWorkspace_sptr; -using Mantid::DataObjects::TableWorkspace_sptr; -using Mantid::DataObjects::TableWorkspace_sptr; using Mantid::Geometry::Goniometer; -using Mantid::Geometry::ParameterMap; -using Mantid::Geometry::IObjComponent_const_sptr; -using Mantid::API::ITableWorkspace; +using Mantid::Geometry::Instrument; +using Mantid::Kernel::Matrix; +using Mantid::Kernel::V3D; + +namespace { +using StringPairVector = std::vector<std::pair<std::string, std::string>>; +using ResultType = std::pair<PeaksWorkspace_sptr, ITableWorkspace_sptr>; + +/** + * Run the OptimizeCrystalPlacement algorithm + * @param peaksWS Passed to the "PeaksWorkspace" property + * @param args Additional vector of arguments. + */ +ResultType +runOptimizePlacement(const PeaksWorkspace_sptr &peaksWS, + const StringPairVector &args = StringPairVector()) { + // This algorithm can not currently run without the ADS so add them and remove + // after + const std::string inputName("__runOptimizePlacement_PeaksWorkspace"), + modPeaksName("__runOptimizePlacement_ModPeaksWorkspace"), + fitTableName("__runOptimizePlacement_FitInfoTable"); + auto &ads = AnalysisDataService::Instance(); + ads.add(inputName, peaksWS); + + OptimizeCrystalPlacement alg; + alg.initialize(); + alg.setPropertyValue("PeaksWorkspace", inputName); + alg.setPropertyValue("ModifiedPeaksWorkspace", modPeaksName); + alg.setPropertyValue("FitInfoTable", fitTableName); + for (const auto &arg : args) { + alg.setPropertyValue(arg.first, arg.second); + } + alg.execute(); + auto modifiedPeaksWS = ads.retrieveWS<PeaksWorkspace>(modPeaksName); + auto fitInfoWS = ads.retrieveWS<ITableWorkspace>(fitTableName); + ads.remove(inputName); + ads.remove(modPeaksName); + ads.remove(fitTableName); + return std::make_pair(modifiedPeaksWS, fitInfoWS); +} +} class OptimizeCrystalPlacementTest : public CxxTest::TestSuite { @@ -47,223 +73,171 @@ public: delete suite; } - OptimizeCrystalPlacementTest() { initted = false; } - - void init() { - if (initted) - return; - initted = true; - LoadIsawPeaks alg; - alg.initialize(); - alg.setProperty("Filename", "TOPAZ_5637_8.peaks"); - alg.setProperty("OutputWorkspace", "abcd"); - alg.execute(); - - alg.setProperty("OutputWorkspace", "abcd"); - API::Workspace_sptr ows = alg.getProperty("OutputWorkspace"); - peaks = boost::dynamic_pointer_cast<DataObjects::PeaksWorkspace>(ows); - - LoadIsawUB loadUB; - loadUB.initialize(); - loadUB.setProperty("InputWorkspace", - alg.getPropertyValue("OutputWorkspace")); - loadUB.setProperty("Filename", "ls5637.mat"); - loadUB.execute(); - - // peaks->setName("abcd"); - } void test_basic() { - init(); - OptimizeCrystalPlacement alg; - alg.initialize(); - alg.setPropertyValue("PeaksWorkspace", "abcd"); - alg.setProperty("PeaksWorkspace", peaks); - alg.setPropertyValue("ModifiedPeaksWorkspace", "ModPeaks"); - alg.setPropertyValue("FitInfoTable", "FitInfoTable"); - alg.execute(); - - alg.setPropertyValue("ModifiedPeaksWorkspace", "ModPeaks"); - peaks1 = alg.getProperty("ModifiedPeaksWorkspace"); - alg.setPropertyValue("FitInfoTable", "FitInfoTable"); - boost::shared_ptr<ITableWorkspace> table = alg.getProperty("FitInfoTable"); - - Kernel::Matrix<double> Gon; - Gon.zeroMatrix(); - Kernel::Matrix<double> ZMat(Gon); - - Kernel::Matrix<double> Rotx = - PeakHKLErrors::RotationMatrixAboutRegAxis(0.5, 'x'); - Kernel::Matrix<double> Roty = - PeakHKLErrors::RotationMatrixAboutRegAxis(-1, 'y'); - - for (int i = 0; i < peaks1->getNumberPeaks(); ++i) - if (peaks1->getPeak(i).getRunNumber() == 5638) { - Geometry::IPeak &peak = peaks1->getPeak(i); - if (Gon == ZMat) { - origGon5638 = peak.getGoniometerMatrix(); - Gon = Rotx * Roty * origGon5638; + auto modPeaksNoFix = calculatePlacementUnfixed(); + Matrix<double> goniomMat, origGon5638Mat; + const auto rotx = PeakHKLErrors::RotationMatrixAboutRegAxis(0.5, 'x'); + const auto roty = PeakHKLErrors::RotationMatrixAboutRegAxis(-1, 'y'); + + const auto npeaks = modPeaksNoFix->getNumberPeaks(); + for (int i = 0; i < npeaks; ++i) { + if (modPeaksNoFix->getPeak(i).getRunNumber() == 5638) { + auto &peak = modPeaksNoFix->getPeak(i); + if (origGon5638Mat.numRows() == 0) { + origGon5638Mat = peak.getGoniometerMatrix(); + goniomMat = rotx * roty * origGon5638Mat; } - peak.setGoniometerMatrix(Gon); + peak.setGoniometerMatrix(goniomMat); } - Geometry::Goniometer GonInst(Gon); - std::vector<double> chiphiOmega = GonInst.getEulerAngles("YZY"); - - OptimizeCrystalPlacement alg1; - alg1.initialize(); - API::AnalysisDataService::Instance().addOrReplace("abcd1", peaks1); - - alg1.setPropertyValue("PeaksWorkspace", "abcd1"); - alg1.setPropertyValue("ModifiedPeaksWorkspace", "ModPeaks"); - alg1.setPropertyValue("FitInfoTable", "FitInfoTable1"); - alg1.setProperty("KeepGoniometerFixedfor", "5637"); - alg1.execute(); - - alg1.setPropertyValue("FitInfoTable", "FitInfoTable1"); - table = alg1.getProperty("FitInfoTable"); - - Geometry::Goniometer IGon(origGon5638); - std::vector<double> GonAngles5638 = IGon.getEulerAngles("YZY"); - - for (size_t i = 0; i < table->rowCount(); ++i) { - std::string nm = table->String(i, 0); + } + Goniometer instGoniom(goniomMat); + std::vector<double> chiphiOmega = instGoniom.getEulerAngles("YZY"); + auto resultsFix5637 = runOptimizePlacement( + modPeaksNoFix, {{"KeepGoniometerFixedfor", "5637"}}); + auto fitInfoFix5637 = resultsFix5637.second; + const auto angles5638 = Goniometer(origGon5638Mat).getEulerAngles("YZY"); + + for (size_t i = 0; i < fitInfoFix5637->rowCount(); ++i) { + std::string nm = fitInfoFix5637->String(i, 0); double d = 0.0; if (nm == "chi5638") - d = GonAngles5638[1] - table->Double(i, 1); + d = angles5638[1] - fitInfoFix5637->Double(i, 1); else if (nm == "phi5638") - d = GonAngles5638[2] - table->Double(i, 1); + d = angles5638[2] - fitInfoFix5637->Double(i, 1); else if (nm == "omega5638") - d = GonAngles5638[0] - table->Double(i, 1); + d = angles5638[0] - fitInfoFix5637->Double(i, 1); TS_ASSERT_DELTA(d, 0, .3); } } void test_tilt() { - init(); - Kernel::Matrix<double> tilt = - PeakHKLErrors::RotationMatrixAboutRegAxis(1, 'x') * - PeakHKLErrors::RotationMatrixAboutRegAxis(-2, 'y') * - PeakHKLErrors::RotationMatrixAboutRegAxis(1.3, 'z'); - - origGon5637.zeroMatrix(); - Kernel::Matrix<double> ZMat; - ZMat.zeroMatrix(); + auto modPeaksNoFix = calculatePlacementUnfixed(); + const auto rotx = PeakHKLErrors::RotationMatrixAboutRegAxis(0.5, 'x'); + const auto roty = PeakHKLErrors::RotationMatrixAboutRegAxis(-1, 'y'); + + Matrix<double> goniomMat, origGon5638Mat; + const auto npeaks = modPeaksNoFix->getNumberPeaks(); + for (int i = 0; i < npeaks; ++i) { + if (modPeaksNoFix->getPeak(i).getRunNumber() == 5638) { + auto &peak = modPeaksNoFix->getPeak(i); + if (origGon5638Mat.numRows() == 0) { + origGon5638Mat = peak.getGoniometerMatrix(); + goniomMat = rotx * roty * origGon5638Mat; + } + peak.setGoniometerMatrix(goniomMat); + } + } - for (int i = 0; i < peaks1->getNumberPeaks(); ++i) { - Geometry::IPeak &peak = peaks1->getPeak(i); - int RunNum = peak.getRunNumber(); + const auto tilt = PeakHKLErrors::RotationMatrixAboutRegAxis(1, 'x') * + PeakHKLErrors::RotationMatrixAboutRegAxis(-2, 'y') * + PeakHKLErrors::RotationMatrixAboutRegAxis(1.3, 'z'); - Kernel::Matrix<double> GG; - if (RunNum == 5637) { - if (origGon5637 == ZMat) - origGon5637 = peak.getGoniometerMatrix(); + Matrix<double> origGon5637Mat; + for (int i = 0; i < npeaks; ++i) { + auto &peak = modPeaksNoFix->getPeak(i); + const int runNum = peak.getRunNumber(); - peak.setGoniometerMatrix(tilt * origGon5637); + Matrix<double> GG; + if (runNum == 5637) { + if (origGon5637Mat.numRows() == 0) { + origGon5637Mat = peak.getGoniometerMatrix(); + } + peak.setGoniometerMatrix(tilt * origGon5637Mat); } else - peak.setGoniometerMatrix(tilt * origGon5638); + peak.setGoniometerMatrix(tilt * origGon5638Mat); } - Geometry::Goniometer IGon(origGon5638); - std::vector<double> GonAngles5638 = IGon.getEulerAngles("YZY"); - - Geometry::Goniometer IGon2(origGon5637); - std::vector<double> GonAngles5637 = IGon2.getEulerAngles("YZY"); - - OptimizeCrystalPlacement alg; - alg.initialize(); - API::AnalysisDataService::Instance().addOrReplace("abcd2", peaks1); - - alg.setPropertyValue("PeaksWorkspace", "abcd2"); - alg.setPropertyValue("ModifiedPeaksWorkspace", "ModPeaks"); - alg.setPropertyValue("FitInfoTable", "FitInfoTable2"); - alg.setProperty("KeepGoniometerFixedfor", "5637,5638"); - alg.setProperty("OptimizeGoniometerTilt", true); - alg.execute(); - alg.setPropertyValue("FitInfoTable", "FitInfoTable2"); - boost::shared_ptr<ITableWorkspace> table = alg.getProperty("FitInfoTable"); + const auto angles5638 = Goniometer(origGon5638Mat).getEulerAngles("YZY"); + const auto angles5637 = Goniometer(origGon5637Mat).getEulerAngles("YZY"); + auto resultsTiltFixBoth = runOptimizePlacement( + modPeaksNoFix, {{"KeepGoniometerFixedfor", "5637, 5638"}, + {"OptimizeGoniometerTilt", "1"}}); - Kernel::V3D Rotxyz; + const auto table = resultsTiltFixBoth.second; + V3D rotXYZ; for (size_t i = 0; i < table->rowCount(); ++i) { - std::string nm = table->String(i, 0); - + const std::string nm = table->String(i, 0); if (nm == "GonRotx") - Rotxyz[0] = table->Double(i, 1); - + rotXYZ[0] = table->Double(i, 1); else if (nm == "GonRoty") - Rotxyz[1] = table->Double(i, 1); - + rotXYZ[1] = table->Double(i, 1); else if (nm == "GonRotz") - Rotxyz[2] = table->Double(i, 1); + rotXYZ[2] = table->Double(i, 1); } - Kernel::Matrix<double> tilt2 = - PeakHKLErrors::RotationMatrixAboutRegAxis(Rotxyz[0], 'x') * - PeakHKLErrors::RotationMatrixAboutRegAxis(Rotxyz[1], 'y') * - PeakHKLErrors::RotationMatrixAboutRegAxis(Rotxyz[2], 'z'); - - Kernel::Matrix<double> Change = tilt2 * tilt; - - Geometry::Goniometer IGon3(Change * origGon5637); - std::vector<double> GonAngles5637a = IGon3.getEulerAngles("YZY"); - - Geometry::Goniometer IGon4(Change * origGon5638); - std::vector<double> GonAngles5638a = IGon4.getEulerAngles("YZY"); + const auto tilt2 = + PeakHKLErrors::RotationMatrixAboutRegAxis(rotXYZ[0], 'x') * + PeakHKLErrors::RotationMatrixAboutRegAxis(rotXYZ[1], 'y') * + PeakHKLErrors::RotationMatrixAboutRegAxis(rotXYZ[2], 'z'); + const auto tiltChange = tilt2 * tilt; + const auto angles5637a = + Goniometer(tiltChange * origGon5637Mat).getEulerAngles("YZY"); + const auto angles5638a = + Goniometer(tiltChange * origGon5638Mat).getEulerAngles("YZY"); for (int i = 0; i < 3; i++) { - TS_ASSERT_DELTA(GonAngles5637[i], GonAngles5637a[i], .2); - TS_ASSERT_DELTA(GonAngles5638[i], GonAngles5638a[i], .15); + TS_ASSERT_DELTA(angles5637[i], angles5637a[i], .2); + TS_ASSERT_DELTA(angles5638[i], angles5638a[i], .15); } } void test_SamplePosition() { - init(); - Geometry::IPeak &peak = peaks1->getPeak(0); - boost::shared_ptr<const Geometry::Instrument> Inst = peak.getInstrument(); - Kernel::V3D SampPos(.0003, -.00025, .00015); - - boost::shared_ptr<Geometry::ParameterMap> pmap = - Inst->getParameterMap(); // check if parameterized. - Geometry::IComponent_const_sptr sample = Inst->getSample(); - - pmap->addPositionCoordinate(sample.get(), "x", SampPos.X()); - pmap->addPositionCoordinate(sample.get(), "y", SampPos.Y()); - pmap->addPositionCoordinate(sample.get(), "z", SampPos.Z()); - boost::shared_ptr<const Geometry::Instrument> newInstr( - new Geometry::Instrument(Inst->baseInstrument(), pmap)); - - for (int i = 0; i < peaks1->getNumberPeaks(); ++i) - peaks1->getPeak(i).setInstrument(newInstr); - OptimizeCrystalPlacement alg; - alg.initialize(); - - API::AnalysisDataService::Instance().addOrReplace("abcd3", peaks1); - - alg.setPropertyValue("PeaksWorkspace", "abcd3"); - alg.setPropertyValue("ModifiedPeaksWorkspace", "ModPeaks"); - alg.setPropertyValue("FitInfoTable", "FitInfoTable2"); - alg.setProperty("KeepGoniometerFixedfor", "5637,5638"); - alg.setProperty("AdjustSampleOffsets", true); - alg.execute(); - - alg.setPropertyValue("FitInfoTable", "FitInfoTable2"); - boost::shared_ptr<ITableWorkspace> table = alg.getProperty("FitInfoTable"); + auto modPeaksNoFix = calculatePlacementUnfixed(); + const auto &peak = modPeaksNoFix->getPeak(0); + auto inst = peak.getInstrument(); + const V3D sampPos(.0003, -.00025, .00015); + + auto pmap = inst->getParameterMap(); + auto sample = inst->getSample(); + pmap->addPositionCoordinate(sample.get(), "x", sampPos.X()); + pmap->addPositionCoordinate(sample.get(), "y", sampPos.Y()); + pmap->addPositionCoordinate(sample.get(), "z", sampPos.Z()); + auto newInst = + boost::make_shared<const Instrument>(inst->baseInstrument(), pmap); + + for (int i = 0; i < modPeaksNoFix->getNumberPeaks(); ++i) { + modPeaksNoFix->getPeak(i).setInstrument(newInst); + } + // optimize + const auto resultsSamplePos = runOptimizePlacement( + modPeaksNoFix, {{"KeepGoniometerFixedfor", "5637, 5638"}, + {"AdjustSampleOffsets", "1"}}); + const auto table = resultsSamplePos.second; TS_ASSERT_DELTA(table->Double(0, 1), 0, .00024); TS_ASSERT_DELTA(table->Double(1, 1), 0, .00024); TS_ASSERT_DELTA(table->Double(2, 1), 0, .00024); - - /* for (size_t i = 0; i < table->rowCount(); ++i) - { - std::string nm = table->String(i, 0); - std::cout<<nm<<","<<table->Double(i,1)<<","<<table->Double(i,2)<<'\n'; - } - */ } private: - DataObjects::PeaksWorkspace_sptr peaks; - DataObjects::PeaksWorkspace_sptr peaks1; - Kernel::Matrix<double> origGon5637, origGon5638; - bool initted; + // Helper method to load peaks and optimize them in the basic + // case + PeaksWorkspace_sptr calculatePlacementUnfixed() { + auto peaksFileWS = loadTestPeaksInput(); + auto resultsNoFix = runOptimizePlacement(peaksFileWS); + return resultsNoFix.first; + } + + PeaksWorkspace_sptr loadTestPeaksInput() { + using Mantid::Crystal::LoadIsawPeaks; + using Mantid::Crystal::LoadIsawUB; + // Load peaks input file and UB + LoadIsawPeaks alg; + alg.setChild(true); + alg.initialize(); + alg.setProperty("Filename", "TOPAZ_5637_8.peaks"); + alg.setProperty("OutputWorkspace", "__"); + alg.execute(); + Workspace_sptr peaksWS = alg.getProperty("OutputWorkspace"); + LoadIsawUB loadUB; + loadUB.setChild(true); + loadUB.initialize(); + loadUB.setProperty("InputWorkspace", peaksWS); + loadUB.setProperty("Filename", "ls5637.mat"); + loadUB.execute(); + Workspace_sptr ows = loadUB.getProperty("InputWorkspace"); + return boost::dynamic_pointer_cast<PeaksWorkspace>(ows); + } }; #endif /* OPTIMIZECRYSTALPLACEMENTTEST_H_ */