diff --git a/Code/Mantid/Framework/API/inc/MantidAPI/IFunction.h b/Code/Mantid/Framework/API/inc/MantidAPI/IFunction.h index ef12bf8e9776e5f7aeda4c2cba2f1df7883976e5..e5640cd2ae8fe53ee88f9a8754898fe962e519c4 100644 --- a/Code/Mantid/Framework/API/inc/MantidAPI/IFunction.h +++ b/Code/Mantid/Framework/API/inc/MantidAPI/IFunction.h @@ -164,10 +164,12 @@ public: /// Convert a value from one unit (inUnit) to unit defined in workspace (ws) double convertValue(double value, Kernel::Unit_sptr& inUnit, boost::shared_ptr<const MatrixWorkspace> ws, - int wsIndex); -void convertValue(std::vector<double>& values, Kernel::Unit_sptr& outUnit, + int wsIndex) const; + + /// Convert a value from one unit (inUnit) to unit defined in workspace (ws) + void convertValue(std::vector<double>& values, Kernel::Unit_sptr& outUnit, boost::shared_ptr<const MatrixWorkspace> ws, - int wsIndex); + int wsIndex) const; protected: diff --git a/Code/Mantid/Framework/API/src/IFunction.cpp b/Code/Mantid/Framework/API/src/IFunction.cpp index 5878dd2d71f2f2972042fc2e12f733448969584b..48f721474767ca7f8ad8b1ec86d71aad062fdd7c 100644 --- a/Code/Mantid/Framework/API/src/IFunction.cpp +++ b/Code/Mantid/Framework/API/src/IFunction.cpp @@ -411,50 +411,15 @@ void IFunction::setMatrixWorkspace(boost::shared_ptr<const API::MatrixWorkspace> */ double IFunction::convertValue(double value, Kernel::Unit_sptr& outUnit, boost::shared_ptr<const MatrixWorkspace> ws, - int wsIndex) + int wsIndex) const { double retVal = value; - Kernel::Unit_sptr wsUnit = ws->getAxis(0)->unit(); - - // if unit required by formula or look-up-table different from ws-unit then - if ( outUnit->unitID().compare(wsUnit->unitID()) != 0 ) - { - // first check if it is possible to do a quick convertion convert - double factor,power; - if (wsUnit->quickConversion(*outUnit,factor,power) ) - { - retVal = factor * std::pow(retVal,power); - } - else - { - double l1,l2,twoTheta; - // Get l1, l2 and theta (see also RemoveBins.calculateDetectorPosition()) - IInstrument_const_sptr instrument = ws->getInstrument(); - Geometry::IObjComponent_const_sptr sample = instrument->getSample(); - l1 = instrument->getSource()->getDistance(*sample); - Geometry::IDetector_const_sptr det = ws->getDetector(wsIndex); - if ( ! det->isMonitor() ) - { - l2 = det->getDistance(*sample); - twoTheta = ws->detectorTwoTheta(det); - } - else // If this is a monitor then make l1+l2 = source-detector distance and twoTheta=0 - { - l2 = det->getDistance(*(instrument->getSource())); - l2 = l2 - l1; - twoTheta = 0.0; - } + std::vector<double> endPoint; + endPoint.push_back(value); + convertValue(endPoint, outUnit, ws, wsIndex); - std::vector<double> endPoint; - endPoint.push_back(retVal); - std::vector<double> emptyVec; - wsUnit->toTOF(endPoint,emptyVec,l1,l2,twoTheta,0,0.0,0.0); - outUnit->fromTOF(endPoint,emptyVec,l1,l2,twoTheta,0,0.0,0.0); - retVal = endPoint[0]; - } - } - return retVal; + return endPoint[0]; } /** Convert values from unit defined in workspace (ws) to outUnit @@ -467,7 +432,7 @@ double IFunction::convertValue(double value, Kernel::Unit_sptr& outUnit, */ void IFunction::convertValue(std::vector<double>& values, Kernel::Unit_sptr& outUnit, boost::shared_ptr<const MatrixWorkspace> ws, - int wsIndex) + int wsIndex) const { Kernel::Unit_sptr wsUnit = ws->getAxis(0)->unit(); @@ -502,8 +467,6 @@ void IFunction::convertValue(std::vector<double>& values, Kernel::Unit_sptr& out twoTheta = 0.0; } - //std::vector<double> endPoint; - //endPoint.push_back(retVal); std::vector<double> emptyVec; wsUnit->toTOF(values,emptyVec,l1,l2,twoTheta,0,0.0,0.0); outUnit->fromTOF(values,emptyVec,l1,l2,twoTheta,0,0.0,0.0); diff --git a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt index f3acc30789aa54320a093f0c20b00cf4848f0dd6..91aac83b0882a3e8d7b4590ac6ed37325b31d809 100644 --- a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt +++ b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt @@ -73,30 +73,30 @@ set ( INC_FILES inc/MantidCurveFitting/BackgroundFunction.h inc/MantidCurveFitting/UserFunction.h ) set ( TEST_FILES #test/BackToBackExponentialTest.h # TODO tests are commented out -# test/BoundaryConstraintTest.h -# test/ChebyshevTest.h -# test/CompositeFunctionTest.h -# test/ConvolutionTest.h -# test/CostFunctionFactoryTest.h -# test/ExpDecayTest.h -# test/FitTest.h -# test/FuncMinimizerFactoryTest.h + test/BoundaryConstraintTest.h + test/ChebyshevTest.h + test/CompositeFunctionTest.h + test/ConvolutionTest.h + test/CostFunctionFactoryTest.h + test/ExpDecayTest.h + test/FitTest.h + test/FuncMinimizerFactoryTest.h test/FunctionFactoryTest.h -# test/FunctionTest.h -# test/Gaussian1DTest.h -# test/GaussianTest.h # TODO tests are commented out -# #test/IkedaCarpenterPVTest.h # TODO tests are commented out -# test/LinearBackgroundTest.h -# test/LinearTest.h -# test/Lorentzian1DTest.h -# test/LorentzianTest.h -# test/PlotPeakByLogValueTest.h -# test/QuadraticTest.h -# test/ResolutionTest.h -# test/SpecialFunctionSupportTest.h -# test/SplineBackgroundTest.h -# test/UserFunction1DTest.h -# test/UserFunctionTest.h + test/FunctionTest.h + test/Gaussian1DTest.h + test/GaussianTest.h + test/IkedaCarpenterPVTest.h + test/LinearBackgroundTest.h + test/LinearTest.h + test/Lorentzian1DTest.h + test/LorentzianTest.h + test/PlotPeakByLogValueTest.h + test/QuadraticTest.h + test/ResolutionTest.h + test/SpecialFunctionSupportTest.h + test/SplineBackgroundTest.h + test/UserFunction1DTest.h + test/UserFunctionTest.h ) # Add the target for this directory diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/IkedaCarpenterPV.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/IkedaCarpenterPV.h index ab6627c2a633355921fb7ced1767eb00a7ffa0b5..8e4cfa8cccabb4aaff72615352fd24142712ff07 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/IkedaCarpenterPV.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/IkedaCarpenterPV.h @@ -88,9 +88,9 @@ namespace Mantid void calWavelengthAtEachDataPoint(const double* xValues, const int& nData) const; /// used for cal m_waveLengthFixed - void convertValue(std::vector<double>& values, Kernel::Unit_sptr& outUnit, - boost::shared_ptr<const API::MatrixWorkspace> ws, - int wsIndex) const; + //void convertValue(std::vector<double>& values, Kernel::Unit_sptr& outUnit, + // boost::shared_ptr<const API::MatrixWorkspace> ws, + // int wsIndex) const; /// convert voigt params to pseudo voigt params void convertVoigtToPseudo(const double& voigtSigmaSq, const double& voigtGamma, diff --git a/Code/Mantid/Framework/CurveFitting/src/IkedaCarpenterPV.cpp b/Code/Mantid/Framework/CurveFitting/src/IkedaCarpenterPV.cpp index 9e208941a130802d6c64650ec0df23602f786b12..31372a86f27dbeb45c15cfb594a5865dd39538e2 100644 --- a/Code/Mantid/Framework/CurveFitting/src/IkedaCarpenterPV.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/IkedaCarpenterPV.cpp @@ -12,7 +12,7 @@ #include <gsl/gsl_multifit_nlin.h> #include <limits> #include "MantidGeometry/Instrument/DetectorGroup.h" - +#include <limits> namespace Mantid { @@ -105,116 +105,26 @@ void IkedaCarpenterPV::calWavelengthAtEachDataPoint(const double* xValues, const // further we make the assumption that no need to recalculate this vector if // it already has the right size - /* if (static_cast<int>(m_waveLength.size()) != nData) { - // This peak shape requires the fit to be done in TOF units - - if ( m_workspace->getAxis(0)->unit()->unitID().compare("TOF") != 0 ) - { - g_log.information() << "IkedaCarpenterPV function is perhaps best used when working with x-axis unit = TOF\n"; - } - m_waveLength.resize(nData); - // Get the geometric information for this detector - - Geometry::IInstrument_const_sptr instrument = m_workspace->getInstrument(); - Geometry::IObjComponent_const_sptr sample = instrument->getSample(); - const double l1 = instrument->getSource()->getDistance(*sample); - Geometry::IDetector_sptr det = m_workspace->getDetector(m_workspaceIndex); // i is the workspace index - const double l2 = det->getDistance(*sample); - const double twoTheta = m_workspace->detectorTwoTheta(det); - - Mantid::Kernel::Unit_const_sptr wavelength = Mantid::Kernel::UnitFactory::Instance().create("Wavelength"); - for (int i = 0; i < nData; i++) - { - m_waveLength[i] = xValues[i]; - } - std::vector<double> y; // Create an empty vector, it's not used in fromTOF - wavelength->fromTOF(m_waveLength,y,l1,l2,twoTheta,0,0.0,0.0); - }*/ - - if (static_cast<int>(m_waveLength.size()) != nData) - { - // This peak shape requires the fit to be done in TOF units - - if ( m_workspace->getAxis(0)->unit()->unitID().compare("TOF") != 0 ) - { - g_log.information() << "IkedaCarpenterPV function is perhaps best used when working with x-axis unit = TOF\n"; - } - - m_waveLength.resize(nData); - - Mantid::Kernel::Unit_sptr wavelength = Mantid::Kernel::UnitFactory::Instance().create("Wavelength"); for (int i = 0; i < nData; i++) { - m_waveLength[i] = xValues[i]; + m_waveLength[i] = xValues[i]; } + + // note if a version of convertValue was added which allows a double* as first argument + // then could avoid copying above plus only have to resize m_wavelength when + // its size smaller than nData convertValue(m_waveLength, wavelength, m_workspace, m_workspaceIndex); + std::cout << "wave " << m_waveLength[0] << "\n"; } } } -/** Convert values from unit defined in workspace (ws) to outUnit - * - * @param values As input: assumed to be in unit of workspace. - * As output: in unit of outUnit - * @param outUnit unit to convert to - * @param ws workspace - * @param wsIndex workspace index - */ -void IkedaCarpenterPV::convertValue(std::vector<double>& values, Kernel::Unit_sptr& outUnit, - boost::shared_ptr<const API::MatrixWorkspace> ws, - int wsIndex) const -{ - Kernel::Unit_sptr wsUnit = ws->getAxis(0)->unit(); - - // if unit required by formula or look-up-table different from ws-unit then - if ( outUnit->unitID().compare(wsUnit->unitID()) != 0 ) - { - // first check if it is possible to do a quick convertion convert - double factor,power; - if (wsUnit->quickConversion(*outUnit,factor,power) ) - { - for (unsigned int i = 0; i < values.size(); i++) - values[i] = factor * std::pow(values[i],power); - } - else - { - double l1,l2,twoTheta; - - // Get l1, l2 and theta (see also RemoveBins.calculateDetectorPosition()) - Geometry::IInstrument_const_sptr instrument = ws->getInstrument(); - Geometry::IObjComponent_const_sptr sample = instrument->getSample(); - l1 = instrument->getSource()->getDistance(*sample); - Geometry::IDetector_sptr det = ws->getDetector(wsIndex); - if ( boost::dynamic_pointer_cast<Geometry::DetectorGroup>(det) ) - { - det = instrument->getDetector(det->getID()); - } - if ( ! det->isMonitor() ) - { - l2 = det->getDistance(*sample); - twoTheta = ws->detectorTwoTheta(det); - } - else // If this is a monitor then make l1+l2 = source-detector distance and twoTheta=0 - { - l2 = det->getDistance(*(instrument->getSource())); - l2 = l2 - l1; - twoTheta = 0.0; - } - - std::vector<double> emptyVec; - wsUnit->toTOF(values,emptyVec,l1,l2,twoTheta,0,0.0,0.0); - outUnit->fromTOF(values,emptyVec,l1,l2,twoTheta,0,0.0,0.0); - } - } -} - - /** convert voigt params to pseudo voigt params * * @param voigtSigmaSq voigt param @@ -235,6 +145,9 @@ void IkedaCarpenterPV::convertVoigtToPseudo(const double& voigtSigmaSq, const do H = pow(fwhmG4*fwhmG+2.69269*fwhmG4*fwhmL+2.42843*fwhmGsq*fwhmG*fwhmLsq +4.47163*fwhmGsq*fwhmLsq*fwhmL+0.07842*fwhmG*fwhmL4+fwhmL4*fwhmL, 0.2); + if (H == 0.0) + H = std::numeric_limits<double>::epsilon()*1000.0; + double tmp = fwhmL/H; eta = 1.36603*tmp - 0.47719*tmp*tmp + 0.11116*tmp*tmp*tmp; diff --git a/Code/Mantid/Framework/CurveFitting/test/IkedaCarpenterPVTest.h b/Code/Mantid/Framework/CurveFitting/test/IkedaCarpenterPVTest.h index 630860bfa166a16471c9fe154b263333fcd3d2d9..988309f4e42471cfedde96e6cf7362e8a642ca8d 100644 --- a/Code/Mantid/Framework/CurveFitting/test/IkedaCarpenterPVTest.h +++ b/Code/Mantid/Framework/CurveFitting/test/IkedaCarpenterPVTest.h @@ -18,10 +18,6 @@ #include "MantidDataHandling/LoadRaw.h" #include "MantidDataHandling/LoadRaw3.h" #include "MantidKernel/Exception.h" -#include "MantidAlgorithms/ConvertUnits.h" -#include "MantidAlgorithms/DiffractionFocussing2.h" -#include "MantidAlgorithms/CheckWorkspacesMatch.h" -#include "MantidAlgorithms/AlignDetectors.h" #include "MantidDataHandling/LoadInstrument.h" #include "MantidNexus/LoadNeXus.h" #include "MantidKernel/ConfigService.h" @@ -32,7 +28,6 @@ using namespace Mantid::API; using namespace Mantid::CurveFitting; using namespace Mantid::DataObjects; using namespace Mantid::DataHandling; -using namespace Mantid::Algorithms; using namespace Mantid::NeXus; class IkedaCarpenterPVTest : public CxxTest::TestSuite @@ -111,7 +106,7 @@ public: e[30] = 1.0112; } - void t1estAgainstMockData() + void testAgainstMockData() { Fit alg2; TS_ASSERT_THROWS_NOTHING(alg2.initialize()); @@ -183,148 +178,6 @@ public: } - // motivation for this test is to figure out way IC function goes absolutely - // nuts when a large data range are selection - void t1estAgainstGEM_dataLargeDataRange() - { - LoadNexus load; - load.initialize(); - load.setPropertyValue("FileName", "../../../../Test/AutoTestData/focussedGEM38370_TOF.nxs"); - std::string wsname = "GEM38370nexus"; - load.setPropertyValue("OutputWorkspace", wsname); - TS_ASSERT_THROWS_NOTHING(load.execute()); - TS_ASSERT( load.isExecuted() ); - - Mantid::DataObjects::Workspace2D_sptr wsToPass = boost::dynamic_pointer_cast<Mantid::DataObjects::Workspace2D>(AnalysisDataService::Instance().retrieve(wsname)); - - Fit alg2; - TS_ASSERT_THROWS_NOTHING(alg2.initialize()); - TS_ASSERT( alg2.isInitialized() ); - - // Set general Fit parameters - alg2.setPropertyValue("InputWorkspace", wsname); - alg2.setPropertyValue("WorkspaceIndex","1"); - alg2.setPropertyValue("StartX","5000"); - alg2.setPropertyValue("EndX","10000"); - - // create function you want to fit against - CompositeFunction *fnWithBk = new CompositeFunction(); - - LinearBackground *bk = new LinearBackground(); - bk->initialize(); - - bk->setParameter("A0",0.0); - bk->setParameter("A1",0.0); - bk->tie("A1", "0.0"); - - // set up fitting function and pass to Fit - IkedaCarpenterPV* icpv = new IkedaCarpenterPV(); - icpv->initialize(); - - icpv->setParameter("I",25094.45); - icpv->setParameter("X0",7316); - //icpv->setParameter("Gamma",1); - //icpv->tie("Gamma", "1"); - - icpv->setMatrixWorkspace(wsToPass, 1,0,1); // for unit testing purpose set workspace here - - TS_ASSERT_DELTA( icpv->getParameter("Alpha0"), 0.734079 ,0.001); - TS_ASSERT_DELTA( icpv->getParameter("Alpha1"), 2.067249 ,0.001); - TS_ASSERT_DELTA( icpv->getParameter("SigmaSquared"), 6403 ,1); - - - std::vector<double> testing; - for (double d=5000; d<=10000; d+=1000) - testing.push_back(d); - - std::vector<double> out=testing; - - icpv->function(&out[0], &testing[0], testing.size()); - - TS_ASSERT_DELTA( out[0], 0.2694,0.001); - - AnalysisDataService::Instance().remove(wsname); - // Append value of date-time tag inside the geometry file to the constructor handle - // for change to LoadInstrument - InstrumentDataService::Instance().remove("GEM_Definition.xml16th Sep 2008"); - } - - - void t1estAgainstGEM_data() - { - LoadNexus load; - load.initialize(); - load.setPropertyValue("FileName", "../../../../Test/AutoTestData/focussedGEM38370_TOF.nxs"); - std::string wsname = "GEM38370nexus"; - load.setPropertyValue("OutputWorkspace", wsname); - TS_ASSERT_THROWS_NOTHING(load.execute()); - TS_ASSERT( load.isExecuted() ); - - Mantid::DataObjects::Workspace2D_sptr wsToPass = boost::dynamic_pointer_cast<Mantid::DataObjects::Workspace2D>(AnalysisDataService::Instance().retrieve(wsname)); - - Fit alg2; - TS_ASSERT_THROWS_NOTHING(alg2.initialize()); - TS_ASSERT( alg2.isInitialized() ); - - // Set general Fit parameters - alg2.setPropertyValue("InputWorkspace", wsname); - alg2.setPropertyValue("WorkspaceIndex","1"); - alg2.setPropertyValue("StartX","6935.79"); - alg2.setPropertyValue("EndX","7682.56"); - - // create function you want to fit against - CompositeFunction *fnWithBk = new CompositeFunction(); - - LinearBackground *bk = new LinearBackground(); - bk->initialize(); - - bk->setParameter("A0",0.0); - bk->setParameter("A1",0.0); - bk->tie("A1", "0.0"); - - // set up fitting function and pass to Fit - IkedaCarpenterPV* icpv = new IkedaCarpenterPV(); - icpv->initialize(); - - icpv->setParameter("I",106860.45); - icpv->setParameter("X0",7326.34); - icpv->setParameter("Gamma",1); - icpv->tie("Gamma", "1"); - - icpv->setMatrixWorkspace(wsToPass, 1,0,1); // for unit testing purpose set workspace here - - TS_ASSERT_DELTA( icpv->getParameter("Alpha0"), 0.734079 ,0.001); - TS_ASSERT_DELTA( icpv->getParameter("Alpha1"), 2.067249 ,0.001); - TS_ASSERT_DELTA( icpv->getParameter("SigmaSquared"), 6422 ,1); - - fnWithBk->addFunction(icpv); - fnWithBk->addFunction(bk); - - alg2.setFunction(fnWithBk); - - // execute fit - TS_ASSERT_THROWS_NOTHING( - TS_ASSERT( alg2.execute() ) - ) - TS_ASSERT( alg2.isExecuted() ); - - // test the output from fit is what you expect - // note this test will never produce good fit because it assumes no background - double dummy = alg2.getProperty("Output Chi^2/DoF"); - TS_ASSERT_DELTA( dummy, 0.831,0.01); - - TS_ASSERT_DELTA( icpv->getParameter("I"), 69562 ,1); - TS_ASSERT_DELTA( icpv->getParameter("Alpha0"), 0.734079 ,0.1); - TS_ASSERT_DELTA( icpv->getParameter("Alpha1"), 2.067249 ,0.1); - TS_ASSERT_DELTA( icpv->getParameter("SigmaSquared"), 3567 ,1); - TS_ASSERT_DELTA( icpv->getParameter("X0"), 7301 ,1); - TS_ASSERT_DELTA( icpv->getParameter("Gamma"), 1 ,0.1); - TS_ASSERT_DELTA( bk->getParameter("A0"), 90.0 ,1); - TS_ASSERT_DELTA( bk->getParameter("A1"), 0.0 ,0.000000001); - - AnalysisDataService::Instance().remove(wsname); - } - }; #endif /*IKEDACARPENTERPVTEST_H_*/ diff --git a/Code/Mantid/Instrument/GEM_Definition.xml b/Code/Mantid/Instrument/GEM_Definition.xml index 96b02d4eedbbcbd864f0c162b852c1d68ddeff57..a564fe6e3354e1ae16a83bb103577d0f7b5209a9 100644 --- a/Code/Mantid/Instrument/GEM_Definition.xml +++ b/Code/Mantid/Instrument/GEM_Definition.xml @@ -49,7 +49,10 @@ <formula eq="2.067249" unit="TOF" result-unit="TOF"/> <fixed /> </parameter> - <parameter name="IkedaCarpenterPV:Kappa" type="fitting"> <value val="48.734158"/> <fixed /> </parameter> + <parameter name="IkedaCarpenterPV:Kappa" type="fitting"> + <value val="48.734158"/> + <fixed /> + </parameter> <!-- LIST OF PHYSICAL COMPONENTS (which the instrument consists of) --> @@ -63,7 +66,7 @@ <component type="bank1" idlist="bank1"> <location r="2.3735" t="0" p="0"> <facing val="none"/> </location> <parameter name="IkedaCarpenterPV:SigmaSquared" type="fitting"> - <formula eq="166.868*centre^2+98.757" unit="dSpacing"/> + <formula eq="166.868*centre^2+98.757" unit="dSpacing" result-unit="TOF^2"/> </parameter> <parameter name="IkedaCarpenterPV:Gamma" type="fitting"> <formula eq="0.277*centre" unit="dSpacing" result-unit="TOF"/>