diff --git a/Framework/API/CMakeLists.txt b/Framework/API/CMakeLists.txt index eb22fa24b759ea59843fe50858efd05cae7987f0..a717a8b0f1c3b2bc549fc91bc6b380f8ef179ce8 100644 --- a/Framework/API/CMakeLists.txt +++ b/Framework/API/CMakeLists.txt @@ -128,6 +128,7 @@ set ( SRC_FILES src/ScriptBuilder.cpp src/ScriptRepository.cpp src/ScriptRepositoryFactory.cpp + src/SingleCountValidator.cpp src/SpectraAxis.cpp src/SpectraAxisValidator.cpp src/SpectrumDetectorMapping.cpp @@ -321,6 +322,7 @@ set ( INC_FILES inc/MantidAPI/ScriptBuilder.h inc/MantidAPI/ScriptRepository.h inc/MantidAPI/ScriptRepositoryFactory.h + inc/MantidAPI/SingleCountValidator.h inc/MantidAPI/SingleValueParameter.h inc/MantidAPI/SingleValueParameterParser.h inc/MantidAPI/SpectraAxis.h @@ -352,10 +354,10 @@ set ( TEST_FILES AlgorithmFactoryTest.h AlgorithmHasPropertyTest.h AlgorithmHistoryTest.h + AlgorithmMPITest.h AlgorithmManagerTest.h AlgorithmPropertyTest.h AlgorithmProxyTest.h - AlgorithmMPITest.h AlgorithmTest.h AnalysisDataServiceTest.h AsynchronousTest.h @@ -439,6 +441,7 @@ set ( TEST_FILES SampleValidatorTest.h ScopedWorkspaceTest.h ScriptBuilderTest.h + SingleCountValidatorTest.h SpectraAxisTest.h SpectraAxisValidatorTest.h SpectrumDetectorMappingTest.h diff --git a/Framework/API/inc/MantidAPI/SingleCountValidator.h b/Framework/API/inc/MantidAPI/SingleCountValidator.h new file mode 100644 index 0000000000000000000000000000000000000000..d8e57f6cba3d944ee2c46930cadcf56305e0a1f7 --- /dev/null +++ b/Framework/API/inc/MantidAPI/SingleCountValidator.h @@ -0,0 +1,55 @@ +#ifndef MANTID_API_SINGLECOUNTVALIDATOR_H_ +#define MANTID_API_SINGLECOUNTVALIDATOR_H_ + +#include "MantidAPI/MatrixWorkspaceValidator.h" + +namespace Mantid { +namespace API { + +/** SingleCountValidator : This validator checks that there is only a single + entry per spectrum, the counts, so no Time-of-Flight data. Warning: only the + first bin of the workspace is checked, for performance reasons. + + Copyright © 2017 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 SingleCountValidator : public MatrixWorkspaceValidator { +public: + explicit SingleCountValidator(const bool &mustBeSingleCounts = true); + + /// Gets the type of the validator + std::string getType() const { return "single_count"; } + /// Clone the current state + Kernel::IValidator_sptr clone() const override; + +private: + /// Check for validity + std::string checkValidity(const MatrixWorkspace_sptr &ws) const override; + + /// A flag indicating whether this validator requires that the workspace be + /// contain only single counts or not + const bool m_mustBeSingleCount; +}; + +} // namespace API +} // namespace Mantid + +#endif /* MANTID_API_SINGLECOUNTVALIDATOR_H_ */ diff --git a/Framework/API/src/SingleCountValidator.cpp b/Framework/API/src/SingleCountValidator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a0f5a7b4d74ddf666e30793a76356cbaebb28417 --- /dev/null +++ b/Framework/API/src/SingleCountValidator.cpp @@ -0,0 +1,45 @@ +#include "MantidAPI/SingleCountValidator.h" + +namespace Mantid { +namespace API { + +/** Constructor + * + * @param mustBeSingleCount :: Flag indicating whether the check is that a + * workspace should contain single counts only (true, default) or should not + * contain single counts (false). + */ +SingleCountValidator::SingleCountValidator(const bool &mustBeSingleCount) + : MatrixWorkspaceValidator(), m_mustBeSingleCount(mustBeSingleCount) {} + +/// Clone the current state +Kernel::IValidator_sptr SingleCountValidator::clone() const { + return boost::make_shared<SingleCountValidator>(*this); +} + +/** Checks if the workspace contains a single counts when it should not and + * vice-versa. For perofrmance reasons this takes the first spectrum size only, + * instead of relying on MatrixWorkspace::blocksize(). + * @param ws The workspace to validate + * @return A user level description if a problem exists, otherwise an empty + * string + */ +std::string +SingleCountValidator::checkValidity(const MatrixWorkspace_sptr &ws) const { + auto blockSize = ws->histogram(0).size(); + + if (m_mustBeSingleCount) { + if (blockSize == 1) + return ""; + else + return "The workspace must contain single counts for all spectra"; + } else { + if (blockSize != 1) + return ""; + else + return "The workspace must not contain single counts"; + } +} + +} // namespace API +} // namespace Mantid diff --git a/Framework/API/test/SingleCountValidatorTest.h b/Framework/API/test/SingleCountValidatorTest.h new file mode 100644 index 0000000000000000000000000000000000000000..4b3aa4c5ec99417edfe34027e639191bafc2652a --- /dev/null +++ b/Framework/API/test/SingleCountValidatorTest.h @@ -0,0 +1,84 @@ +#ifndef MANTID_API_SINGLECOUNTVALIDATORTEST_H_ +#define MANTID_API_SINGLECOUNTVALIDATORTEST_H_ + +#include <cxxtest/TestSuite.h> + +#include "MantidAPI/SingleCountValidator.h" + +#include "MantidHistogramData/Histogram.h" +#include "MantidTestHelpers/FakeObjects.h" +#include "MantidAPI/WorkspaceFactory.h" + +using namespace Mantid::HistogramData; +using Mantid::API::SingleCountValidator; + +class SingleCountValidatorTest : 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 SingleCountValidatorTest *createSuite() { + return new SingleCountValidatorTest(); + } + static void destroySuite(SingleCountValidatorTest *suite) { delete suite; } + + void test_single_count_workspace_success() { + auto ws = boost::make_shared<WorkspaceTester>(); + ws->initialize(2, 2, 1); + SingleCountValidator validator(true); + TS_ASSERT_EQUALS(validator.isValid(ws), ""); + } + + void test_non_single_count_workspace_failure() { + auto ws = boost::make_shared<WorkspaceTester>(); + ws->initialize(2, 3, 2); + SingleCountValidator validator(true); + TS_ASSERT_EQUALS( + validator.isValid(ws), + "The workspace must contain single counts for all spectra"); + } + + void test_single_count_workspace_failure() { + auto ws = boost::make_shared<WorkspaceTester>(); + ws->initialize(2, 2, 1); + SingleCountValidator validator(false); + TS_ASSERT_EQUALS(validator.isValid(ws), + "The workspace must not contain single counts"); + } + + void test_non_single_count_workspace_success() { + auto ws = boost::make_shared<WorkspaceTester>(); + ws->initialize(2, 3, 2); + SingleCountValidator validator(false); + TS_ASSERT_EQUALS(validator.isValid(ws), ""); + } + + // The next two tests serve as a warning - only the first bin is checked! + void test_variable_bin_workspace_actually_succeeds() { + auto ws = boost::make_shared<VariableBinThrowingTester>(); + ws->initialize(2, 3, 2); + BinEdges bins{-1.0, 1.0}; + Counts counts{1.0}; + Histogram hist(bins, counts); + + ws->setHistogram(0, hist); + + SingleCountValidator validator(true); + TS_ASSERT_EQUALS(validator.isValid(ws), ""); + } + + void test_variable_bin_workspace_actually_fails() { + auto ws = boost::make_shared<VariableBinThrowingTester>(); + ws->initialize(2, 3, 2); + BinEdges bins{-1.0, 1.0}; + Counts counts{1.0}; + Histogram hist(bins, counts); + + ws->setHistogram(0, hist); + + SingleCountValidator validator(false); + TS_ASSERT_EQUALS(validator.isValid(ws), + "The workspace must not contain single counts"); + } +}; + +#endif /* MANTID_API_SINGLECOUNTVALIDATORTEST_H_ */ diff --git a/Framework/Algorithms/inc/MantidAlgorithms/NormaliseToMonitor.h b/Framework/Algorithms/inc/MantidAlgorithms/NormaliseToMonitor.h index 1a378e320c9d833aa625ba396f69d02ec03dccbd..0d62e5d949d7d81b28b2dcc154dbafa34ad152bf 100644 --- a/Framework/Algorithms/inc/MantidAlgorithms/NormaliseToMonitor.h +++ b/Framework/Algorithms/inc/MantidAlgorithms/NormaliseToMonitor.h @@ -102,19 +102,22 @@ private: protected: // for testing void checkProperties(const API::MatrixWorkspace_sptr &inputWorkspace); API::MatrixWorkspace_sptr - getInWSMonitorSpectrum(const API::MatrixWorkspace_sptr &inputWorkspace, - int &spectra_num); + getInWSMonitorSpectrum(const API::MatrixWorkspace_sptr &inputWorkspace); + size_t getInWSMonitorIndex(const API::MatrixWorkspace_sptr &inputWorkspace); API::MatrixWorkspace_sptr - getMonitorWorkspace(const API::MatrixWorkspace_sptr &inputWorkspace, - int &wsID); + getMonitorWorkspace(const API::MatrixWorkspace_sptr &inputWorkspace); API::MatrixWorkspace_sptr - extractMonitorSpectrum(const API::MatrixWorkspace_sptr &WS, - std::size_t index); - bool setIntegrationProps(); + extractMonitorSpectra(const API::MatrixWorkspace_sptr &ws, + const std::vector<size_t> &workspaceIndexes); + bool setIntegrationProps(const bool isSingleCountWorkspace); void normaliseByIntegratedCount(const API::MatrixWorkspace_sptr &inputWorkspace, - API::MatrixWorkspace_sptr &outputWorkspace); + API::MatrixWorkspace_sptr &outputWorkspace, + const bool isSingleCountWorkspace); + + void performHistogramDivision(const API::MatrixWorkspace_sptr &inputWorkspace, + API::MatrixWorkspace_sptr &outputWorkspace); void normaliseBinByBin(const API::MatrixWorkspace_sptr &inputWorkspace, API::MatrixWorkspace_sptr &outputWorkspace); @@ -131,6 +134,8 @@ private: double m_integrationMin = EMPTY_DBL(); /// The upper bound of the integration range double m_integrationMax = EMPTY_DBL(); + bool m_syncScanInput; + std::vector<size_t> m_workspaceIndexes; }; // the internal class to verify and modify interconnected properties affecting diff --git a/Framework/Algorithms/src/NormaliseToMonitor.cpp b/Framework/Algorithms/src/NormaliseToMonitor.cpp index af549f5054909b6c524e19fea0e9c389129c2e1d..9fa0f3816e43bf9f40a1723db4efb95292414308 100644 --- a/Framework/Algorithms/src/NormaliseToMonitor.cpp +++ b/Framework/Algorithms/src/NormaliseToMonitor.cpp @@ -1,21 +1,25 @@ #include "MantidAlgorithms/NormaliseToMonitor.h" #include "MantidAPI/HistogramValidator.h" #include "MantidAPI/RawCountValidator.h" +#include "MantidAPI/SingleCountValidator.h" #include "MantidAPI/SpectraAxis.h" #include "MantidAPI/SpectrumInfo.h" #include "MantidAPI/WorkspaceFactory.h" #include "MantidAPI/WorkspaceOpOverloads.h" #include "MantidDataObjects/EventWorkspace.h" #include "MantidGeometry/Instrument.h" +#include "MantidGeometry/Instrument/DetectorInfo.h" #include "MantidKernel/BoundedValidator.h" #include "MantidKernel/CompositeValidator.h" #include "MantidKernel/EnabledWhenProperty.h" #include "MantidKernel/ListValidator.h" #include "MantidKernel/VectorHelper.h" +#include "MantidTypes/SpectrumDefinition.h" #include <cfloat> #include <numeric> +using namespace Mantid::API; using namespace Mantid::DataObjects; using namespace Mantid::HistogramData; using Mantid::Kernel::IPropertyManager; @@ -40,7 +44,7 @@ bool MonIDPropChanger::isEnabled(const IPropertyManager *algo) const { // is there the ws property, which describes monitors ws. It also disables the // monitors ID property - API::MatrixWorkspace_const_sptr monitorsWS = + MatrixWorkspace_const_sptr monitorsWS = algo->getProperty(MonitorWorkspaceProp); if (monitorsWS) { is_enabled = false; @@ -57,9 +61,8 @@ bool MonIDPropChanger::isConditionChanged(const IPropertyManager *algo) const { if (!is_enabled) return false; // read monitors list from the input workspace - API::MatrixWorkspace_const_sptr inputWS = algo->getProperty(hostWSname); - bool monitors_changed = monitorIdReader(inputWS); - return monitors_changed; + MatrixWorkspace_const_sptr inputWS = algo->getProperty(hostWSname); + return monitorIdReader(inputWS); } // function which modifies the allowed values for the list of monitors. @@ -73,7 +76,7 @@ void MonIDPropChanger::applyChanges(const IPropertyManager *algo, } if (iExistingAllowedValues.empty()) { - API::MatrixWorkspace_const_sptr inputWS = algo->getProperty(hostWSname); + MatrixWorkspace_const_sptr inputWS = algo->getProperty(hostWSname); int spectra_max(-1); if (inputWS) { // let's assume that detectors IDs correspond to spectraID -- @@ -91,7 +94,7 @@ void MonIDPropChanger::applyChanges(const IPropertyManager *algo, // read the monitors list from the workspace and try to do it once for any // particular ws; bool MonIDPropChanger::monitorIdReader( - API::MatrixWorkspace_const_sptr inputWS) const { + MatrixWorkspace_const_sptr inputWS) const { // no workspace if (!inputWS) return false; @@ -121,33 +124,39 @@ bool MonIDPropChanger::monitorIdReader( return true; } } - // index list can be less or equal to the mon list size (some monitors do not - // have spectra) + // index list can be less or equal to the mon list size (some monitors do + // not have spectra) size_t mon_count = (mon.size() < indexList.size()) ? mon.size() : indexList.size(); - std::vector<int> allowed_values(mon_count); - for (size_t i = 0; i < mon_count; i++) { - allowed_values[i] = mon[i]; - } + mon.resize(mon_count); // are known values the same as the values we have just identified? - if (iExistingAllowedValues.size() != mon_count) { + if (iExistingAllowedValues.size() != mon.size()) { iExistingAllowedValues.clear(); - iExistingAllowedValues.assign(allowed_values.begin(), allowed_values.end()); + iExistingAllowedValues.assign(mon.begin(), mon.end()); return true; } // the monitor list has the same size as before. Is it equivalent to the // existing one? bool values_redefined = false; - for (size_t i = 0; i < mon_count; i++) { - if (iExistingAllowedValues[i] != allowed_values[i]) { + for (size_t i = 0; i < mon.size(); i++) { + if (iExistingAllowedValues[i] != mon[i]) { values_redefined = true; - iExistingAllowedValues[i] = allowed_values[i]; + iExistingAllowedValues[i] = mon[i]; } } return values_redefined; } +bool spectrumDefinitionsMatchTimeIndex(const SpectrumDefinition &specDef, + const size_t timeIndex) { + + for (const auto &spec : specDef) + if (spec.second != timeIndex) + return false; + return true; +} + // Register with the algorithm factory DECLARE_ALGORITHM(NormaliseToMonitor) @@ -156,14 +165,19 @@ using namespace API; using std::size_t; void NormaliseToMonitor::init() { - auto val = boost::make_shared<CompositeValidator>(); - val->add<HistogramValidator>(); - val->add<RawCountValidator>(); - // It's been said that we should restrict the unit to being wavelength, but - // I'm not sure about that... + // Must be histograms OR one count per bin + // Must be raw counts + auto validatorHistSingle = + boost::make_shared<CompositeValidator>(CompositeRelation::OR); + validatorHistSingle->add<HistogramValidator>(); + validatorHistSingle->add<SingleCountValidator>(); + auto validator = boost::make_shared<CompositeValidator>(); + validator->add(validatorHistSingle); + validator->add<RawCountValidator>(); + declareProperty( make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, - val), + validator), "Name of the input workspace. Must be a non-distribution histogram."); declareProperty(make_unique<WorkspaceProperty<>>("OutputWorkspace", "", @@ -194,11 +208,11 @@ void NormaliseToMonitor::init() { "MonitorWorkspace")); // ...or provide it in a separate workspace (note: optional WorkspaceProperty) - declareProperty(make_unique<WorkspaceProperty<>>("MonitorWorkspace", "", - Direction::Input, - PropertyMode::Optional, val), - "A workspace containing one or more spectra to normalize the " - "InputWorkspace by."); + declareProperty( + make_unique<WorkspaceProperty<>>("MonitorWorkspace", "", Direction::Input, + PropertyMode::Optional, validator), + "A workspace containing one or more spectra to normalize the " + "InputWorkspace by."); setPropertySettings("MonitorWorkspace", Kernel::make_unique<Kernel::EnabledWhenProperty>( "MonitorSpectrum", IS_DEFAULT)); @@ -245,16 +259,24 @@ void NormaliseToMonitor::exec() { const MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace"); MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace"); // First check the inputs - this->checkProperties(inputWS); + checkProperties(inputWS); + + bool isSingleCountWorkspace = false; + try { + isSingleCountWorkspace = + (!inputWS->isHistogramData()) && (inputWS->blocksize() == 1); + } catch (std::length_error &) { + // inconsistent bin size, not a single count workspace + isSingleCountWorkspace = false; + } // See if the normalization with integration properties are set. - const bool integrate = this->setIntegrationProps(); + const bool integrate = setIntegrationProps(isSingleCountWorkspace); - if (integrate) { - this->normaliseByIntegratedCount(inputWS, outputWS); - } else { - this->normaliseBinByBin(inputWS, outputWS); - } + if (integrate) + normaliseByIntegratedCount(inputWS, outputWS, isSingleCountWorkspace); + else + normaliseBinByBin(inputWS, outputWS); setProperty("OutputWorkspace", outputWS); std::string norm_ws_name = getPropertyValue("NormFactorWS"); @@ -267,10 +289,29 @@ void NormaliseToMonitor::exec() { auto pProp = (this->getPointerToProperty("NormFactorWS")); pProp->setValue(norm_ws_name); } + if (!integrate) + m_monitor = extractMonitorSpectra(m_monitor, m_workspaceIndexes); setProperty("NormFactorWS", m_monitor); } } +/** + * Pulls the monitor spectra out of a larger workspace + * @param ws + * @param workspaceIndexes The indexes of the spectra to extract + * @return A workspace containing the spectra requested + */ +MatrixWorkspace_sptr NormaliseToMonitor::extractMonitorSpectra( + const MatrixWorkspace_sptr &ws, + const std::vector<std::size_t> &workspaceIndexes) { + IAlgorithm_sptr childAlg = createChildAlgorithm("ExtractSpectra"); + childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", ws); + childAlg->setProperty("WorkspaceIndexList", workspaceIndexes); + childAlg->executeAsChildAlg(); + MatrixWorkspace_sptr outWS = childAlg->getProperty("OutputWorkspace"); + return outWS; +} + /** Validates input properties. * @return A map of input properties as keys and (error) messages as values. */ @@ -327,7 +368,7 @@ std::map<std::string, std::string> NormaliseToMonitor::validateInputs() { * @param inputWorkspace The input workspace */ void NormaliseToMonitor::checkProperties( - const API::MatrixWorkspace_sptr &inputWorkspace) { + const MatrixWorkspace_sptr &inputWorkspace) { // Check where the monitor spectrum should come from Property *monSpec = getProperty("MonitorSpectrum"); @@ -335,6 +376,7 @@ void NormaliseToMonitor::checkProperties( Property *monID = getProperty("MonitorID"); // Is the monitor spectrum within the main input workspace const bool inWS = !monSpec->isDefault(); + m_syncScanInput = inputWorkspace->detectorInfo().isSyncScan(); // Or is it in a separate workspace bool sepWS{monWS}; // or monitor ID @@ -361,38 +403,42 @@ void NormaliseToMonitor::checkProperties( } // Do a check for common binning and store - m_commonBins = API::WorkspaceHelpers::commonBoundaries(*inputWorkspace); + m_commonBins = WorkspaceHelpers::commonBoundaries(*inputWorkspace); - int spec_num(-1); // Check the monitor spectrum or workspace and extract into new workspace - m_monitor = sepWS ? this->getMonitorWorkspace(inputWorkspace, spec_num) - : this->getInWSMonitorSpectrum(inputWorkspace, spec_num); + m_monitor = sepWS ? getMonitorWorkspace(inputWorkspace) + : getInWSMonitorSpectrum(inputWorkspace); // Check that the 'monitor' spectrum actually relates to a monitor - warn if // not try { - if (!m_monitor->spectrumInfo().isMonitor(0)) { - g_log.warning() << "The spectrum N: " << spec_num - << " in MonitorWorkspace does not refer to a monitor.\n" - << "Continuing with normalization regardless."; - } - } catch (Kernel::Exception::NotFoundError &) { - g_log.warning( - "Unable to check if the spectrum provided relates to a monitor - " - "the instrument is not fully specified.\n" - "Continuing with normalization regardless."); + const auto &monitorSpecInfo = m_monitor->spectrumInfo(); + for (const auto workspaceIndex : m_workspaceIndexes) + if (!monitorSpecInfo.isMonitor(workspaceIndex)) + g_log.warning() << "The spectrum N: " << workspaceIndex + << " in MonitorWorkspace does not refer to a monitor.\n" + << "Continuing with normalization regardless."; + } catch (Kernel::Exception::NotFoundError &e) { + g_log.warning("Unable to check if the spectrum provided relates to a " + "monitor - the instrument is not fully specified.\n " + "Continuing with normalization regardless."); + g_log.warning() << "Error was: " << e.what() << "\n"; + if (m_syncScanInput) + throw std::runtime_error("Can not continue, spectrum can not be obtained " + "for monitor workspace, but the input workspace " + "has a detector scan."); } } /** Checks and retrieves the requested spectrum out of the input workspace - * @param inputWorkspace The input workspace. - * @param spectra_num The spectra number. - * @returns A workspace containing the monitor spectrum only. - * @returns spectra number (WS ID) which is used to normalize by. + * + * @param inputWorkspace The input workspace + * @returns The unchanged input workspace (so that signature is the same as + *getMonitorWorkspace) * @throw std::runtime_error If the properties are invalid */ -API::MatrixWorkspace_sptr NormaliseToMonitor::getInWSMonitorSpectrum( - const API::MatrixWorkspace_sptr &inputWorkspace, int &spectra_num) { +MatrixWorkspace_sptr NormaliseToMonitor::getInWSMonitorSpectrum( + const MatrixWorkspace_sptr &inputWorkspace) { // this is the index of the spectra within the workspace and we need to // identify it either from DetID or from SpecID // size_t spectra_num(-1); @@ -413,12 +459,17 @@ API::MatrixWorkspace_sptr NormaliseToMonitor::getInWSMonitorSpectrum( throw std::runtime_error( "Can not find spectra, corresponding to the requested monitor ID"); } - if (indexList.size() > 1) { - throw std::runtime_error("More then one spectra corresponds to the " - "requested monitor ID, which is unheard of"); + if (indexList.size() > 1 && !m_syncScanInput) { + throw std::runtime_error("More then one spectrum corresponds to the " + "requested monitor ID. This is unexpected in a " + "non-scanning workspace."); } - spectra_num = static_cast<int>(indexList[0]); + m_workspaceIndexes = indexList; } else { // monitor spectrum is specified. + if (m_syncScanInput) + throw std::runtime_error("For a sync-scan input workspace the monitor ID " + "must be provided. Normalisation can not be " + "performed to a spectrum."); const SpectraAxis *axis = dynamic_cast<const SpectraAxis *>(inputWorkspace->getAxis(1)); if (!axis) { @@ -430,62 +481,50 @@ API::MatrixWorkspace_sptr NormaliseToMonitor::getInWSMonitorSpectrum( throw std::runtime_error("Input workspace does not contain spectrum " "number given for MonitorSpectrum"); } - spectra_num = static_cast<int>(specs[monitorSpec]); + m_workspaceIndexes = std::vector<size_t>(1, specs[monitorSpec]); } - return this->extractMonitorSpectrum(inputWorkspace, spectra_num); + return inputWorkspace; } /** Checks and retrieves the monitor spectrum out of the input workspace - * @param inputWorkspace The input workspace. - * @param wsID The workspace ID. + * @param inputWorkspace The input workspace * @returns A workspace containing the monitor spectrum only */ -API::MatrixWorkspace_sptr NormaliseToMonitor::getMonitorWorkspace( - const API::MatrixWorkspace_sptr &inputWorkspace, int &wsID) { +MatrixWorkspace_sptr NormaliseToMonitor::getMonitorWorkspace( + const MatrixWorkspace_sptr &inputWorkspace) { MatrixWorkspace_sptr monitorWS = getProperty("MonitorWorkspace"); - wsID = getProperty("MonitorWorkspaceIndex"); + const int wsID = getProperty("MonitorWorkspaceIndex"); + m_workspaceIndexes = std::vector<size_t>(1, wsID); // In this case we need to test whether the bins in the monitor workspace // match - m_commonBins = (m_commonBins && API::WorkspaceHelpers::matchingBins( + m_commonBins = (m_commonBins && WorkspaceHelpers::matchingBins( *inputWorkspace, *monitorWS, true)); // Copy the monitor spectrum because it will get changed - return this->extractMonitorSpectrum(monitorWS, wsID); + return monitorWS; } -/** Pulls the monitor spectrum out of a larger workspace - * @param WS :: The workspace containing the spectrum to extract - * @param index :: The index of the spectrum to extract - * @returns A workspace containing the single spectrum requested +/** + * @return True if the maximum or minimum values are set */ -API::MatrixWorkspace_sptr -NormaliseToMonitor::extractMonitorSpectrum(const API::MatrixWorkspace_sptr &WS, - const std::size_t index) { - IAlgorithm_sptr childAlg = createChildAlgorithm("ExtractSingleSpectrum"); - childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", WS); - childAlg->setProperty<int>("WorkspaceIndex", static_cast<int>(index)); - childAlg->executeAsChildAlg(); - MatrixWorkspace_sptr outWS = childAlg->getProperty("OutputWorkspace"); - - IAlgorithm_sptr alg = createChildAlgorithm("ConvertToMatrixWorkspace"); - alg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", outWS); - alg->executeAsChildAlg(); - outWS = alg->getProperty("OutputWorkspace"); - // Only get to here if successful - return outWS; -} - -/** Sets the maximum and minimum X values of the monitor spectrum to use for +/** + * Sets the maximum and minimum X values of the monitor spectrum to use for * integration - * @return True if the maximum or minimum values are set + * + * @param isSingleCountWorkspace Whether or not the input workspace is point + *data with single counts per spectrum + * @return True if the maximum or minimum values are set, or it is a single + *count workspace */ -bool NormaliseToMonitor::setIntegrationProps() { +bool NormaliseToMonitor::setIntegrationProps( + const bool isSingleCountWorkspace) { m_integrationMin = getProperty("IntegrationRangeMin"); m_integrationMax = getProperty("IntegrationRangeMax"); // Check if neither of these have been changed from their defaults // (EMPTY_DBL()) - if (isEmpty(m_integrationMin) && isEmpty(m_integrationMax)) { + if ((isEmpty(m_integrationMin) && isEmpty(m_integrationMax)) && + !isSingleCountWorkspace) { // Nothing has been set so the user doesn't want to use integration so let's // move on return false; @@ -510,35 +549,103 @@ bool NormaliseToMonitor::setIntegrationProps() { /** Carries out a normalization based on the integrated count of the monitor * over a range - * @param inputWorkspace The input workspace - * @param outputWorkspace The result workspace + * @param inputWorkspace The input workspace + * @param outputWorkspace The result workspace + * @param isSingleCountWorkspace Whether or not the input workspace is point + *data with single counts per spectrum */ void NormaliseToMonitor::normaliseByIntegratedCount( - const API::MatrixWorkspace_sptr &inputWorkspace, - API::MatrixWorkspace_sptr &outputWorkspace) { - // Add up all the bins so it's just effectively a single value with an error - IAlgorithm_sptr integrate = createChildAlgorithm("Integration"); - integrate->setProperty<MatrixWorkspace_sptr>("InputWorkspace", m_monitor); - integrate->setProperty("RangeLower", m_integrationMin); - integrate->setProperty("RangeUpper", m_integrationMax); - integrate->setProperty<bool>("IncludePartialBins", - getProperty("IncludePartialBins")); - - integrate->executeAsChildAlg(); - - // Get back the result - m_monitor = integrate->getProperty("OutputWorkspace"); - - // Run the divide algorithm explicitly to enable progress reporting - IAlgorithm_sptr divide = createChildAlgorithm("Divide", 0.0, 1.0); - divide->setProperty<MatrixWorkspace_sptr>("LHSWorkspace", inputWorkspace); - divide->setProperty<MatrixWorkspace_sptr>("RHSWorkspace", m_monitor); - divide->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", outputWorkspace); - - divide->executeAsChildAlg(); - - // Get back the result - outputWorkspace = divide->getProperty("OutputWorkspace"); + const MatrixWorkspace_sptr &inputWorkspace, + MatrixWorkspace_sptr &outputWorkspace, const bool isSingleCountWorkspace) { + m_monitor = extractMonitorSpectra(m_monitor, m_workspaceIndexes); + + // If single counting no need to integrate, monitor already guaranteed to be a + // single count + if (!isSingleCountWorkspace) { + // Add up all the bins so it's just effectively a series of values with + // errors + IAlgorithm_sptr integrate = createChildAlgorithm("Integration"); + integrate->setProperty<MatrixWorkspace_sptr>("InputWorkspace", m_monitor); + integrate->setProperty("RangeLower", m_integrationMin); + integrate->setProperty("RangeUpper", m_integrationMax); + integrate->setProperty<bool>("IncludePartialBins", + getProperty("IncludePartialBins")); + integrate->executeAsChildAlg(); + m_monitor = integrate->getProperty("OutputWorkspace"); + } + + EventWorkspace_sptr inputEvent = + boost::dynamic_pointer_cast<EventWorkspace>(inputWorkspace); + + if (inputEvent) { + // Run the divide algorithm explicitly to enable progress reporting + IAlgorithm_sptr divide = createChildAlgorithm("Divide", 0.0, 1.0); + divide->setProperty<MatrixWorkspace_sptr>("LHSWorkspace", inputWorkspace); + divide->setProperty<MatrixWorkspace_sptr>("RHSWorkspace", m_monitor); + divide->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", + outputWorkspace); + divide->executeAsChildAlg(); + + // Get back the result + outputWorkspace = divide->getProperty("OutputWorkspace"); + } else { + performHistogramDivision(inputWorkspace, outputWorkspace); + } +} + +/** + * This performs a similar operation to divide, but is a separate algorithm so + *that the correct spectra are used in the case of detector scans. This + *currently does not support event workspaces properly, but should be made to in + *the future. + * + * @param inputWorkspace The workspace with the spectra to divide by the monitor + * @param outputWorkspace The resulting workspace + */ +void NormaliseToMonitor::performHistogramDivision( + const MatrixWorkspace_sptr &inputWorkspace, + MatrixWorkspace_sptr &outputWorkspace) { + if (outputWorkspace != inputWorkspace) + outputWorkspace = inputWorkspace->clone(); + + size_t monitorWorkspaceIndex = 0; + + const auto &specInfo = inputWorkspace->spectrumInfo(); + for (const auto workspaceIndex : m_workspaceIndexes) { + // Errors propagated according to + // http://docs.mantidproject.org/nightly/concepts/ErrorPropagation.html#error-propagation + // This is similar to that in MantidAlgorithms::Divide + + size_t timeIndex = 0; + if (m_syncScanInput) + timeIndex = specInfo.spectrumDefinition(workspaceIndex)[0].second; + + const auto newYFactor = + 1.0 / m_monitor->histogram(monitorWorkspaceIndex).y()[0]; + const auto divisorError = + m_monitor->histogram(monitorWorkspaceIndex).e()[0]; + const double yErrorFactor = pow(divisorError * newYFactor, 2); + monitorWorkspaceIndex++; + + for (size_t i = 0; i < outputWorkspace->getNumberHistograms(); ++i) { + const auto &specDef = specInfo.spectrumDefinition(i); + + if (!spectrumDefinitionsMatchTimeIndex(specDef, timeIndex)) + continue; + + auto hist = outputWorkspace->histogram(i); + auto &yValues = hist.mutableY(); + auto &eValues = hist.mutableE(); + + for (size_t j = 0; j < yValues.size(); ++j) { + eValues[j] = newYFactor * sqrt(eValues[j] * eValues[j] + + yValues[j] * yValues[j] * yErrorFactor); + yValues[j] *= newYFactor; + } + + outputWorkspace->setHistogram(i, hist); + } + } } /** Carries out the bin-by-bin normalization @@ -546,8 +653,8 @@ void NormaliseToMonitor::normaliseByIntegratedCount( * @param outputWorkspace The result workspace */ void NormaliseToMonitor::normaliseBinByBin( - const API::MatrixWorkspace_sptr &inputWorkspace, - API::MatrixWorkspace_sptr &outputWorkspace) { + const MatrixWorkspace_sptr &inputWorkspace, + MatrixWorkspace_sptr &outputWorkspace) { EventWorkspace_sptr inputEvent = boost::dynamic_pointer_cast<EventWorkspace>(inputWorkspace); @@ -561,99 +668,116 @@ void NormaliseToMonitor::normaliseBinByBin( auto outputEvent = boost::dynamic_pointer_cast<EventWorkspace>(outputWorkspace); - // Get hold of the monitor spectrum - const auto &monX = m_monitor->binEdges(0); - auto monY = m_monitor->counts(0); - auto monE = m_monitor->countStandardDeviations(0); - // Calculate the overall normalization just the once if bins are all matching - if (m_commonBins) - this->normalisationFactor(monX, monY, monE); - - const size_t numHists = inputWorkspace->getNumberHistograms(); - auto specLength = inputWorkspace->blocksize(); - // Flag set when a division by 0 is found - bool hasZeroDivision = false; - Progress prog(this, 0.0, 1.0, numHists); - // Loop over spectra - PARALLEL_FOR_IF( - Kernel::threadSafe(*inputWorkspace, *outputWorkspace, *m_monitor)) - for (int64_t i = 0; i < int64_t(numHists); ++i) { - PARALLEL_START_INTERUPT_REGION - prog.report(); - - const auto &X = inputWorkspace->binEdges(i); - // If not rebinning, just point to our monitor spectra, otherwise create new - // vectors - - auto Y = (m_commonBins ? monY : Counts(specLength)); - auto E = (m_commonBins ? monE : CountStandardDeviations(specLength)); - - if (!m_commonBins) { - // ConvertUnits can give X vectors of all zeros - skip these, they cause - // problems - if (X.back() == 0.0 && X.front() == 0.0) + const auto &inputSpecInfo = inputWorkspace->spectrumInfo(); + const auto &monitorSpecInfo = m_monitor->spectrumInfo(); + + for (auto &workspaceIndex : m_workspaceIndexes) { + // Get hold of the monitor spectrum + const auto &monX = m_monitor->binEdges(workspaceIndex); + auto monY = m_monitor->counts(workspaceIndex); + auto monE = m_monitor->countStandardDeviations(workspaceIndex); + size_t timeIndex = 0; + if (m_syncScanInput) + timeIndex = monitorSpecInfo.spectrumDefinition(workspaceIndex)[0].second; + // Calculate the overall normalization just the once if bins are all + // matching + if (m_commonBins) + this->normalisationFactor(monX, monY, monE); + + const size_t numHists = inputWorkspace->getNumberHistograms(); + auto specLength = inputWorkspace->blocksize(); + // Flag set when a division by 0 is found + bool hasZeroDivision = false; + Progress prog(this, 0.0, 1.0, numHists); + // Loop over spectra + PARALLEL_FOR_IF( + Kernel::threadSafe(*inputWorkspace, *outputWorkspace, *m_monitor)) + for (int64_t i = 0; i < int64_t(numHists); ++i) { + PARALLEL_START_INTERUPT_REGION + prog.report(); + + const auto &specDef = inputSpecInfo.spectrumDefinition(i); + if (!spectrumDefinitionsMatchTimeIndex(specDef, timeIndex)) continue; - // Rebin the monitor spectrum to match the binning of the current data - // spectrum - VectorHelper::rebinHistogram( - monX.rawData(), monY.mutableRawData(), monE.mutableRawData(), - X.rawData(), Y.mutableRawData(), E.mutableRawData(), false); - // Recalculate the overall normalization factor - this->normalisationFactor(X, Y, E); - } - if (inputEvent) { - // --- EventWorkspace --- - EventList &outEL = outputEvent->getSpectrum(i); - outEL.divide(X.rawData(), Y.mutableRawData(), E.mutableRawData()); - } else { - // --- Workspace2D --- - auto &YOut = outputWorkspace->mutableY(i); - auto &EOut = outputWorkspace->mutableE(i); - const auto &inY = inputWorkspace->y(i); - const auto &inE = inputWorkspace->e(i); - outputWorkspace->mutableX(i) = inputWorkspace->x(i); - - // The code below comes more or less straight out of Divide.cpp - for (size_t k = 0; k < specLength; ++k) { - // Get the input Y's - const double leftY = inY[k]; - const double rightY = Y[k]; - - if (rightY == 0.0) { - hasZeroDivision = true; - } - - // Calculate result and store in local variable to avoid overwriting - // original data if output workspace is same as one of the input ones - const double newY = leftY / rightY; - - if (fabs(rightY) > 1.0e-12 && fabs(newY) > 1.0e-12) { - const double lhsFactor = (inE[k] < 1.0e-12 || fabs(leftY) < 1.0e-12) - ? 0.0 - : pow((inE[k] / leftY), 2); - const double rhsFactor = - E[k] < 1.0e-12 ? 0.0 : pow((E[k] / rightY), 2); - EOut[k] = std::abs(newY) * sqrt(lhsFactor + rhsFactor); - } - - // Now store the result - YOut[k] = newY; - } // end Workspace2D case - } // end loop over current spectrum - - PARALLEL_END_INTERUPT_REGION - } // end loop over spectra - PARALLEL_CHECK_INTERUPT_REGION - if (hasZeroDivision) { - g_log.warning() << "Division by zero in some of the bins.\n"; + const auto &X = inputWorkspace->binEdges(i); + // If not rebinning, just point to our monitor spectra, otherwise create + // new vectors + + auto Y = (m_commonBins ? monY : Counts(specLength)); + auto E = (m_commonBins ? monE : CountStandardDeviations(specLength)); + + if (!m_commonBins) { + // ConvertUnits can give X vectors of all zeros - skip these, they + // cause + // problems + if (X.back() == 0.0 && X.front() == 0.0) + continue; + // Rebin the monitor spectrum to match the binning of the current data + // spectrum + VectorHelper::rebinHistogram( + monX.rawData(), monY.mutableRawData(), monE.mutableRawData(), + X.rawData(), Y.mutableRawData(), E.mutableRawData(), false); + // Recalculate the overall normalization factor + this->normalisationFactor(X, Y, E); + } + + if (inputEvent) { + // --- EventWorkspace --- + EventList &outEL = outputEvent->getSpectrum(i); + outEL.divide(X.rawData(), Y.mutableRawData(), E.mutableRawData()); + } else { + // --- Workspace2D --- + auto &YOut = outputWorkspace->mutableY(i); + auto &EOut = outputWorkspace->mutableE(i); + const auto &inY = inputWorkspace->y(i); + const auto &inE = inputWorkspace->e(i); + outputWorkspace->setSharedX(i, inputWorkspace->sharedX(i)); + + // The code below comes more or less straight out of Divide.cpp + for (size_t k = 0; k < specLength; ++k) { + // Get the input Y's + const double leftY = inY[k]; + const double rightY = Y[k]; + + if (rightY == 0.0) { + hasZeroDivision = true; + } + + // Calculate result and store in local variable to avoid overwriting + // original data if output workspace is same as one of the input + // ones + const double newY = leftY / rightY; + + if (fabs(rightY) > 1.0e-12 && fabs(newY) > 1.0e-12) { + const double lhsFactor = (inE[k] < 1.0e-12 || fabs(leftY) < 1.0e-12) + ? 0.0 + : pow((inE[k] / leftY), 2); + const double rhsFactor = + E[k] < 1.0e-12 ? 0.0 : pow((E[k] / rightY), 2); + EOut[k] = std::abs(newY) * sqrt(lhsFactor + rhsFactor); + } + + // Now store the result + YOut[k] = newY; + } // end Workspace2D case + } // end loop over current spectrum + + PARALLEL_END_INTERUPT_REGION + } // end loop over spectra + PARALLEL_CHECK_INTERUPT_REGION + + if (hasZeroDivision) { + g_log.warning() << "Division by zero in some of the bins.\n"; + } + if (inputEvent) + outputEvent->clearMRU(); } - if (inputEvent) - outputEvent->clearMRU(); } /** Calculates the overall normalization factor. - * This multiplies result by (bin width * sum of monitor counts) / total frame + * This multiplies result by (bin width * sum of monitor counts) / total + * frame * width. * @param X The BinEdges of the workspace * @param Y The Counts of the workspace diff --git a/Framework/Algorithms/test/NormaliseToMonitorTest.h b/Framework/Algorithms/test/NormaliseToMonitorTest.h index 8f9ff5cdde6f6d36cdd5bc3812ea5309fe12fda2..83d65fb877586dbd6e9675c815591ebbdd82b5b9 100644 --- a/Framework/Algorithms/test/NormaliseToMonitorTest.h +++ b/Framework/Algorithms/test/NormaliseToMonitorTest.h @@ -7,16 +7,22 @@ #include "MantidDataObjects/WorkspaceCreation.h" #include "MantidAPI/Axis.h" #include "MantidAPI/FrameworkManager.h" +#include "MantidAPI/SpectrumInfo.h" #include "MantidGeometry/Instrument.h" +#include "MantidGeometry/Instrument/DetectorInfo.h" #include "MantidHistogramData/BinEdges.h" #include "MantidHistogramData/Counts.h" +#include "MantidHistogramData/LinearGenerator.h" +#include "MantidIndexing/IndexInfo.h" #include "MantidKernel/UnitFactory.h" #include "MantidTestHelpers/WorkspaceCreationHelper.h" +using namespace Mantid; using namespace Mantid::Kernel; using namespace Mantid::API; using namespace Mantid::Algorithms; using namespace Mantid::DataObjects; +using namespace Mantid::HistogramData; using Mantid::Geometry::Instrument; // Anonymous namespace for shared methods between unit and performance test @@ -423,8 +429,6 @@ public: } void testMonitorWorkspaceNotInADSWorks() { - using namespace Mantid::DataObjects; - using namespace Mantid::HistogramData; BinEdges xs{-1.0, 1.0}; Counts ys{1.0}; Histogram h(xs, ys); @@ -441,6 +445,172 @@ public: TS_ASSERT_THROWS_NOTHING(alg.execute()) TS_ASSERT(alg.isExecuted()) } + + void test_with_scanning_workspace_bin_by_bin() { + auto testWS = makeTestDetectorScanWorkspace(); + + NormaliseToMonitor alg; + alg.setChild(true); + TS_ASSERT_THROWS_NOTHING(alg.initialize()) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", testWS)) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("OutputWorkspace", "outputWS")) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("MonitorID", "9")) + TS_ASSERT_THROWS_NOTHING(alg.execute()) + TS_ASSERT(alg.isExecuted()) + + MatrixWorkspace_sptr outWS = alg.getProperty("OutputWorkspace"); + + const auto &specOutInfo = outWS->spectrumInfo(); + for (size_t i = 0; i < specOutInfo.size(); ++i) { + const auto &yValues = outWS->histogram(i).y(); + for (size_t j = 0; j < yValues.size(); ++j) { + if (specOutInfo.isMonitor(i)) + TS_ASSERT_DELTA(yValues[j], 3.0, 1e-12) + else + TS_ASSERT_DELTA(yValues[j], 6.0 / double(j + 1), 1e-12) + } + } + } + + void test_with_scanning_workspace_integration_range() { + auto testWS = makeTestDetectorScanWorkspace(); + + NormaliseToMonitor alg; + alg.setChild(true); + TS_ASSERT_THROWS_NOTHING(alg.initialize()) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", testWS)) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("OutputWorkspace", "outputWS")) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("MonitorID", "9")) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("IntegrationRangeMin", "-1")) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("IntegrationRangeMax", "99")) + TS_ASSERT_THROWS_NOTHING(alg.execute()) + TS_ASSERT(alg.isExecuted()) + + MatrixWorkspace_sptr outWS = alg.getProperty("OutputWorkspace"); + + const auto &specOutInfo = outWS->spectrumInfo(); + for (size_t i = 0; i < specOutInfo.size(); ++i) { + const auto &yValues = outWS->histogram(i).y(); + for (size_t j = 0; j < yValues.size(); ++j) { + if (specOutInfo.isMonitor(i)) + TS_ASSERT_DELTA(yValues[j], double(j + 1) / 15.0, 1e-12) + else + TS_ASSERT_DELTA(yValues[j], 2.0 / 15.0, 1e-12) + } + } + } + + void test_with_async_scan_workspace_throws() { + const size_t N_DET = 10; + const size_t N_BINS = 5; + + // Set up 2 workspaces to be merged + auto ws1 = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument( + N_DET, N_BINS, true); + auto ws2 = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument( + N_DET, N_BINS, true); + auto &detInfo1 = ws1->mutableDetectorInfo(); + auto &detInfo2 = ws2->mutableDetectorInfo(); + for (size_t i = 0; i < N_DET; ++i) { + detInfo1.setScanInterval(i, {10, 20}); + detInfo2.setScanInterval(i, {20, 30}); + } + // Merge + auto merged = WorkspaceFactory::Instance().create(ws1, 2 * N_DET); + auto &detInfo = merged->mutableDetectorInfo(); + detInfo.merge(detInfo2); + merged->setIndexInfo(Indexing::IndexInfo(merged->getNumberHistograms())); + + NormaliseToMonitor alg; + alg.setChild(true); + alg.setRethrows(true); + TS_ASSERT_THROWS_NOTHING(alg.initialize()) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", merged)) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("OutputWorkspace", "outputWS")) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("MonitorID", "9")) + TS_ASSERT_THROWS_EQUALS( + alg.execute(), std::runtime_error & e, std::string(e.what()), + "More then one spectrum corresponds to the requested monitor ID. This " + "is unexpected in a non-scanning workspace.") + } + + void test_with_single_count_point_data_workspace() { + const size_t N_DET = 10; + const size_t N_BINS = 1; + auto testWS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument( + N_DET, N_BINS, true); + + Points points({0.0}); + + const auto &specInfo = testWS->spectrumInfo(); + for (size_t i = 0; i < specInfo.size(); ++i) { + auto hist = testWS->histogram(i); + auto &yValues = hist.mutableY(); + for (size_t j = 0; j < yValues.size(); ++j) { + if (specInfo.isMonitor(i)) + yValues[j] = 3.0; + else + TS_ASSERT_EQUALS(yValues[j], 2.0) + } + testWS->setHistogram(i, hist); + testWS->setPoints(i, points); + } + + NormaliseToMonitor alg; + alg.setChild(true); + TS_ASSERT_THROWS_NOTHING(alg.initialize()) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", testWS)) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("OutputWorkspace", "outputWS")) + TS_ASSERT_THROWS_NOTHING(alg.setProperty("MonitorID", "9")) + TS_ASSERT_THROWS_NOTHING(alg.execute()) + TS_ASSERT(alg.isExecuted()) + + MatrixWorkspace_sptr outWS = alg.getProperty("OutputWorkspace"); + + const auto &specOutInfo = outWS->spectrumInfo(); + for (size_t i = 0; i < specOutInfo.size(); ++i) { + const auto &yValues = outWS->histogram(i).y(); + for (size_t j = 0; j < yValues.size(); ++j) { + if (specOutInfo.isMonitor(i)) + TS_ASSERT_DELTA(yValues[j], 1.0, 1e-12) + else + TS_ASSERT_DELTA(yValues[j], 2.0 / 3.0, 1e-12) + } + } + } + +private: + MatrixWorkspace_sptr makeTestDetectorScanWorkspace() { + const size_t N_DET = 10; + const size_t N_BINS = 5; + const size_t N_TIMEINDEXES = 5; + + MatrixWorkspace_sptr testWS = WorkspaceCreationHelper:: + create2DDetectorScanWorkspaceWithFullInstrument( + N_DET, N_BINS, N_TIMEINDEXES, 0, 1, true); + + double x0 = 0.0; + double deltax = 1.0; + + BinEdges x(N_BINS + 1, LinearGenerator(x0, deltax)); + Counts y(N_BINS, 2.0); + Counts yMonitor(N_BINS, 2.0); + for (size_t j = 0; j < y.size(); ++j) { + y.mutableRawData()[j] = 2.0; + yMonitor.mutableRawData()[j] = double(j + 1); + } + CountStandardDeviations e(N_BINS, M_SQRT2); + const auto &specInfo = testWS->spectrumInfo(); + for (size_t i = 0; i < N_DET * N_TIMEINDEXES; i++) { + testWS->setBinEdges(i, x); + if (specInfo.isMonitor(i)) + testWS->setCounts(i, yMonitor); + else + testWS->setCounts(i, y); + testWS->setCountStandardDeviations(i, e); + } + return testWS; + } }; class NormaliseToMonitorTestPerformance : public CxxTest::TestSuite { diff --git a/Framework/Beamline/inc/MantidBeamline/DetectorInfo.h b/Framework/Beamline/inc/MantidBeamline/DetectorInfo.h index ea98ed3a73e38f9e53f18622eda671fa30266140..2406255d6ae3d23abe7a93b8d3eb150e9e16823a 100644 --- a/Framework/Beamline/inc/MantidBeamline/DetectorInfo.h +++ b/Framework/Beamline/inc/MantidBeamline/DetectorInfo.h @@ -71,6 +71,7 @@ public: size_t size() const; size_t scanSize() const; bool isScanning() const; + bool isSyncScan() const; bool isMonitor(const size_t index) const; bool isMonitor(const std::pair<size_t, size_t> &index) const; diff --git a/Framework/Beamline/src/DetectorInfo.cpp b/Framework/Beamline/src/DetectorInfo.cpp index b321392485ce7d7d34727a8de98ab0a498049460..ccac4e77ab7d4fae1f8d1b9e681742b947018480 100644 --- a/Framework/Beamline/src/DetectorInfo.cpp +++ b/Framework/Beamline/src/DetectorInfo.cpp @@ -102,6 +102,12 @@ size_t DetectorInfo::scanSize() const { return m_positions->size(); } +/** + * Returns true if all of the detectors all have the same scan interval. Will + * return false if DetectorInfo is not scanning. + */ +bool DetectorInfo::isSyncScan() const { return isScanning() && m_isSyncScan; } + /// Returns true if the detector with given detector index is a monitor. bool DetectorInfo::isMonitor(const size_t index) const { // No check for time dependence since monitor flags are not time dependent. diff --git a/Framework/Beamline/test/DetectorInfoTest.h b/Framework/Beamline/test/DetectorInfoTest.h index 958a55c458f3fb937192e26f4365991cc9a03f1f..971d886b7d13501a94a35be29c99baebd4f256a3 100644 --- a/Framework/Beamline/test/DetectorInfoTest.h +++ b/Framework/Beamline/test/DetectorInfoTest.h @@ -545,6 +545,7 @@ public: b.setScanInterval(0, interval2); TS_ASSERT_THROWS_NOTHING(a.merge(b)); TS_ASSERT(a.isScanning()); + TS_ASSERT(!a.isSyncScan()); TS_ASSERT(!a.isEquivalent(b)); TS_ASSERT_EQUALS(a.size(), 2); TS_ASSERT_EQUALS(a.scanSize(), 3); @@ -572,6 +573,7 @@ public: b.setScanInterval(interval2); TS_ASSERT_THROWS_NOTHING(a.merge(b)); TS_ASSERT(a.isScanning()); + TS_ASSERT(a.isSyncScan()); TS_ASSERT(!a.isEquivalent(b)); TS_ASSERT_EQUALS(a.size(), 2); TS_ASSERT_EQUALS(a.scanSize(), 4); @@ -632,6 +634,7 @@ public: TS_ASSERT_THROWS_NOTHING(a.merge(b)); TS_ASSERT_THROWS_NOTHING(a.merge(c)); TS_ASSERT(a.isScanning()); + TS_ASSERT(!a.isSyncScan()); TS_ASSERT(!a.isEquivalent(b)); TS_ASSERT(!a.isEquivalent(c)); TS_ASSERT_EQUALS(a.size(), 2); @@ -666,6 +669,7 @@ public: TS_ASSERT_THROWS_NOTHING(a.merge(b)); TS_ASSERT_THROWS_NOTHING(a.merge(c)); TS_ASSERT(a.isScanning()); + TS_ASSERT(a.isSyncScan()); TS_ASSERT(!a.isEquivalent(b)); TS_ASSERT(!a.isEquivalent(c)); TS_ASSERT_EQUALS(a.size(), 2); diff --git a/Framework/Geometry/inc/MantidGeometry/Instrument/DetectorInfo.h b/Framework/Geometry/inc/MantidGeometry/Instrument/DetectorInfo.h index 006bf2878dd37fae9abf1a30a8bc2f734a3c655e..18924dc6f9c4c58bd5cfba333c228106e53d708f 100644 --- a/Framework/Geometry/inc/MantidGeometry/Instrument/DetectorInfo.h +++ b/Framework/Geometry/inc/MantidGeometry/Instrument/DetectorInfo.h @@ -78,6 +78,7 @@ public: size_t size() const; size_t scanSize() const; bool isScanning() const; + bool isSyncScan() const; bool isMonitor(const size_t index) const; bool isMonitor(const std::pair<size_t, size_t> &index) const; diff --git a/Framework/Geometry/src/Instrument/DetectorInfo.cpp b/Framework/Geometry/src/Instrument/DetectorInfo.cpp index 6c24ac4f7a624165dde3f7750604df09bdf38cef..d2b98326e484949b50555dc7ac2456254804b2e6 100644 --- a/Framework/Geometry/src/Instrument/DetectorInfo.cpp +++ b/Framework/Geometry/src/Instrument/DetectorInfo.cpp @@ -90,6 +90,10 @@ size_t DetectorInfo::scanSize() const { return m_detectorInfo->scanSize(); } /// Returns true if the beamline has scanning detectors. bool DetectorInfo::isScanning() const { return m_detectorInfo->isScanning(); } +/// Returns true if the beamline has scanning detectors and they have all the +/// same scan intervals. +bool DetectorInfo::isSyncScan() const { return m_detectorInfo->isSyncScan(); } + /// Returns true if the detector is a monitor. bool DetectorInfo::isMonitor(const size_t index) const { return m_detectorInfo->isMonitor(index); diff --git a/Framework/TestHelpers/inc/MantidTestHelpers/FakeObjects.h b/Framework/TestHelpers/inc/MantidTestHelpers/FakeObjects.h index a1e3bac4969d22e40f1744bb21a7b6fd1037c9ac..ba21bb218f3ae99c49d057867abbdaa64c82af4d 100644 --- a/Framework/TestHelpers/inc/MantidTestHelpers/FakeObjects.h +++ b/Framework/TestHelpers/inc/MantidTestHelpers/FakeObjects.h @@ -398,4 +398,15 @@ protected: throw std::runtime_error("void_pointer const not implemented"); } }; + +class VariableBinThrowingTester : public AxeslessWorkspaceTester { + size_t blocksize() const override { + if (getSpectrum(0).dataY().size() == getSpectrum(1).dataY().size()) + return getSpectrum(0).dataY().size(); + else + throw std::length_error("Mismatched bins sizes"); + + return 0; + } +}; #endif /* FAKEOBJECTS_H_ */ diff --git a/docs/source/algorithms/NormaliseToMonitor-v1.rst b/docs/source/algorithms/NormaliseToMonitor-v1.rst index 1580d16a7d168417a5a084a853e279b5e64c475c..520ea336220f699f86a6e843f472952ba5049163 100644 --- a/docs/source/algorithms/NormaliseToMonitor-v1.rst +++ b/docs/source/algorithms/NormaliseToMonitor-v1.rst @@ -55,20 +55,26 @@ technically dimensionless. Restrictions on the input workspace ################################### -The data must be histogram, non-distribution data. +The data must be histogram, non-distribution data. The exception to the histogram requirement is for workspaces that contain point data with +a single count per spectrum. In this case the normalisation is performed by dividing every spectrum by the monitor counts, taking into +account the error on the monitor counts. + +Detector Scan Workspaces +######################## + +Workspaces that have scanning detectors are supported by this algorithm, both for bin-by-bin mode and normlisation by integrated +count. The only option for specifying the monitor is by 'MonitorID', attempting to use 'MonitorSpectrum' or MonitorWorkspaceIndex' +will throw an error. In this case the 'NormFactorWS' output will contain a monitor spectrum for each time index. Child Algorithms used ##################### -The :ref:`algm-ExtractSingleSpectrum` algorithm is used +The :ref:`algm-ExtractSpectra` algorithm is used to pull out the monitor spectrum if it's part of the InputWorkspace or MonitorWorkspace. For the 'integrated range' option, the :ref:`algm-Integration` algorithm is used to integrate the monitor spectrum. -In both cases, the :ref:`algm-Divide` algorithm is used to perform the -normalisation. - Usage ----- .. include:: ../usagedata-note.txt diff --git a/docs/source/release/v3.12.0/framework.rst b/docs/source/release/v3.12.0/framework.rst index 08092d3c6c9c816a07a58c11c19f788593d238f7..e9096578e0629614f938efc44c79fe1148ae399e 100644 --- a/docs/source/release/v3.12.0/framework.rst +++ b/docs/source/release/v3.12.0/framework.rst @@ -22,6 +22,7 @@ Please contact the Mantid Team if you experience any further problems as a resul Algorithms ---------- +:ref:`NormaliseToMonitor <algm-NormaliseToMonitor>` now supports workspaces with detector scans and workspaces with single-count point data. - :ref:`CreateWorkspace <algm-CreateWorkspace>` will no longer create a default (and potentially wrong) mapping from spectra to detectors, unless a parent workspace is given. This change ensures that accidental bad mappings that could lead to corrupted data are not created silently anymore. This change does *not* affect the use of this algorithm if: (1) a parent workspace is given, or (2) no instrument is loaded into to workspace at a later point, or (3) an instrument is loaded at a later point but ``LoadInstrument`` is used with ``RewriteSpectraMapping=True``. See also the algorithm documentation for details. - :ref:`Fit <algm-Fit>` will now respect excluded ranges when ``CostFunction = 'Unweighted least squares'``.