diff --git a/Framework/API/CMakeLists.txt b/Framework/API/CMakeLists.txt index 4abf2d8319c87a4da5a48594c9337a0d8881f83c..b9e08b45bc5f22b14eb9498ab48112fcc8cbc34e 100644 --- a/Framework/API/CMakeLists.txt +++ b/Framework/API/CMakeLists.txt @@ -46,6 +46,7 @@ set ( SRC_FILES src/FunctionValues.cpp src/GridDomain.cpp src/GridDomain1D.cpp + src/GroupingLoader.cpp src/HistogramValidator.cpp src/HistoryItem.cpp src/HistoryView.cpp @@ -198,6 +199,7 @@ set ( INC_FILES inc/MantidAPI/FunctionValues.h inc/MantidAPI/GridDomain.h inc/MantidAPI/GridDomain1D.h + inc/MantidAPI/GroupingLoader.h inc/MantidAPI/HistogramValidator.h inc/MantidAPI/HistoryItem.h inc/MantidAPI/HistoryView.h @@ -356,6 +358,7 @@ set ( TEST_FILES FunctionPropertyTest.h FunctionTest.h FunctionValuesTest.h + GroupingLoaderTest.h HistogramValidatorTest.h HistoryItemTest.h HistoryViewTest.h diff --git a/Framework/API/inc/MantidAPI/GroupingLoader.h b/Framework/API/inc/MantidAPI/GroupingLoader.h new file mode 100644 index 0000000000000000000000000000000000000000..aad5270e0c67208b006cf0a63b0278d2abab4181 --- /dev/null +++ b/Framework/API/inc/MantidAPI/GroupingLoader.h @@ -0,0 +1,68 @@ +#ifndef MANTID_API_GROUPINGLOADER_H_ +#define MANTID_API_GROUPINGLOADER_H_ + +#include "MantidAPI/DllConfig.h" +#include "MantidGeometry/Instrument.h" + +namespace Mantid { +namespace API { + +/// Structure to represent grouping information +struct Grouping { + std::vector<std::string> groupNames; + std::vector<std::string> groups; // Range strings, e.g. "1-32" + + std::vector<std::string> pairNames; + std::vector<std::pair<size_t, size_t>> pairs; // Pairs of group ids + std::vector<double> pairAlphas; + + std::string description; + std::string defaultName; // Not storing id because can be either group or pair +}; + +/** GroupingLoader : Loads instrument grouping from IDF file + + Copyright © 2016 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_API_DLL GroupingLoader { +public: + explicit GroupingLoader(Geometry::Instrument_const_sptr instrument); + GroupingLoader(Geometry::Instrument_const_sptr instrument, + const std::string &mainFieldDirection); + virtual ~GroupingLoader(); + /// Load the grouping from the instrument's IDF + boost::shared_ptr<Grouping> getGroupingFromIDF() const; + /// Loads grouping from the XML file specified + static void loadGroupingFromXML(const std::string &filename, + Grouping &grouping); + +private: + /// Instrument to load grouping from + const Geometry::Instrument_const_sptr m_instrument; + /// Orientation of instrument (e.g. for MUSR) + const std::string m_mainFieldDirection; +}; + +} // namespace API +} // namespace Mantid + +#endif /* MANTID_API_GROUPINGLOADER_H_ */ \ No newline at end of file diff --git a/Framework/API/src/GroupingLoader.cpp b/Framework/API/src/GroupingLoader.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7a005a92523c3d0198fa9bb72e270bab3b5d17c7 --- /dev/null +++ b/Framework/API/src/GroupingLoader.cpp @@ -0,0 +1,220 @@ +#include "MantidAPI/GroupingLoader.h" +#include "MantidKernel/ConfigService.h" +#include "MantidKernel/Exception.h" +#include "MantidKernel/SingletonHolder.h" +#include <boost/make_shared.hpp> +#include <Poco/DOM/Document.h> +#include <Poco/DOM/DOMParser.h> +#include <Poco/DOM/NodeList.h> + +using namespace Poco::XML; + +namespace Mantid { +namespace API { + +//---------------------------------------------------------------------------------------------- +/** Constructor without field direction +* @param instrument :: [input] Instrument +*/ +GroupingLoader::GroupingLoader(Geometry::Instrument_const_sptr instrument) + : m_instrument(instrument) {} + +/** Constructor with field direction +* @param instrument :: [input] Instrument +* @param mainFieldDirection :: [input] Direction of main field (for MUSR) +*/ +GroupingLoader::GroupingLoader(Geometry::Instrument_const_sptr instrument, + const std::string &mainFieldDirection) + : m_instrument(instrument), m_mainFieldDirection(mainFieldDirection) {} + +//---------------------------------------------------------------------------------------------- +/** Destructor + */ +GroupingLoader::~GroupingLoader() {} + +/** + * Attempts to load a grouping information referenced by IDF. + * @return Grouping information + */ +boost::shared_ptr<Grouping> GroupingLoader::getGroupingFromIDF() const { + std::string parameterName = "Default grouping file"; + auto loadedGrouping = boost::make_shared<Grouping>(); + + // Special case for MUSR, because it has two possible groupings + if (m_instrument->getName() == "MUSR") { + parameterName.append(" - " + m_mainFieldDirection); + } + + std::vector<std::string> groupingFiles = + m_instrument->getStringParameter(parameterName); + + if (groupingFiles.size() == 1) { + const std::string groupingFile = groupingFiles[0]; + + // Get search directory for XML instrument definition files (IDFs) + std::string directoryName = + Kernel::ConfigService::Instance().getInstrumentDirectory(); + + loadGroupingFromXML(directoryName + groupingFile, *loadedGrouping); + } else { + throw std::runtime_error("Multiple groupings specified for the instrument"); + } + return loadedGrouping; +} + +/** + * Loads grouping from the XML file specified. + * + * @param filename :: XML filename to load grouping information from + * @param grouping :: Struct to store grouping information to + */ +void GroupingLoader::loadGroupingFromXML(const std::string &filename, + Grouping &grouping) { + // Set up the DOM parser and parse xml file + DOMParser pParser; + Poco::AutoPtr<Document> pDoc; + try { + pDoc = pParser.parse(filename); + } catch (...) { + throw Mantid::Kernel::Exception::FileError("Unable to parse File", + filename); + } + + // Get pointer to root element + Element *pRootElem = pDoc->documentElement(); + if (!pRootElem->hasChildNodes()) + throw Mantid::Kernel::Exception::FileError( + "No root element in XML grouping file", filename); + + // Parse information for groups + Poco::AutoPtr<NodeList> groups = pRootElem->getElementsByTagName("group"); + if (groups->length() == 0) + throw Mantid::Kernel::Exception::FileError( + "No groups specified in XML grouping file", filename); + + // Resize vectors + grouping.groupNames.resize(groups->length()); + grouping.groups.resize(groups->length()); + + for (size_t ig = 0; ig < groups->length(); ig++) { + Element *pGroupElem = + static_cast<Element *>(groups->item(static_cast<long>(ig))); + + if (!pGroupElem->hasAttribute("name")) + throw Mantid::Kernel::Exception::FileError("Group element without name", + filename); + + grouping.groupNames[ig] = pGroupElem->getAttribute("name"); + + Element *idlistElement = pGroupElem->getChildElement("ids"); + if (!idlistElement) + throw Mantid::Kernel::Exception::FileError("Group element without <ids>", + filename); + + grouping.groups[ig] = idlistElement->getAttribute("val"); + } + + // Parse information for pairs + Poco::AutoPtr<NodeList> pairs = pRootElem->getElementsByTagName("pair"); + + // Resize vectors + grouping.pairNames.resize(pairs->length()); + grouping.pairs.resize(pairs->length()); + grouping.pairAlphas.resize(pairs->length()); + + for (size_t ip = 0; ip < pairs->length(); ip++) { + Element *pPairElem = + static_cast<Element *>(pairs->item(static_cast<long>(ip))); + + if (!pPairElem->hasAttribute("name")) + throw Mantid::Kernel::Exception::FileError("Pair element without name", + filename); + + grouping.pairNames[ip] = pPairElem->getAttribute("name"); + + size_t fwdGroupId, bwdGroupId; // Ids of forward/backward groups + + // Try to get id of the first group + if (Element *fwdElement = pPairElem->getChildElement("forward-group")) { + if (!fwdElement->hasAttribute("val")) + throw Mantid::Kernel::Exception::FileError( + "Pair forward-group without <val>", filename); + + // Find the group with the given name + auto it = + std::find(grouping.groupNames.begin(), grouping.groupNames.end(), + fwdElement->getAttribute("val")); + + if (it == grouping.groupNames.end()) + throw Mantid::Kernel::Exception::FileError( + "Pair forward-group name not recognized", filename); + + // Get index of the iterator + fwdGroupId = it - grouping.groupNames.begin(); + } else { + throw Mantid::Kernel::Exception::FileError( + "Pair element without <forward-group>", filename); + } + + // Try to get id of the second group + if (Element *bwdElement = pPairElem->getChildElement("backward-group")) { + if (!bwdElement->hasAttribute("val")) + throw Mantid::Kernel::Exception::FileError( + "Pair backward-group without <val>", filename); + + // Find the group with the given name + auto it = + std::find(grouping.groupNames.begin(), grouping.groupNames.end(), + bwdElement->getAttribute("val")); + + if (it == grouping.groupNames.end()) + throw Mantid::Kernel::Exception::FileError( + "Pair backward-group name not recognized", filename); + + // Get index of the iterator + bwdGroupId = it - grouping.groupNames.begin(); + } else { + throw Mantid::Kernel::Exception::FileError( + "Pair element without <backward-group>", filename); + } + + grouping.pairs[ip] = std::make_pair(fwdGroupId, bwdGroupId); + + // Try to get alpha element + if (Element *aElement = pPairElem->getChildElement("alpha")) { + if (!aElement->hasAttribute("val")) + throw Mantid::Kernel::Exception::FileError( + "Pair alpha element with no <val>", filename); + + try // ... to convert value to double + { + grouping.pairAlphas[ip] = + boost::lexical_cast<double>(aElement->getAttribute("val")); + } catch (boost::bad_lexical_cast &) { + throw Mantid::Kernel::Exception::FileError( + "Pair alpha value is not a number", filename); + } + } + // If alpha element not there, default it to 1.0 + else { + grouping.pairAlphas[ip] = 1.0; + } + } + + // Try to get description + if (pRootElem->hasAttribute("description")) { + grouping.description = pRootElem->getAttribute("description"); + } + + // Try to get default group/pair name + if (Element *defaultElement = pRootElem->getChildElement("default")) { + if (!defaultElement->hasAttribute("name")) + throw Mantid::Kernel::Exception::FileError( + "Default element with no <name>", filename); + + grouping.defaultName = defaultElement->getAttribute("name"); + } +} + +} // namespace API +} // namespace Mantid diff --git a/Framework/API/test/GroupingLoaderTest.h b/Framework/API/test/GroupingLoaderTest.h new file mode 100644 index 0000000000000000000000000000000000000000..c906f15ac912f299a9a062c879865109803521d9 --- /dev/null +++ b/Framework/API/test/GroupingLoaderTest.h @@ -0,0 +1,78 @@ +#ifndef MANTID_API_GROUPINGLOADERTEST_H_ +#define MANTID_API_GROUPINGLOADERTEST_H_ + +#include <cxxtest/TestSuite.h> +#include <Poco/Path.h> +#include <Poco/File.h> +#include "MantidAPI/FrameworkManager.h" +#include "MantidAPI/GroupingLoader.h" + +using Mantid::API::GroupingLoader; +using Mantid::API::Grouping; + +class GroupingLoaderTest : 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 GroupingLoaderTest *createSuite() { return new GroupingLoaderTest(); } + static void destroySuite(GroupingLoaderTest *suite) { delete suite; } + + /// Constructor + GroupingLoaderTest() { + using Mantid::Kernel::ConfigService; + + auto dataPaths = ConfigService::Instance().getDataSearchDirs(); + + // Find the path of AutoTestData + for (auto it = dataPaths.begin(); it != dataPaths.end(); ++it) { + Poco::Path path(*it); + + if (path.directory(path.depth() - 1) == "UnitTest") { + m_testDataDir = *it; + break; + } + } + + TSM_ASSERT("Unable to find UnitTest data directory", + !m_testDataDir.empty()); + + m_tmpDir = ConfigService::Instance().getTempDir(); + + // To make sure API is initialized properly + Mantid::API::FrameworkManager::Instance(); + } + + void test_loadGroupingFromXML() { + Grouping g; + + TS_ASSERT_THROWS_NOTHING(GroupingLoader::loadGroupingFromXML( + m_testDataDir + "MUSRGrouping.xml", g)); + + TS_ASSERT_EQUALS(g.groupNames.size(), 2); + TS_ASSERT_EQUALS(g.groupNames[0], "fwd"); + TS_ASSERT_EQUALS(g.groupNames[1], "bwd"); + + TS_ASSERT_EQUALS(g.groups.size(), 2); + TS_ASSERT_EQUALS(g.groups[0], "33-64"); + TS_ASSERT_EQUALS(g.groups[1], "1-32"); + + TS_ASSERT_EQUALS(g.pairNames.size(), 1); + TS_ASSERT_EQUALS(g.pairNames[0], "long"); + + TS_ASSERT_EQUALS(g.pairs.size(), 1); + TS_ASSERT_EQUALS(g.pairs[0].first, 0); + TS_ASSERT_EQUALS(g.pairs[0].second, 1); + + TS_ASSERT_EQUALS(g.pairAlphas.size(), 1); + TS_ASSERT_EQUALS(g.pairAlphas[0], 1); + + TS_ASSERT_EQUALS(g.description, "musr longitudinal (64 detectors)"); + TS_ASSERT_EQUALS(g.defaultName, "long"); + } + +private: + std::string m_testDataDir; + std::string m_tmpDir; +}; + +#endif /* MANTID_API_GROUPINGLOADERTEST_H_ */ \ No newline at end of file diff --git a/Framework/Algorithms/inc/MantidAlgorithms/CalMuonDetectorPhases.h b/Framework/Algorithms/inc/MantidAlgorithms/CalMuonDetectorPhases.h index f62cd5d280d49af8e1906b7f25a079cad8151c46..a1ba06f68a8f6f70a6c48e2c461be6a07c855ec0 100644 --- a/Framework/Algorithms/inc/MantidAlgorithms/CalMuonDetectorPhases.h +++ b/Framework/Algorithms/inc/MantidAlgorithms/CalMuonDetectorPhases.h @@ -58,19 +58,49 @@ private: void exec() override; /// Validate the inputs std::map<std::string, std::string> validateInputs() override; - /// Prepare workspace for fit + /// Prepare workspace for fit by extracting data + API::MatrixWorkspace_sptr extractDataFromWorkspace(double startTime, + double endTime); + /// Remove exponential data from workspace API::MatrixWorkspace_sptr - prepareWorkspace(const API::MatrixWorkspace_sptr &ws, double startTime, - double endTime); + removeExpDecay(const API::MatrixWorkspace_sptr &wsInput); /// Fit the workspace void fitWorkspace(const API::MatrixWorkspace_sptr &ws, double freq, std::string groupName, API::ITableWorkspace_sptr &resTab, API::WorkspaceGroup_sptr &resGroup); /// Create the fitting function as string - std::string createFittingFunction(int nspec, double freq); + std::string createFittingFunction(double freq, bool fixFreq); /// Extract asymmetry and phase from fitting results - API::ITableWorkspace_sptr - extractDetectorInfo(const API::ITableWorkspace_sptr ¶mTab, size_t nspec); + void extractDetectorInfo(const API::ITableWorkspace_sptr ¶mTab, + const API::ITableWorkspace_sptr &resultsTab, + const int ispec); + /// Find frequency to use in sequential fit + double getFrequency(const API::MatrixWorkspace_sptr &ws); + /// Get frequency hint to use when finding frequency + double getFrequencyHint() const; + /// Get start time for fit + double getStartTime() const; + /// Get end time for fit + double getEndTime() const; + /// Calculate detector efficiency (alpha) + double getAlpha(const API::MatrixWorkspace_sptr &ws, + const std::vector<int> &forward, + const std::vector<int> &backward); + /// Calculate asymmetry + API::MatrixWorkspace_sptr getAsymmetry(const API::MatrixWorkspace_sptr &ws, + const std::vector<int> &forward, + const std::vector<int> &backward, + const double alpha); + /// Fit asymmetry to get frequency + double fitFrequencyFromAsymmetry(const API::MatrixWorkspace_sptr &wsAsym); + /// Find the grouping from the instrument + void getGroupingFromInstrument(const API::MatrixWorkspace_sptr &ws, + std::vector<int> &forward, + std::vector<int> &backward); + /// Report progress in GUI + void reportProgress(const int thisSpectrum, const int totalSpectra); + /// Pointer to input workspace + API::MatrixWorkspace_sptr m_inputWS; }; } // namespace Algorithms } // namespace Mantid diff --git a/Framework/Algorithms/src/CalMuonDetectorPhases.cpp b/Framework/Algorithms/src/CalMuonDetectorPhases.cpp index f1e84df270601a669b8805733be7b2920ca8c53e..8615381ca2bf56f84ff006d32e2bf095ec6a21f2 100644 --- a/Framework/Algorithms/src/CalMuonDetectorPhases.cpp +++ b/Framework/Algorithms/src/CalMuonDetectorPhases.cpp @@ -2,11 +2,14 @@ #include "MantidAPI/Axis.h" #include "MantidAPI/FunctionFactory.h" +#include "MantidAPI/GroupingLoader.h" #include "MantidAPI/ITableWorkspace.h" #include "MantidAPI/IFunction.h" #include "MantidAPI/MatrixWorkspace.h" #include "MantidAPI/MultiDomainFunction.h" #include "MantidAPI/TableRow.h" +#include "MantidKernel/ArrayProperty.h" +#include "MantidKernel/PhysicalConstants.h" #include "MantidAPI/WorkspaceFactory.h" namespace Mantid { @@ -46,6 +49,14 @@ void CalMuonDetectorPhases::init() { declareProperty(new API::WorkspaceProperty<API::WorkspaceGroup>( "DataFitted", "", Direction::Output), "Name of the output workspace holding fitting results"); + + declareProperty(new ArrayProperty<int>("ForwardSpectra", Direction::Input), + "The spectra numbers of the forward group. If not specified " + "will read from file."); + + declareProperty(new ArrayProperty<int>("BackwardSpectra", Direction::Input), + "The spectra numbers of the backward group. If not specified " + "will read from file."); } //---------------------------------------------------------------------------------------------- @@ -64,6 +75,21 @@ std::map<std::string, std::string> CalMuonDetectorPhases::validateInputs() { result["InputWorkspace"] = "InputWorkspace units must be microseconds"; } + // Check spectra numbers are valid, if specified + int nspec = static_cast<int>(inputWS->getNumberHistograms()); + std::vector<int> forward = getProperty("ForwardSpectra"); + std::vector<int> backward = getProperty("BackwardSpectra"); + for (int spec : forward) { + if (spec < 1 || spec > nspec) { + result["ForwardSpectra"] = "Invalid spectrum numbers in ForwardSpectra"; + } + } + for (int spec : backward) { + if (spec < 1 || spec > nspec) { + result["BackwardSpectra"] = "Invalid spectrum numbers in BackwardSpectra"; + } + } + return result; } //---------------------------------------------------------------------------------------------- @@ -72,62 +98,29 @@ std::map<std::string, std::string> CalMuonDetectorPhases::validateInputs() { void CalMuonDetectorPhases::exec() { // Get the input ws - API::MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace"); + m_inputWS = getProperty("InputWorkspace"); // Get start and end time - double startTime = getProperty("FirstGoodData"); - double endTime = getProperty("LastGoodData"); + double startTime = getStartTime(); + double endTime = getEndTime(); - // If startTime is EMPTY_DBL(): - if (startTime == EMPTY_DBL()) { - try { - // Read FirstGoodData from workspace logs if possible - double firstGoodData = - inputWS->run().getLogAsSingleValue("FirstGoodData"); - startTime = firstGoodData; - } catch (...) { - g_log.warning("Couldn't read FirstGoodData, setting to 0"); - // Set startTime to 0 - startTime = 0.; - } - } - // If endTime is EMPTY_DBL(): - if (endTime == EMPTY_DBL()) { - // Last available time - endTime = inputWS->readX(0).back(); - } + // Prepares the workspaces: extracts data from [startTime, endTime] + API::MatrixWorkspace_sptr tempWS = + extractDataFromWorkspace(startTime, endTime); // Get the frequency - double freq = getProperty("Frequency"); - - // If frequency is EMPTY_DBL(): - if (freq == EMPTY_DBL()) { - try { - // Read sample_magn_field from workspace logs - freq = inputWS->run().getLogAsSingleValue("sample_magn_field"); - // Multiply by muon gyromagnetic ratio: 0.01355 MHz/G - freq *= 2 * M_PI * 0.01355; - } catch (...) { - throw std::runtime_error( - "Couldn't read sample_magn_field. Please provide a value for " - "the frequency"); - } - } - - // Prepares the workspaces: extracts data from [startTime, endTime] and - // removes exponential decay - API::MatrixWorkspace_sptr tempWS = - prepareWorkspace(inputWS, startTime, endTime); + double freq = getFrequency(tempWS); // Create the output workspaces auto tab = API::WorkspaceFactory::Instance().createTable("TableWorkspace"); - auto group = API::WorkspaceGroup_sptr(new API::WorkspaceGroup()); + auto group = boost::make_shared<API::WorkspaceGroup>(); // Get the name of 'DataFitted' std::string groupName = getPropertyValue("DataFitted"); - // Fit the workspace - fitWorkspace(tempWS, freq, groupName, tab, group); + // Remove exponential decay and fit the workspace + auto wsToFit = removeExpDecay(tempWS); + fitWorkspace(wsToFit, freq, groupName, tab, group); // Set the table setProperty("DetectorTable", tab); @@ -150,162 +143,388 @@ void CalMuonDetectorPhases::fitWorkspace(const API::MatrixWorkspace_sptr &ws, int nspec = static_cast<int>(ws->getNumberHistograms()); // Create the fitting function f(x) = A * sin ( w * x + p ) - std::string funcStr = createFittingFunction(nspec, freq); - // Create the function from string - API::IFunction_sptr func = - API::FunctionFactory::Instance().createInitialized(funcStr); - // Create the multi domain function - boost::shared_ptr<API::MultiDomainFunction> multi = - boost::dynamic_pointer_cast<API::MultiDomainFunction>(func); - // Set the domain indices - for (int i = 0; i < nspec; ++i) { - multi->setDomainIndex(i, i); - } + // The same function and initial parameters are used for each fit + std::string funcStr = createFittingFunction(freq, true); + + // Set up results table + resTab->addColumn("int", "Detector"); + resTab->addColumn("double", "Asymmetry"); + resTab->addColumn("double", "Phase"); + + // Loop through fitting all spectra individually + const static std::string success = "success"; + for (int ispec = 0; ispec < nspec; ispec++) { + reportProgress(ispec, nspec); + auto fit = createChildAlgorithm("Fit"); + fit->initialize(); + fit->setPropertyValue("Function", funcStr); + fit->setProperty("InputWorkspace", ws); + fit->setProperty("WorkspaceIndex", ispec); + fit->setProperty("CreateOutput", true); + fit->setPropertyValue("Output", groupName); + fit->execute(); + + std::string status = fit->getProperty("OutputStatus"); + if (!fit->isExecuted() || status != success) { + std::ostringstream error; + error << "Fit failed for spectrum " << ispec; + throw std::runtime_error(error.str()); + } - API::IAlgorithm_sptr fit = createChildAlgorithm("Fit"); - fit->initialize(); - fit->addObserver(this->progressObserver()); - setChildStartProgress(0.); - setChildEndProgress(1.); - fit->setProperty("Function", - boost::dynamic_pointer_cast<API::IFunction>(multi)); - fit->setProperty("InputWorkspace", ws); - fit->setProperty("WorkspaceIndex", 0); - for (int s = 1; s < nspec; s++) { - std::string suffix = boost::lexical_cast<std::string>(s); - fit->setProperty("InputWorkspace_" + suffix, ws); - fit->setProperty("WorkspaceIndex_" + suffix, s); + API::MatrixWorkspace_sptr fitOut = fit->getProperty("OutputWorkspace"); + resGroup->addWorkspace(fitOut); + API::ITableWorkspace_sptr tab = fit->getProperty("OutputParameters"); + // Now we have our fitting results stored in tab + // but we need to extract the relevant information, i.e. + // the detector phases (parameter 'p') and asymmetries ('A') + extractDetectorInfo(tab, resTab, ispec); } - fit->setProperty("CreateOutput", true); - fit->setProperty("Output", groupName); - fit->execute(); - - // Get the parameter table - API::ITableWorkspace_sptr tab = fit->getProperty("OutputParameters"); - // Now we have our fitting results stored in tab - // but we need to extract the relevant information, i.e. - // the detector phases (parameter 'p') and asymmetries ('A') - resTab = extractDetectorInfo(tab, static_cast<size_t>(nspec)); - - // Get the fitting results - resGroup = fit->getProperty("OutputWorkspace"); } /** Extracts detector asymmetries and phases from fitting results +* and adds a new row to the results table with them * @param paramTab :: [input] Output parameter table resulting from the fit -* @param nspec :: [input] Number of detectors/spectra -* @return :: A new table workspace storing the asymmetries and phases +* @param resultsTab :: [input] Results table to update with a new row +* @param ispec :: [input] Spectrum number */ -API::ITableWorkspace_sptr CalMuonDetectorPhases::extractDetectorInfo( - const API::ITableWorkspace_sptr ¶mTab, size_t nspec) { - - // Make sure paramTable is the right size - // It should contain three parameters per detector/spectrum plus 'const - // function value' - if (paramTab->rowCount() != nspec * 3 + 1) { - throw std::invalid_argument( - "Can't extract detector parameters from fit results"); +void CalMuonDetectorPhases::extractDetectorInfo( + const API::ITableWorkspace_sptr ¶mTab, + const API::ITableWorkspace_sptr &resultsTab, const int ispec) { + + double asym = paramTab->Double(0, 1); + double phase = paramTab->Double(2, 1); + // If asym<0, take the absolute value and add \pi to phase + // f(x) = A * sin( w * x + p) = -A * sin( w * x + p + PI) + if (asym < 0) { + asym = -asym; + phase = phase + M_PI; } - - // Create the table to store detector info - auto tab = API::WorkspaceFactory::Instance().createTable("TableWorkspace"); - tab->addColumn("int", "Detector"); - tab->addColumn("double", "Asymmetry"); - tab->addColumn("double", "Phase"); - - // Reference frequency, all w values should be the same - double omegaRef = paramTab->Double(1, 1); - - for (size_t s = 0; s < nspec; s++) { - // The following '3' factor corresponds to the number of function params - size_t specRow = s * 3; - double asym = paramTab->Double(specRow, 1); - double omega = paramTab->Double(specRow + 1, 1); - double phase = paramTab->Double(specRow + 2, 1); - // If omega != omegaRef something went wrong with the fit - if (omega != omegaRef) { - throw std::runtime_error("Fit failed"); - } - // If asym<0, take the absolute value and add \pi to phase - // f(x) = A * sin( w * x + p) = -A * sin( w * x + p + PI) - if (asym < 0) { - asym = -asym; - phase = phase + M_PI; - } - // Now convert phases to interval [0, 2PI) - int factor = static_cast<int>(floor(phase / 2 / M_PI)); - if (factor) { - phase = phase - factor * 2 * M_PI; - } - // Copy parameters to new table - API::TableRow row = tab->appendRow(); - row << static_cast<int>(s) << asym << phase; + // Now convert phases to interval [0, 2PI) + int factor = static_cast<int>(floor(phase / 2 / M_PI)); + if (factor) { + phase = phase - factor * 2 * M_PI; } - - return tab; + // Copy parameters to new row in results table + API::TableRow row = resultsTab->appendRow(); + row << ispec << asym << phase; } -/** Creates the fitting function f(x) = A * sin( w*x + p) as string -* @param nspec :: [input] The number of domains (spectra) -* @param freq :: [input] Hint for the frequency (w) +/** Creates the fitting function f(x) = A * sin( w*x + p) + B as string +* Two modes: +* 1) Fixed frequency, no background - for main sequential fit +* 2) Varying frequency, flat background - for finding frequency from asymmetry +* @param freq :: [input] Value for the frequency (w) +* @param fixFreq :: [input] True: fixed frequency, no background. False: varying +* frequency with flat background. * @returns :: The fitting function as a string */ -std::string CalMuonDetectorPhases::createFittingFunction(int nspec, - double freq) { - +std::string CalMuonDetectorPhases::createFittingFunction(double freq, + bool fixFreq) { // The fitting function is: - // f(x) = A * sin ( w * x + p ) - // where w is shared across workspaces + // f(x) = A * sin ( w * x + p ) [+ B] std::ostringstream ss; - ss << "composite=MultiDomainFunction,NumDeriv=true;"; - for (int s = 0; s < nspec; s++) { - ss << "name=UserFunction,"; + ss << "name=UserFunction,"; + if (fixFreq) { + // no background ss << "Formula=A*sin(w*x+p),"; - ss << "A=0.5,"; - ss << "w=" << freq << ","; - ss << "p=0.;"; + } else { + // flat background + ss << "Formula=A*sin(w*x+p)+B,"; + ss << "B=0.5,"; } - ss << "ties=("; - for (int s = 1; s < nspec - 1; s++) { - ss << "f"; - ss << boost::lexical_cast<std::string>(s); - ss << ".w=f0.w,"; + ss << "A=0.5,"; + ss << "w=" << freq << ","; + ss << "p=0.5;"; + if (fixFreq) { + // w is shared across workspaces + ss << "ties=(f0.w=" << freq << ")"; } - ss << "f"; - ss << boost::lexical_cast<std::string>(nspec - 1); - ss << ".w=f0.w)"; return ss.str(); } -/** Extracts relevant data from a workspace and removes exponential decay -* @param ws :: [input] The input workspace +/** Extracts relevant data from a workspace * @param startTime :: [input] First X value to consider * @param endTime :: [input] Last X value to consider * @return :: Pre-processed workspace to fit */ API::MatrixWorkspace_sptr -CalMuonDetectorPhases::prepareWorkspace(const API::MatrixWorkspace_sptr &ws, - double startTime, double endTime) { - +CalMuonDetectorPhases::extractDataFromWorkspace(double startTime, + double endTime) { // Extract counts from startTime to endTime API::IAlgorithm_sptr crop = createChildAlgorithm("CropWorkspace"); - crop->setProperty("InputWorkspace", ws); + crop->setProperty("InputWorkspace", m_inputWS); crop->setProperty("XMin", startTime); crop->setProperty("XMax", endTime); crop->executeAsChildAlg(); boost::shared_ptr<API::MatrixWorkspace> wsCrop = crop->getProperty("OutputWorkspace"); + return wsCrop; +} - // Remove exponential decay +/** + * Removes exponential decay from a workspace + * @param wsInput :: [input] Workspace to work on + * @return :: Workspace with decay removed + */ +API::MatrixWorkspace_sptr CalMuonDetectorPhases::removeExpDecay( + const API::MatrixWorkspace_sptr &wsInput) { API::IAlgorithm_sptr remove = createChildAlgorithm("RemoveExpDecay"); - remove->setProperty("InputWorkspace", wsCrop); + remove->setProperty("InputWorkspace", wsInput); remove->executeAsChildAlg(); - boost::shared_ptr<API::MatrixWorkspace> wsRem = - remove->getProperty("OutputWorkspace"); - + API::MatrixWorkspace_sptr wsRem = remove->getProperty("OutputWorkspace"); return wsRem; } +/** + * Returns the frequency hint to use as a starting point for finding the + * frequency. + * + * If user has provided a frequency as input, use that. + * Otherwise, use 2*pi*g_mu*(sample_magn_field) + * + * @return :: Frequency hint to use + */ +double CalMuonDetectorPhases::getFrequencyHint() const { + double freq = getProperty("Frequency"); + + if (freq == EMPTY_DBL()) { + try { + // Read sample_magn_field from workspace logs + freq = m_inputWS->run().getLogAsSingleValue("sample_magn_field"); + // Multiply by muon gyromagnetic ratio: 0.01355 MHz/G + freq *= 2 * M_PI * PhysicalConstants::MuonGyromagneticRatio; + } catch (...) { + throw std::runtime_error( + "Couldn't read sample_magn_field. Please provide a value for " + "the frequency"); + } + } + return freq; +} + +/** + * Returns the frequency to use in the sequential fit. + * + * Finds this by grouping the spectra and calculating the asymmetry, then + * fitting this to get the frequency. + * The starting value for this fit is taken from the frequency hint or logs. + * @param ws :: [input] Pointer to cropped workspace with exp decay removed + * @return :: Fixed frequency value to use in the sequential fit + */ +double +CalMuonDetectorPhases::getFrequency(const API::MatrixWorkspace_sptr &ws) { + std::vector<int> forward = getProperty("ForwardSpectra"); + std::vector<int> backward = getProperty("BackwardSpectra"); + + // If grouping not provided, read it from the instrument + if (forward.empty() || backward.empty()) { + getGroupingFromInstrument(ws, forward, backward); + } + + // Calculate asymmetry + const double alpha = getAlpha(ws, forward, backward); + const API::MatrixWorkspace_sptr wsAsym = + getAsymmetry(ws, forward, backward, alpha); + + // Fit an oscillating function, allowing frequency to vary + double frequency = fitFrequencyFromAsymmetry(wsAsym); + + return frequency; +} + +/** + * If grouping was not provided, find the instrument from the input workspace + * and read the default grouping from its IDF. Returns the forward and backward + * groupings as arrays of integers. + * @param ws :: [input] Workspace to find grouping for + * @param forward :: [output] Forward spectrum indices for given instrument + * @param backward :: [output] Backward spectrum indices for given instrument + */ +void CalMuonDetectorPhases::getGroupingFromInstrument( + const API::MatrixWorkspace_sptr &ws, std::vector<int> &forward, + std::vector<int> &backward) { + // make sure both arrays are empty + forward.clear(); + backward.clear(); + + const auto instrument = ws->getInstrument(); + if (instrument->getName() == "MUSR") { + // Two possibilities for grouping, but we have no way of knowing which + throw new std::invalid_argument( + "Cannot use default instrument grouping for MUSR " + "as main field direction is unknown"); + } + + // Load grouping and find forward, backward groups + std::string fwdRange, bwdRange; + API::GroupingLoader loader(instrument); + const auto grouping = loader.getGroupingFromIDF(); + size_t nGroups = grouping->groups.size(); + for (size_t iGroup = 0; iGroup < nGroups; iGroup++) { + const std::string name = grouping->groupNames[iGroup]; + if (name == "fwd") { + fwdRange = grouping->groups[iGroup]; + } else if (name == "bwd" || name == "bkwd") { + bwdRange = grouping->groups[iGroup]; + } + } + + // Use ArrayProperty's functionality to convert string ranges to groups + this->setProperty("ForwardSpectra", fwdRange); + this->setProperty("BackwardSpectra", bwdRange); + forward = getProperty("ForwardSpectra"); + backward = getProperty("BackwardSpectra"); +} + +/** + * Get start time for fit + * If not provided as input, try to read from workspace logs. + * If it's not there either, set to 0 and warn user. + * @return :: Start time for fit + */ +double CalMuonDetectorPhases::getStartTime() const { + double startTime = getProperty("FirstGoodData"); + if (startTime == EMPTY_DBL()) { + try { + // Read FirstGoodData from workspace logs if possible + double firstGoodData = + m_inputWS->run().getLogAsSingleValue("FirstGoodData"); + startTime = firstGoodData; + } catch (...) { + g_log.warning("Couldn't read FirstGoodData, setting to 0"); + startTime = 0.; + } + } + return startTime; +} + +/** + * Get end time for fit + * If it's not there, use the last available time in the workspace. + * @return :: End time for fit + */ +double CalMuonDetectorPhases::getEndTime() const { + double endTime = getProperty("LastGoodData"); + if (endTime == EMPTY_DBL()) { + // Last available time + endTime = m_inputWS->readX(0).back(); + } + return endTime; +} + +/** + * Calculate alpha (detector efficiency) from the given workspace + * If calculation fails, returns default 1.0 + * @param ws :: [input] Workspace to calculate alpha from + * @param forward :: [input] Forward group spectra numbers + * @param backward :: [input] Backward group spectra numbers + * @return :: Alpha, or 1.0 if calculation failed + */ +double CalMuonDetectorPhases::getAlpha(const API::MatrixWorkspace_sptr &ws, + const std::vector<int> &forward, + const std::vector<int> &backward) { + double alpha = 1.0; + try { + auto alphaAlg = createChildAlgorithm("AlphaCalc"); + alphaAlg->setProperty("InputWorkspace", ws); + alphaAlg->setProperty("ForwardSpectra", forward); + alphaAlg->setProperty("BackwardSpectra", backward); + alphaAlg->executeAsChildAlg(); + alpha = alphaAlg->getProperty("Alpha"); + } catch (const std::exception &e) { + // Eat the error and return default 1.0 so algorithm can continue. + // Warn the user that calculating alpha failed + std::ostringstream message; + message << "Calculating alpha failed, default to 1.0: " << e.what(); + g_log.error(message.str()); + } + return alpha; +} + +/** + * Calculate asymmetry for the given workspace + * @param ws :: [input] Workspace to calculate asymmetry from + * @param forward :: [input] Forward group spectra numbers + * @param backward :: [input] Backward group spectra numbers + * @param alpha :: [input] Detector efficiency + * @return :: Asymmetry for workspace + */ +API::MatrixWorkspace_sptr CalMuonDetectorPhases::getAsymmetry( + const API::MatrixWorkspace_sptr &ws, const std::vector<int> &forward, + const std::vector<int> &backward, const double alpha) { + auto alg = createChildAlgorithm("AsymmetryCalc"); + alg->setProperty("InputWorkspace", ws); + alg->setProperty("OutputWorkspace", "__NotUsed"); + alg->setProperty("ForwardSpectra", forward); + alg->setProperty("BackwardSpectra", backward); + alg->setProperty("Alpha", alpha); + alg->executeAsChildAlg(); + API::MatrixWorkspace_sptr wsAsym = alg->getProperty("OutputWorkspace"); + return wsAsym; +} + +/** + * Fit the asymmetry and return the frequency found. + * Starting value for the frequency is taken from the hint. + * If the fit fails, return the initial hint. + * @param wsAsym :: [input] Workspace with asymmetry to fit + * @return :: Frequency found from fit + */ +double CalMuonDetectorPhases::fitFrequencyFromAsymmetry( + const API::MatrixWorkspace_sptr &wsAsym) { + // Starting value for frequency is hint + double hint = getFrequencyHint(); + std::string funcStr = createFittingFunction(hint, false); + double frequency = hint; + + std::string fitStatus = "success"; + try { + auto func = API::FunctionFactory::Instance().createInitialized(funcStr); + auto fit = createChildAlgorithm("Fit"); + fit->setProperty("Function", func); + fit->setProperty("InputWorkspace", wsAsym); + fit->setProperty("WorkspaceIndex", 0); + fit->setProperty("CreateOutput", true); + fit->setProperty("OutputParametersOnly", true); + fit->setProperty("Output", "__Invisible"); + fit->executeAsChildAlg(); + fitStatus = fit->getPropertyValue("OutputStatus"); + if (fitStatus == "success") { + API::ITableWorkspace_sptr params = fit->getProperty("OutputParameters"); + const size_t rows = params->rowCount(); + static size_t colName(0), colValue(1); + for (size_t iRow = 0; iRow < rows; iRow++) { + if (params->cell<std::string>(iRow, colName) == "w") { + frequency = params->cell<double>(iRow, colValue); + break; + } + } + } + } catch (const std::exception &e) { + // Report fit failure to user + fitStatus = e.what(); + } + if (fitStatus != "success") { // Either failed, or threw an exception + std::ostringstream message; + message << "Fit failed (" << fitStatus << "), using omega hint = " << hint; + g_log.error(message.str()); + } + return frequency; +} + +/** + * Updates the algorithm progress + * @param thisSpectrum :: [input] Spectrum number currently being fitted + * @param totalSpectra :: [input] Total number of spectra to fit + */ +void CalMuonDetectorPhases::reportProgress(const int thisSpectrum, + const int totalSpectra) { + double proportionDone = (double)thisSpectrum / (double)totalSpectra; + std::ostringstream progMessage; + progMessage << "Fitting " << thisSpectrum + 1 << " of " << totalSpectra; + this->progress(proportionDone, progMessage.str()); +} + } // namespace Algorithms } // namespace Mantid diff --git a/Framework/Algorithms/test/CalMuonDetectorPhasesTest.h b/Framework/Algorithms/test/CalMuonDetectorPhasesTest.h index d8184ef6afa689ebcd4f44f1fd7b92f672aa6c48..bfe88c7c86739e497f7aa98d9405afa923897de2 100644 --- a/Framework/Algorithms/test/CalMuonDetectorPhasesTest.h +++ b/Framework/Algorithms/test/CalMuonDetectorPhasesTest.h @@ -43,6 +43,8 @@ public: calc->setPropertyValue("Frequency", "25"); calc->setPropertyValue("DataFitted", "fit"); calc->setPropertyValue("DetectorTable", "tab"); + calc->setProperty("ForwardSpectra", std::vector<int>{1, 2}); + calc->setProperty("BackwardSpectra", std::vector<int>{3, 4}); TS_ASSERT_THROWS_NOTHING(calc->execute()); @@ -54,19 +56,19 @@ public: TS_ASSERT_EQUALS(tab->columnCount(), 3); // Test asymmetries TS_ASSERT_DELTA(tab->Double(0, 1), 0.099, 0.001); - TS_ASSERT_DELTA(tab->Double(1, 1), 0.099, 0.001); - TS_ASSERT_DELTA(tab->Double(2, 1), 0.099, 0.001); + TS_ASSERT_DELTA(tab->Double(1, 1), 0.100, 0.001); + TS_ASSERT_DELTA(tab->Double(2, 1), 0.100, 0.001); TS_ASSERT_DELTA(tab->Double(3, 1), 0.100, 0.001); // Test phases - TS_ASSERT_DELTA(tab->Double(0, 2), 6.281, 0.001); - TS_ASSERT_DELTA(tab->Double(1, 2), 0.785, 0.001); - TS_ASSERT_DELTA(tab->Double(2, 2), 1.570, 0.001); - TS_ASSERT_DELTA(tab->Double(3, 2), 2.354, 0.001); + TS_ASSERT_DELTA(tab->Double(0, 2), 6.278, 0.001); + TS_ASSERT_DELTA(tab->Double(1, 2), 0.781, 0.001); + TS_ASSERT_DELTA(tab->Double(2, 2), 1.566, 0.001); + TS_ASSERT_DELTA(tab->Double(3, 2), 2.350, 0.001); } void testBadWorkspaceUnits() { - auto ws = createWorkspace(1, 4, "Wavelength"); + auto ws = createWorkspace(2, 4, "Wavelength"); auto calc = AlgorithmManager::Instance().create("CalMuonDetectorPhases"); calc->initialize(); calc->setChild(true); @@ -74,6 +76,8 @@ public: calc->setPropertyValue("Frequency", "25"); calc->setPropertyValue("DataFitted", "fit"); calc->setPropertyValue("DetectorTable", "tab"); + calc->setProperty("ForwardSpectra", std::vector<int>{1}); + calc->setProperty("BackwardSpectra", std::vector<int>{2}); TS_ASSERT_THROWS(calc->execute(), std::runtime_error); TS_ASSERT(!calc->isExecuted()); @@ -81,13 +85,15 @@ public: void testNoFrequencySupplied() { - auto ws = createWorkspace(1, 4, "Microseconds"); + auto ws = createWorkspace(2, 4, "Microseconds"); auto calc = AlgorithmManager::Instance().create("CalMuonDetectorPhases"); calc->initialize(); calc->setChild(true); calc->setProperty("InputWorkspace", ws); calc->setPropertyValue("DataFitted", "fit"); calc->setPropertyValue("DetectorTable", "tab"); + calc->setProperty("ForwardSpectra", std::vector<int>{1}); + calc->setProperty("BackwardSpectra", std::vector<int>{2}); TS_ASSERT_THROWS(calc->execute(), std::runtime_error); TS_ASSERT(!calc->isExecuted()); diff --git a/Framework/CurveFitting/src/Functions/DynamicKuboToyabe.cpp b/Framework/CurveFitting/src/Functions/DynamicKuboToyabe.cpp index e4c5fff1e176627241c6c75c4867d7311d801184..514424fc0f21c449a68658dc3afd1d0c201f3918 100644 --- a/Framework/CurveFitting/src/Functions/DynamicKuboToyabe.cpp +++ b/Framework/CurveFitting/src/Functions/DynamicKuboToyabe.cpp @@ -5,6 +5,7 @@ #include "MantidAPI/Jacobian.h" #include "MantidAPI/FunctionFactory.h" #include "MantidAPI/MatrixWorkspace.h" +#include "MantidKernel/PhysicalConstants.h" #include <vector> namespace Mantid { @@ -141,7 +142,8 @@ double ZFKT(const double x, const double G) { double HKT(const double x, const double G, const double F) { const double q = G * G * x * x; - const double gm = 2 * M_PI * 0.01355342; // Muon gyromagnetic ratio * 2 * PI + // Muon gyromagnetic ratio * 2 * PI + const double gm = 2 * M_PI * PhysicalConstants::MuonGyromagneticRatio; double w; if (F > 2 * G) { diff --git a/Framework/Kernel/inc/MantidKernel/PhysicalConstants.h b/Framework/Kernel/inc/MantidKernel/PhysicalConstants.h index 3ec37fc890bc097a63fdce12b6403a650fae609d..c384229be657f2d1b4f26a801f8f853f5391f619 100644 --- a/Framework/Kernel/inc/MantidKernel/PhysicalConstants.h +++ b/Framework/Kernel/inc/MantidKernel/PhysicalConstants.h @@ -94,6 +94,11 @@ static const double StandardAtmosphere = 101.325; */ static const double BoltzmannConstant = 8.6173324e-02; +/** Muon gyromagnetic ratio in MHz/G + * Taken from CalMuonDetectorPhases and DynamicKuboToyabe on 02/02/2016 + */ +static const double MuonGyromagneticRatio = 0.01355342; + } // namespace PhysicalConstants } // namespace Mantid diff --git a/MantidQt/CustomInterfaces/inc/MantidQtCustomInterfaces/Muon/IO_MuonGrouping.h b/MantidQt/CustomInterfaces/inc/MantidQtCustomInterfaces/Muon/IO_MuonGrouping.h index 4c83b29c945c066dc6cd36224741b204eecdad7d..1249bde21142ba1e580e8fd21d322059e52ecbf7 100644 --- a/MantidQt/CustomInterfaces/inc/MantidQtCustomInterfaces/Muon/IO_MuonGrouping.h +++ b/MantidQt/CustomInterfaces/inc/MantidQtCustomInterfaces/Muon/IO_MuonGrouping.h @@ -45,38 +45,21 @@ Code Documentation is available at: <http://doxygen.mantidproject.org> using namespace Mantid; using namespace Mantid::API; -/// Structure to represent grouping information for Muon Analysis -struct Grouping { - std::vector<std::string> groupNames; - std::vector<std::string> groups; // Range strings, e.g. "1-32" - - std::vector<std::string> pairNames; - std::vector<std::pair<size_t, size_t> > pairs; // Pairs of group ids - std::vector<double> pairAlphas; - - std::string description; - std::string defaultName; // Not storing id because can be either group or pair -}; - /// Saves grouping to the XML file specified -void MANTIDQT_CUSTOMINTERFACES_DLL saveGroupingToXML(const Grouping& grouping, - const std::string& filename); - -/// Loads grouping from the XML file specified -void MANTIDQT_CUSTOMINTERFACES_DLL loadGroupingFromXML(const std::string& filename, - Grouping& grouping); +void MANTIDQT_CUSTOMINTERFACES_DLL saveGroupingToXML( + const Mantid::API::Grouping &grouping, const std::string &filename); /// Parses information from the grouping table and saves to Grouping struct -void MANTIDQT_CUSTOMINTERFACES_DLL parseGroupingTable(const Ui::MuonAnalysis& form, - Grouping& grouping); +void MANTIDQT_CUSTOMINTERFACES_DLL parseGroupingTable( + const Ui::MuonAnalysis &form, Mantid::API::Grouping &grouping); /// Fills in the grouping table using information from provided Grouping struct -void MANTIDQT_CUSTOMINTERFACES_DLL fillGroupingTable(const Grouping& grouping, - Ui::MuonAnalysis& form); +void MANTIDQT_CUSTOMINTERFACES_DLL fillGroupingTable( + const Mantid::API::Grouping &grouping, Ui::MuonAnalysis &form); /// Groups the workspace according to grouping provided -MatrixWorkspace_sptr MANTIDQT_CUSTOMINTERFACES_DLL groupWorkspace(MatrixWorkspace_const_sptr ws, - const Grouping& g); +MatrixWorkspace_sptr MANTIDQT_CUSTOMINTERFACES_DLL +groupWorkspace(MatrixWorkspace_const_sptr ws, const Mantid::API::Grouping &g); /// create 'map' relating group number to row number in group table void MANTIDQT_CUSTOMINTERFACES_DLL whichGroupToWhichRow(const Ui::MuonAnalysis& m_uiForm, @@ -91,19 +74,16 @@ void MANTIDQT_CUSTOMINTERFACES_DLL setGroupGroupPair(Ui::MuonAnalysis& m_uiForm, const std::string& name); /// Convert a grouping table to a grouping struct -boost::shared_ptr<Grouping> MANTIDQT_CUSTOMINTERFACES_DLL tableToGrouping(ITableWorkspace_sptr table); +boost::shared_ptr<Mantid::API::Grouping> + MANTIDQT_CUSTOMINTERFACES_DLL tableToGrouping(ITableWorkspace_sptr table); /// Converts a grouping information to a grouping table -ITableWorkspace_sptr MANTIDQT_CUSTOMINTERFACES_DLL groupingToTable(boost::shared_ptr<Grouping> grouping); +ITableWorkspace_sptr MANTIDQT_CUSTOMINTERFACES_DLL +groupingToTable(boost::shared_ptr<Mantid::API::Grouping> grouping); /// Returns a "dummy" grouping which a single group with all the detectors in it -boost::shared_ptr<Grouping> getDummyGrouping(Instrument_const_sptr instrument); - -/// Attempts to load a grouping information referenced by IDF -boost::shared_ptr<Grouping> getGroupingFromIDF(Instrument_const_sptr instrument, - const std::string& mainFieldDirection); - - +boost::shared_ptr<Mantid::API::Grouping> +getDummyGrouping(Instrument_const_sptr instrument); } } } diff --git a/MantidQt/CustomInterfaces/inc/MantidQtCustomInterfaces/Muon/MuonAnalysis.h b/MantidQt/CustomInterfaces/inc/MantidQtCustomInterfaces/Muon/MuonAnalysis.h index 33abadbed6d0c2b14585f3ac783e09dd39f9b587..c7f3377920b85d733f5224d1f26bbccdf544ee6e 100644 --- a/MantidQt/CustomInterfaces/inc/MantidQtCustomInterfaces/Muon/MuonAnalysis.h +++ b/MantidQt/CustomInterfaces/inc/MantidQtCustomInterfaces/Muon/MuonAnalysis.h @@ -9,6 +9,7 @@ #include "MantidAPI/AnalysisDataService.h" #include "MantidAPI/MatrixWorkspace_fwd.h" #include "MantidAPI/ITableWorkspace_fwd.h" +#include "MantidAPI/GroupingLoader.h" #include "MantidGeometry/Instrument.h" @@ -34,7 +35,6 @@ namespace Muon class MuonAnalysisOptionTab; class MuonAnalysisFitDataTab; class MuonAnalysisResultTableTab; - struct Grouping; struct LoadResult { Workspace_sptr loadedWorkspace; @@ -48,7 +48,7 @@ namespace Muon struct GroupResult { bool usedExistGrouping; - boost::shared_ptr<Grouping> groupingUsed; + boost::shared_ptr<Mantid::API::Grouping> groupingUsed; }; } diff --git a/MantidQt/CustomInterfaces/src/Muon/IO_MuonGrouping.cpp b/MantidQt/CustomInterfaces/src/Muon/IO_MuonGrouping.cpp index 422feabe5fb8a858d4d2f62eef8dcc5b85519ff7..d888b4d5fe2dbba068b10e3b089f84f3effb128c 100644 --- a/MantidQt/CustomInterfaces/src/Muon/IO_MuonGrouping.cpp +++ b/MantidQt/CustomInterfaces/src/Muon/IO_MuonGrouping.cpp @@ -49,8 +49,8 @@ using namespace MantidQt::API; * @param g :: Struct with grouping information * @param filename :: XML filename where information will be saved */ -void saveGroupingToXML(const Grouping& g, const std::string& filename) -{ +void saveGroupingToXML(const Mantid::API::Grouping &g, + const std::string &filename) { std::ofstream outFile(filename.c_str()); if (!outFile) throw Mantid::Kernel::Exception::FileError("Unable to open output file", filename); @@ -106,165 +106,14 @@ void saveGroupingToXML(const Grouping& g, const std::string& filename) writer.writeNode(outFile, mDoc); } -/** - * Loads grouping from the XML file specified. - * - * @param filename :: XML filename to load grouping information from - * @param g :: Struct to store grouping information to - */ -void loadGroupingFromXML(const std::string& filename, Grouping& g) -{ - // Set up the DOM parser and parse xml file - DOMParser pParser; - Poco::AutoPtr<Document> pDoc; - try - { - pDoc = pParser.parse(filename); - } - catch(...) - { - throw Mantid::Kernel::Exception::FileError("Unable to parse File" , filename); - } - - // Get pointer to root element - Element* pRootElem = pDoc->documentElement(); - if (!pRootElem->hasChildNodes()) - throw Mantid::Kernel::Exception::FileError("No root element in XML grouping file" , filename); - - // Parse information for groups - Poco::AutoPtr<NodeList> groups = pRootElem->getElementsByTagName("group"); - if (groups->length() == 0) - throw Mantid::Kernel::Exception::FileError("No groups specified in XML grouping file" , filename); - - // Resize vectors - g.groupNames.resize(groups->length()); - g.groups.resize(groups->length()); - - for (size_t ig = 0; ig < groups->length(); ig++) - { - Element* pGroupElem = static_cast<Element*>(groups->item(static_cast<long>(ig))); - - if (!pGroupElem->hasAttribute("name")) - throw Mantid::Kernel::Exception::FileError("Group element without name" , filename); - - g.groupNames[ig] = pGroupElem->getAttribute("name"); - - Element* idlistElement = pGroupElem->getChildElement("ids"); - if (!idlistElement) - throw Mantid::Kernel::Exception::FileError("Group element without <ids>" , filename); - - g.groups[ig] = idlistElement->getAttribute("val"); - } - - - // Parse information for pairs - Poco::AutoPtr<NodeList> pairs = pRootElem->getElementsByTagName("pair"); - - // Resize vectors - g.pairNames.resize(pairs->length()); - g.pairs.resize(pairs->length()); - g.pairAlphas.resize(pairs->length()); - - for (size_t ip = 0; ip < pairs->length(); ip++) - { - Element* pPairElem = static_cast<Element*>(pairs->item(static_cast<long>(ip))); - - if ( !pPairElem->hasAttribute("name") ) - throw Mantid::Kernel::Exception::FileError("Pair element without name" , filename); - - g.pairNames[ip] = pPairElem->getAttribute("name"); - - size_t fwdGroupId, bwdGroupId; // Ids of forward/backward groups - - // Try to get id of the first group - if (Element* fwdElement = pPairElem->getChildElement("forward-group")) - { - if(!fwdElement->hasAttribute("val")) - throw Mantid::Kernel::Exception::FileError("Pair forward-group without <val>" , filename); - - // Find the group with the given name - auto it = std::find(g.groupNames.begin(), g.groupNames.end(), fwdElement->getAttribute("val")); - - if(it == g.groupNames.end()) - throw Mantid::Kernel::Exception::FileError("Pair forward-group name not recognized" , filename); - - // Get index of the iterator - fwdGroupId = it - g.groupNames.begin(); - } - else - { - throw Mantid::Kernel::Exception::FileError("Pair element without <forward-group>" , filename); - } - - // Try to get id of the second group - if(Element* bwdElement = pPairElem->getChildElement("backward-group")) - { - if(!bwdElement->hasAttribute("val")) - throw Mantid::Kernel::Exception::FileError("Pair backward-group without <val>" , filename); - - // Find the group with the given name - auto it = std::find(g.groupNames.begin(), g.groupNames.end(), bwdElement->getAttribute("val")); - - if(it == g.groupNames.end()) - throw Mantid::Kernel::Exception::FileError("Pair backward-group name not recognized" , filename); - - // Get index of the iterator - bwdGroupId = it - g.groupNames.begin(); - } - else - { - throw Mantid::Kernel::Exception::FileError("Pair element without <backward-group>" , filename); - } - - g.pairs[ip] = std::make_pair(fwdGroupId, bwdGroupId); - - // Try to get alpha element - if (Element* aElement = pPairElem->getChildElement("alpha")) - { - if (!aElement->hasAttribute("val") ) - throw Mantid::Kernel::Exception::FileError("Pair alpha element with no <val>" , filename); - - try // ... to convert value to double - { - g.pairAlphas[ip] = boost::lexical_cast<double>(aElement->getAttribute("val")); - } - catch(boost::bad_lexical_cast&) - { - throw Mantid::Kernel::Exception::FileError("Pair alpha value is not a number" , filename); - } - } - // If alpha element not there, default it to 1.0 - else - { - g.pairAlphas[ip] = 1.0; - } - - } - - // Try to get description - if (pRootElem->hasAttribute("description")) - { - g.description = pRootElem->getAttribute("description"); - } - - // Try to get default group/pair name - if(Element* defaultElement = pRootElem->getChildElement("default")) - { - if(!defaultElement->hasAttribute("name")) - throw Mantid::Kernel::Exception::FileError("Default element with no <name>" , filename); - - g.defaultName = defaultElement->getAttribute("name"); - } -} - /** * Parses information from the grouping table and saves to Grouping struct. * * @param form :: Muon Analysis UI containing table widgets * @param g :: Grouping struct to store parsed info to */ -void parseGroupingTable(const Ui::MuonAnalysis& form, Grouping& g) -{ +void parseGroupingTable(const Ui::MuonAnalysis &form, + Mantid::API::Grouping &g) { // Parse description g.description = form.groupDescription->text().toStdString(); @@ -315,7 +164,7 @@ void parseGroupingTable(const Ui::MuonAnalysis& form, Grouping& g) * @param g :: Grouping struct to use for filling the table * @param form :: Muon Analysis UI containing table widgets */ -void fillGroupingTable(const Grouping& g, Ui::MuonAnalysis& form) +void fillGroupingTable(const Mantid::API::Grouping& g, Ui::MuonAnalysis& form) { // Add groups to a table for(int gi = 0; gi < static_cast<int>(g.groups.size()); gi++) @@ -374,8 +223,8 @@ void setGroupGroupPair(Ui::MuonAnalysis& m_uiForm, const std::string& name) * @param g :: The grouping information * @return Sptr to created grouped workspace */ -MatrixWorkspace_sptr groupWorkspace(MatrixWorkspace_const_sptr ws, const Grouping& g) -{ +MatrixWorkspace_sptr groupWorkspace(MatrixWorkspace_const_sptr ws, + const Mantid::API::Grouping &g) { // As I couldn't specify multiple groups for GroupDetectors, I am going down quite a complicated // route - for every group distinct grouped workspace is created using GroupDetectors. These // workspaces are then merged into the output workspace. @@ -495,9 +344,9 @@ void whichPairToWhichRow(const Ui::MuonAnalysis& m_uiForm, std::vector<int>& pai * @param table :: A table to convert * @return Grouping info */ -boost::shared_ptr<Grouping> tableToGrouping(ITableWorkspace_sptr table) -{ - auto grouping = boost::make_shared<Grouping>(); +boost::shared_ptr<Mantid::API::Grouping> +tableToGrouping(ITableWorkspace_sptr table) { + auto grouping = boost::make_shared<Mantid::API::Grouping>(); for ( size_t row = 0; row < table->rowCount(); ++row ) { @@ -530,8 +379,8 @@ boost::shared_ptr<Grouping> tableToGrouping(ITableWorkspace_sptr table) * @param grouping :: Grouping information to convert * @return A grouping table as accepted by MuonGroupDetectors */ -ITableWorkspace_sptr groupingToTable(boost::shared_ptr<Grouping> grouping) -{ +ITableWorkspace_sptr +groupingToTable(boost::shared_ptr<Mantid::API::Grouping> grouping) { auto newTable = boost::dynamic_pointer_cast<ITableWorkspace>( WorkspaceFactory::Instance().createTable("TableWorkspace") ); @@ -551,56 +400,19 @@ ITableWorkspace_sptr groupingToTable(boost::shared_ptr<Grouping> grouping) * @param instrument :: Instrument we want a dummy grouping for * @return Grouping information */ -boost::shared_ptr<Grouping> getDummyGrouping(Instrument_const_sptr instrument) -{ +boost::shared_ptr<Mantid::API::Grouping> +getDummyGrouping(Instrument_const_sptr instrument) { // Group with all the detectors std::ostringstream all; all << "1-" << instrument->getNumberDetectors(); - auto dummyGrouping = boost::make_shared<Grouping>(); + auto dummyGrouping = boost::make_shared<Mantid::API::Grouping>(); dummyGrouping->description = "Dummy grouping"; dummyGrouping->groupNames.emplace_back("all"); dummyGrouping->groups.push_back(all.str()); return dummyGrouping; } -/** - * Attempts to load a grouping information referenced by IDF. - * @param instrument :: Intrument which we went the grouping for - * @param mainFieldDirection :: (MUSR) orientation of the instrument - * @return Grouping information - */ -boost::shared_ptr<Grouping> getGroupingFromIDF(Instrument_const_sptr instrument, - const std::string& mainFieldDirection) -{ - std::string parameterName = "Default grouping file"; - - // Special case for MUSR, because it has two possible groupings - if (instrument->getName() == "MUSR") - { - parameterName.append(" - " + mainFieldDirection); - } - - std::vector<std::string> groupingFiles = instrument->getStringParameter(parameterName); - - if ( groupingFiles.size() == 1 ) - { - const std::string groupingFile = groupingFiles[0]; - - // Get search directory for XML instrument definition files (IDFs) - std::string directoryName = ConfigService::Instance().getInstrumentDirectory(); - - auto loadedGrouping = boost::make_shared<Grouping>(); - loadGroupingFromXML(directoryName + groupingFile, *loadedGrouping); - - return loadedGrouping; - } - else - { - throw std::runtime_error("Multiple groupings specified for the instrument"); - } -} - } } } diff --git a/MantidQt/CustomInterfaces/src/Muon/MuonAnalysis.cpp b/MantidQt/CustomInterfaces/src/Muon/MuonAnalysis.cpp index c31424a27f84d9f780171b5777479f32524a2207..476864aac116dfad108ddc6a91d9a5a7e31853c2 100644 --- a/MantidQt/CustomInterfaces/src/Muon/MuonAnalysis.cpp +++ b/MantidQt/CustomInterfaces/src/Muon/MuonAnalysis.cpp @@ -73,6 +73,7 @@ using namespace MantidQt::CustomInterfaces; using namespace MantidQt::CustomInterfaces::Muon; using namespace Mantid::Geometry; using Mantid::API::Workspace_sptr; +using Mantid::API::Grouping; namespace { @@ -662,7 +663,7 @@ void MuonAnalysis::runSaveGroupButton() if( ! groupingFile.isEmpty() ) { - Grouping groupingToSave; + Mantid::API::Grouping groupingToSave; parseGroupingTable(m_uiForm, groupingToSave); saveGroupingToXML(groupingToSave, groupingFile.toStdString()); @@ -698,11 +699,12 @@ void MuonAnalysis::runLoadGroupButton() QString directory = QFileInfo(groupingFile).path(); prevValues.setValue("dir", directory); - Grouping loadedGrouping; + Mantid::API::Grouping loadedGrouping; try { - loadGroupingFromXML(groupingFile.toStdString(), loadedGrouping); + Mantid::API::GroupingLoader::loadGroupingFromXML(groupingFile.toStdString(), + loadedGrouping); } catch (Exception::FileError& e) { @@ -1353,7 +1355,7 @@ boost::shared_ptr<GroupResult> MuonAnalysis::getGrouping(boost::shared_ptr<LoadResult> loadResult) const { auto result = boost::make_shared<GroupResult>(); - boost::shared_ptr<Grouping> groupingToUse; + boost::shared_ptr<Mantid::API::Grouping> groupingToUse; Instrument_const_sptr instr = firstPeriod(loadResult->loadedWorkspace)->getInstrument(); @@ -1378,7 +1380,7 @@ MuonAnalysis::getGrouping(boost::shared_ptr<LoadResult> loadResult) const { { // Use grouping currently set result->usedExistGrouping = true; - groupingToUse = boost::make_shared<Grouping>(); + groupingToUse = boost::make_shared<Mantid::API::Grouping>(); parseGroupingTable(m_uiForm, *groupingToUse); } else @@ -1388,7 +1390,8 @@ MuonAnalysis::getGrouping(boost::shared_ptr<LoadResult> loadResult) const { try // to get grouping from IDF { - groupingToUse = getGroupingFromIDF(instr, loadResult->mainFieldDirection); + Mantid::API::GroupingLoader loader(instr, loadResult->mainFieldDirection); + groupingToUse = loader.getGroupingFromIDF(); } catch(std::runtime_error& e) { @@ -3019,7 +3022,7 @@ void MuonAnalysis::groupLoadedWorkspace() */ ITableWorkspace_sptr MuonAnalysis::parseGrouping() { - auto grouping = boost::make_shared<Grouping>(); + auto grouping = boost::make_shared<Mantid::API::Grouping>(); parseGroupingTable(m_uiForm, *grouping); return groupingToTable(grouping); } diff --git a/MantidQt/CustomInterfaces/test/IO_MuonGroupingTest.h b/MantidQt/CustomInterfaces/test/IO_MuonGroupingTest.h index 322c75cf05bffa8781ecab3268eb04c91b1a9425..fe299f6600ef8abeec9d57b1eda0cfde53041205 100644 --- a/MantidQt/CustomInterfaces/test/IO_MuonGroupingTest.h +++ b/MantidQt/CustomInterfaces/test/IO_MuonGroupingTest.h @@ -50,34 +50,6 @@ public: FrameworkManager::Instance(); } - void test_loadGroupingFromXML() - { - Grouping g; - - TS_ASSERT_THROWS_NOTHING(loadGroupingFromXML(m_testDataDir + "MUSRGrouping.xml", g)); - - TS_ASSERT_EQUALS(g.groupNames.size(), 2); - TS_ASSERT_EQUALS(g.groupNames[0], "fwd"); - TS_ASSERT_EQUALS(g.groupNames[1], "bwd"); - - TS_ASSERT_EQUALS(g.groups.size(), 2); - TS_ASSERT_EQUALS(g.groups[0], "33-64"); - TS_ASSERT_EQUALS(g.groups[1], "1-32"); - - TS_ASSERT_EQUALS(g.pairNames.size(), 1); - TS_ASSERT_EQUALS(g.pairNames[0], "long"); - - TS_ASSERT_EQUALS(g.pairs.size(), 1); - TS_ASSERT_EQUALS(g.pairs[0].first, 0); - TS_ASSERT_EQUALS(g.pairs[0].second, 1); - - TS_ASSERT_EQUALS(g.pairAlphas.size(), 1); - TS_ASSERT_EQUALS(g.pairAlphas[0], 1); - - TS_ASSERT_EQUALS(g.description, "musr longitudinal (64 detectors)"); - TS_ASSERT_EQUALS(g.defaultName, "long"); - } - void test_saveGroupingToXML() { Grouping g, lg; @@ -85,13 +57,15 @@ public: std::string tmpFile = m_tmpDir + "tmp_MUSRGrouping.xml"; // Load grouping first - TS_ASSERT_THROWS_NOTHING(loadGroupingFromXML(m_testDataDir + "MUSRGrouping.xml", g)); + TS_ASSERT_THROWS_NOTHING(API::GroupingLoader::loadGroupingFromXML( + m_testDataDir + "MUSRGrouping.xml", g)); // Then save it TS_ASSERT_THROWS_NOTHING(saveGroupingToXML(g, tmpFile)); // And load it again - TS_ASSERT_THROWS_NOTHING(loadGroupingFromXML(tmpFile, lg)); + TS_ASSERT_THROWS_NOTHING( + API::GroupingLoader::loadGroupingFromXML(tmpFile, lg)); // Check that all the information was saved TS_ASSERT_EQUALS(lg.groupNames.size(), 2); @@ -123,7 +97,8 @@ public: { // Load grouping for MUSR Grouping g; - TS_ASSERT_THROWS_NOTHING(loadGroupingFromXML(m_testDataDir + "MUSRGrouping.xml", g)); + TS_ASSERT_THROWS_NOTHING(API::GroupingLoader::loadGroupingFromXML( + m_testDataDir + "MUSRGrouping.xml", g)); // Load MUSR data file IAlgorithm_sptr loadAlg = AlgorithmManager::Instance().create("LoadMuonNexus"); diff --git a/docs/source/algorithms/CalMuonDetectorPhases-v1.rst b/docs/source/algorithms/CalMuonDetectorPhases-v1.rst index 96ad027fdc42a3d10ddf9d286661a55928534a7b..a90bf7b0e8041d988774e1bc09f41503ef55283e 100644 --- a/docs/source/algorithms/CalMuonDetectorPhases-v1.rst +++ b/docs/source/algorithms/CalMuonDetectorPhases-v1.rst @@ -16,19 +16,28 @@ the spectra in the input workspace to: .. math:: f_i(t) = A_i \sin\left(\omega t + \phi_i\right) where :math:`\omega` is shared across spectra and :math:`A_i` and :math:`\phi_i` are -detector-dependent. The algorithm outputs a table workspace containing the detector ID (i.e. the +detector-dependent. + +Before the spectra are fitted, :math:`\omega` is determined by grouping the detectors, +calculating the asymmetry and fitting this to get the frequency. This value of :math:`\omega` +is then treated as a fixed constant when fitting the spectra to the function above. + +The algorithm outputs a table workspace containing the detector ID (i.e. the spectrum index), the asymmetry and the phase. This table is intended to be used as the input *PhaseTable* to :ref:`PhaseQuad <algm-PhaseQuad>`. In addition, the fitting results are returned in a workspace group, where each of the items stores the original data (after removing the exponential decay), the data simulated with the fitting function and the difference between data and fit as spectra 0, 1 and 2 respectively. -There are three optional input properties: *FirstGoodData* and *LastGoodData* define the fitting range. +There are five optional input properties: *FirstGoodData* and *LastGoodData* define the fitting range. When left blank, *FirstGoodData* is set to the value stored in the input workspace and *LastGoodData* is set to the last available bin. The optional property *Frequency* allows the user to select an -initial value for :math:`\omega`. If this property is not supplied, the algortihm takes this -value from the *sample_magn_field* log multiplied by :math:`2\pi\cdot g_\mu`, where :math:`g_\mu` is -the muon gyromagnetic ratio (0.01355 MHz/G). +initial value for :math:`\omega` (a starting value for the fit). If this property is not supplied, the +algorithm takes this value from the *sample_magn_field* log multiplied by :math:`2\pi\cdot g_\mu`, where +:math:`g_\mu` is the muon gyromagnetic ratio (0.01355 MHz/G). +Finally, the optional properties *ForwardSpectra* and *BackwardSpectra* are the sets of spectra in the +forward and backward groups. If these are not supplied, the algorithm will find the instrument from the +input workspace and use the default grouping for this instrument. Usage ----- @@ -42,7 +51,7 @@ Usage # Load four spectra from a muon nexus file ws = Load(Filename='MUSR00022725.nxs', SpectrumMin=1, SpectrumMax=4) # Calibrate the phases and amplituds - detectorTable, fittingResults = CalMuonDetectorPhases(InputWorkspace='ws', LastGoodData=4) + detectorTable, fittingResults = CalMuonDetectorPhases(InputWorkspace='ws', LastGoodData=4, ForwardSpectra="1,2", BackwardSpectra="3,4") # Print the result print "Detector 1 has phase %f and amplitude %f" % (detectorTable.cell(0,2), detectorTable.cell(0,1)) @@ -54,10 +63,10 @@ Output: .. testoutput:: CalMuonDetectorPhasesExample - Detector 1 has phase 0.673861 and amplitude 0.133419 - Detector 2 has phase 0.452013 and amplitude 0.134742 - Detector 3 has phase 0.269103 and amplitude 0.149764 - Detector 4 has phase 0.140418 and amplitude 0.153004 + Detector 1 has phase 0.620299 and amplitude 0.133113 + Detector 2 has phase 0.399003 and amplitude 0.134679 + Detector 3 has phase 0.214079 and amplitude 0.149431 + Detector 4 has phase 0.086315 and amplitude 0.152870 .. categories::