diff --git a/Framework/Muon/CMakeLists.txt b/Framework/Muon/CMakeLists.txt index 9607b7420874cd78da2e550732eff8b1f8e976af..7d495a27ea01718a488e5f5a80c767afe4c3ce5e 100644 --- a/Framework/Muon/CMakeLists.txt +++ b/Framework/Muon/CMakeLists.txt @@ -13,6 +13,7 @@ set ( SRC_FILES src/MuonAlgorithmHelper.cpp src/MuonAsymmetryHelper.cpp src/MuonGroupDetectors.cpp + src/MuonPreProcess.cpp src/PhaseQuadMuon.cpp src/PlotAsymmetryByLogValue.cpp src/RemoveExpDecay.cpp @@ -34,6 +35,7 @@ set ( INC_FILES inc/MantidMuon/MuonAlgorithmHelper.h inc/MantidMuon/MuonAsymmetryHelper.h inc/MantidMuon/MuonGroupDetectors.h + inc/MantidMuon/MuonPreProcess.h inc/MantidMuon/PhaseQuadMuon.h inc/MantidMuon/PlotAsymmetryByLogValue.h inc/MantidMuon/RemoveExpDecay.h @@ -54,6 +56,7 @@ set ( TEST_FILES MuonAlgorithmHelperTest.h EstimateMuonAsymmetryFromCountsTest.h MuonGroupDetectorsTest.h + MuonPreProcessTest.h PhaseQuadMuonTest.h PlotAsymmetryByLogValueTest.h RemoveExpDecayTest.h diff --git a/Framework/Muon/inc/MantidMuon/MuonPreProcess.h b/Framework/Muon/inc/MantidMuon/MuonPreProcess.h new file mode 100644 index 0000000000000000000000000000000000000000..5d4e721e17ae89baabae146f7786f4608dd334bb --- /dev/null +++ b/Framework/Muon/inc/MantidMuon/MuonPreProcess.h @@ -0,0 +1,72 @@ +// Mantid Repository : https://github.com/mantidproject/mantid +// +// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI, +// NScD Oak Ridge National Laboratory, European Spallation Source +// & Institut Laue - Langevin +// SPDX - License - Identifier: GPL - 3.0 + +#ifndef MANTID_MUON_MUONPREPROCESS_H_ +#define MANTID_MUON_MUONPREPROCESS_H_ + +#include "MantidAPI/Algorithm.h" +#include "MantidAPI/WorkspaceGroup.h" +#include "MantidDataObjects/TableWorkspace.h" + +using namespace Mantid::API; +using namespace Mantid::DataObjects; + +namespace Mantid { +namespace Muon { + +class DLLExport MuonPreProcess : public API::Algorithm { +public: + MuonPreProcess() : API::Algorithm() {} + ~MuonPreProcess() {} + + const std::string name() const override { return "MuonPreProcess"; } + int version() const override { return (1); } + const std::string category() const override { return "Muon\\DataHandling"; } + const std::string summary() const override { + return "Perform a series of common analysis pre-processing steps on Muon " + "data. Sample logs are modified to record the input parameters."; + } + const std::vector<std::string> seeAlso() const override { + return {"MuonProcess", "ApplyDeadTimeCorr", "ChangeBinOffset", + "CropWorkspace", "Rebin"}; + } + +private: + void init() override; + void exec() override; + + /// Apply a series of corrections ; DTC, offset, rebin, crop + WorkspaceGroup_sptr correctWorkspaces(WorkspaceGroup_sptr wsGroup); + MatrixWorkspace_sptr correctWorkspace(MatrixWorkspace_sptr ws); + + MatrixWorkspace_sptr applyDTC(MatrixWorkspace_sptr ws, + TableWorkspace_sptr dt); + + MatrixWorkspace_sptr applyTimeOffset(MatrixWorkspace_sptr ws, + const double &offset); + + MatrixWorkspace_sptr applyCropping(MatrixWorkspace_sptr ws, + const double &xMin, const double &xMax); + + MatrixWorkspace_sptr applyRebinning(MatrixWorkspace_sptr ws, + const std::vector<double> &rebinArgs); + + MatrixWorkspace_sptr cloneWorkspace(MatrixWorkspace_sptr ws); + + /// Add the correction inputs into the logs + void addPreProcessSampleLogs(WorkspaceGroup_sptr group); + + /// Perform validation of inputs to the algorithm + std::map<std::string, std::string> validateInputs() override; + + /// Allow WorkspaceGroup property to function correctly. + bool checkGroups() override; +}; + +} // namespace Muon +} // namespace Mantid + +#endif /* MANTID_MUON_MUONPREPROCESS_H_ */ diff --git a/Framework/Muon/src/MuonPreProcess.cpp b/Framework/Muon/src/MuonPreProcess.cpp new file mode 100644 index 0000000000000000000000000000000000000000..2d483498f0ce14e5861f66494befa245e484fef2 --- /dev/null +++ b/Framework/Muon/src/MuonPreProcess.cpp @@ -0,0 +1,318 @@ +// Mantid Repository : https://github.com/mantidproject/mantid +// +// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI, +// NScD Oak Ridge National Laboratory, European Spallation Source +// & Institut Laue - Langevin +// SPDX - License - Identifier: GPL - 3.0 + +#include "MantidMuon/MuonPreProcess.h" +#include "MantidAPI/Algorithm.h" +#include "MantidAPI/AlgorithmManager.h" +#include "MantidAPI/DataProcessorAlgorithm.h" +#include "MantidAPI/FileProperty.h" +#include "MantidAPI/MatrixWorkspace.h" +#include "MantidAPI/WorkspaceFactory.h" +#include "MantidAPI/WorkspaceGroup.h" +#include "MantidAPI/WorkspaceGroup_fwd.h" +#include "MantidDataObjects/TableWorkspace.h" +#include "MantidHistogramData/HistogramMath.h" +#include "MantidKernel/ArrayProperty.h" +#include "MantidKernel/CompositeValidator.h" +#include "MantidKernel/ListValidator.h" +#include "MantidKernel/MandatoryValidator.h" +#include "MantidKernel/System.h" +#include "MantidKernel/make_unique.h" +#include "MantidMuon/MuonAlgorithmHelper.h" + +using namespace Mantid::API; +using namespace Mantid::DataObjects; +using namespace Mantid::Kernel; + +namespace {} // namespace + +namespace Mantid { +namespace Muon { + +// Register the algorithm into the AlgorithmFactory +DECLARE_ALGORITHM(MuonPreProcess) + +void MuonPreProcess::init() { + + declareProperty( + make_unique<WorkspaceProperty<Workspace>>( + "InputWorkspace", "", Direction::Input, PropertyMode::Mandatory), + "Input workspace containing data from detectors that the " + "grouping/pairing will be applied to."); + + declareProperty( + Mantid::Kernel::make_unique<WorkspaceProperty<WorkspaceGroup>>( + "OutputWorkspace", "", Direction::Output), + "The output workspace group with all corrections applied. For single " + "period data, a group is returned with a single workspace."); + + declareProperty("TimeMin", EMPTY_DBL(), + "Start time for the data in micro seconds.", + Direction::Input); + + declareProperty("TimeMax", EMPTY_DBL(), + "End time for the data in micro seconds.", Direction::Input); + + declareProperty( + make_unique<ArrayProperty<double>>("RebinArgs", Direction::Input), + "Parameters used for rebinning. If empty - rebinning is not done."); + + declareProperty("TimeOffset", EMPTY_DBL(), + "Shift the times of all data by a fixed amount (in micro " + "seconds). The value given corresponds to the bin that will " + "become 0.0 seconds.", + Direction::Input); + + declareProperty( + make_unique<WorkspaceProperty<TableWorkspace>>( + "DeadTimeTable", "", Direction::Input, PropertyMode::Optional), + "TableWorkspace with dead time information, used to apply dead time " + "correction."); + + std::string analysisGrp("Analysis Options"); + setPropertyGroup("TimeMin", analysisGrp); + setPropertyGroup("TimeMax", analysisGrp); + setPropertyGroup("RebinArgs", analysisGrp); + setPropertyGroup("TimeOffset", analysisGrp); + setPropertyGroup("DeadTimeTable", analysisGrp); +} + +std::map<std::string, std::string> MuonPreProcess::validateInputs() { + std::map<std::string, std::string> errors; + + double tmin = this->getProperty("TimeMin"); + double tmax = this->getProperty("TimeMax"); + if (tmin != EMPTY_DBL() && tmax != EMPTY_DBL()) { + if (tmin > tmax) { + errors["TimeMin"] = "TimeMin > TimeMax"; + } + if (tmin != EMPTY_DBL() && tmin == tmax) { + errors["TimeMin"] = "TimeMin and TimeMax must be different"; + } + } + + // Checks for dead time table + Workspace_sptr inputWS = this->getProperty("InputWorkspace"); + if (auto ws = boost::dynamic_pointer_cast<MatrixWorkspace>(inputWS)) { + TableWorkspace_sptr deadTimeTable = this->getProperty("DeadTimeTable"); + if (deadTimeTable) { + if (deadTimeTable->rowCount() > ws->getNumberHistograms()) { + errors["DeadTimeTable"] = "DeadTimeTable must have as many rows as " + "there are spectra in InputWorkspace."; + } + } + } + + if (auto ws = boost::dynamic_pointer_cast<WorkspaceGroup>(inputWS)) { + if (ws->getNumberOfEntries() == 0) { + errors["InputWorkspace"] = "Input WorkspaceGroup is empty."; + } else { + auto nSpectra = + boost::dynamic_pointer_cast<MatrixWorkspace>(ws->getItem(0)) + ->getNumberHistograms(); + for (int index = 1; index < ws->getNumberOfEntries(); index++) { + if (boost::dynamic_pointer_cast<MatrixWorkspace>(ws->getItem(index)) + ->getNumberHistograms() != nSpectra) { + errors["InputWorkspace"] = + "Numbers of spectra should be identical across all workspaces in " + "the workspace group."; + } + } + } + } + + return errors; +} + +void MuonPreProcess::exec() { + this->setRethrows(true); + + Workspace_sptr inputWS = getProperty("InputWorkspace"); + + // If single period, add workspace to a group + auto allPeriodsWS = boost::make_shared<WorkspaceGroup>(); + if (auto ws = boost::dynamic_pointer_cast<MatrixWorkspace>(inputWS)) { + allPeriodsWS->addWorkspace(ws); + } else if (auto group = + boost::dynamic_pointer_cast<WorkspaceGroup>(inputWS)) { + allPeriodsWS = group; + } + + allPeriodsWS = correctWorkspaces(allPeriodsWS); + + addPreProcessSampleLogs(allPeriodsWS); + + setProperty("OutputWorkspace", allPeriodsWS); +} + +/** + * Applies offset, crops and rebin the workspaces in the group according to + * specified params. + * @param wsGroup :: Workspaces to correct + * @return Corrected workspaces + */ +WorkspaceGroup_sptr +MuonPreProcess::correctWorkspaces(WorkspaceGroup_sptr wsGroup) { + WorkspaceGroup_sptr outWS = boost::make_shared<WorkspaceGroup>(); + for (auto &&workspace : *wsGroup) { + if (auto ws = boost::dynamic_pointer_cast<MatrixWorkspace>(workspace)) { + outWS->addWorkspace(correctWorkspace(ws)); + } + } + return outWS; +} + +/** + * Applies offset, crops and rebin the workspace according to specified params. + * @param ws :: Workspace to correct + * @return Corrected workspace + */ +MatrixWorkspace_sptr MuonPreProcess::correctWorkspace(MatrixWorkspace_sptr ws) { + double offset = getProperty("TimeOffset"); + double xMin = getProperty("TimeMin"); + double xMax = getProperty("TimeMax"); + std::vector<double> rebinParams = getProperty("RebinArgs"); + TableWorkspace_sptr deadTimes = getProperty("DeadTimeTable"); + + ws = applyDTC(ws, deadTimes); + ws = applyTimeOffset(ws, offset); + ws = applyCropping(ws, xMin, xMax); + ws = applyRebinning(ws, rebinParams); + + if (deadTimes == nullptr && offset == EMPTY_DBL() && + (xMin == EMPTY_DBL() || xMax == EMPTY_DBL()) && rebinParams.empty()) { + ws = cloneWorkspace(ws); + } + + return ws; +} + +MatrixWorkspace_sptr MuonPreProcess::applyDTC(MatrixWorkspace_sptr ws, + TableWorkspace_sptr dt) { + if (dt != nullptr) { + IAlgorithm_sptr dtc = this->createChildAlgorithm("ApplyDeadTimeCorr"); + dtc->setProperty("InputWorkspace", ws); + dtc->setProperty("DeadTimeTable", dt); + dtc->execute(); + return dtc->getProperty("OutputWorkspace"); + } else { + return ws; + } +} + +MatrixWorkspace_sptr MuonPreProcess::applyTimeOffset(MatrixWorkspace_sptr ws, + const double &offset) { + if (offset != EMPTY_DBL()) { + IAlgorithm_sptr changeOffset = createChildAlgorithm("ChangeBinOffset"); + changeOffset->setProperty("InputWorkspace", ws); + changeOffset->setProperty("Offset", offset); + changeOffset->execute(); + return changeOffset->getProperty("OutputWorkspace"); + } else { + return ws; + } +} + +MatrixWorkspace_sptr MuonPreProcess::applyCropping(MatrixWorkspace_sptr ws, + const double &xMin, + const double &xMax) { + if (xMin != EMPTY_DBL() || xMax != EMPTY_DBL()) { + IAlgorithm_sptr crop = createChildAlgorithm("CropWorkspace"); + crop->setProperty("InputWorkspace", ws); + if (xMin != EMPTY_DBL()) + crop->setProperty("Xmin", xMin); + if (xMax != EMPTY_DBL()) + crop->setProperty("Xmax", xMax); + crop->execute(); + return crop->getProperty("OutputWorkspace"); + } else { + return ws; + } +} + +MatrixWorkspace_sptr +MuonPreProcess::applyRebinning(MatrixWorkspace_sptr ws, + const std::vector<double> &rebinArgs) { + if (!rebinArgs.empty()) { + IAlgorithm_sptr rebin = createChildAlgorithm("Rebin"); + rebin->setProperty("InputWorkspace", ws); + rebin->setProperty("Params", rebinArgs); + rebin->setProperty("FullBinsOnly", false); + rebin->execute(); + return rebin->getProperty("OutputWorkspace"); + } else { + return ws; + } +} + +MatrixWorkspace_sptr MuonPreProcess::cloneWorkspace(MatrixWorkspace_sptr ws) { + IAlgorithm_sptr cloneWorkspace = this->createChildAlgorithm("CloneWorkspace"); + cloneWorkspace->setProperty("InputWorkspace", ws); + cloneWorkspace->execute(); + Workspace_sptr wsClone = cloneWorkspace->getProperty("OutputWorkspace"); + return boost::dynamic_pointer_cast<MatrixWorkspace>(wsClone); +} + +void MuonPreProcess::addPreProcessSampleLogs(WorkspaceGroup_sptr group) { + const std::string numPeriods = std::to_string(group->getNumberOfEntries()); + + for (auto &&workspace : *group) { + + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_periods", numPeriods); + + std::vector<double> rebinArgs = getProperty("RebinArgs"); + if (rebinArgs.empty()) { + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_rebin_args", ""); + } else { + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_rebin_args", getPropertyValue("RebinArgs")); + } + + double xmin = getProperty("TimeMin"); + if (xmin == EMPTY_DBL()) { + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_crop_x_min", ""); + } else { + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_crop_x_min", getPropertyValue("TimeMin")); + } + + double xmax = getProperty("TimeMax"); + if (xmax == EMPTY_DBL()) { + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_crop_x_max", ""); + } else { + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_crop_x_max", getPropertyValue("TimeMax")); + } + + double offset = getProperty("TimeOffset"); + if (offset == EMPTY_DBL()) { + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_time_offset", ""); + } else { + MuonAlgorithmHelper::addSampleLog( + boost::dynamic_pointer_cast<MatrixWorkspace>(workspace), + "analysis_time_offset", getPropertyValue("TimeOffset")); + } + } +} + +// Allow WorkspaceGroup property to function correctly. +bool MuonPreProcess::checkGroups() { return false; } + +} // namespace Muon +} // namespace Mantid diff --git a/Framework/Muon/test/MuonPreProcessTest.h b/Framework/Muon/test/MuonPreProcessTest.h new file mode 100644 index 0000000000000000000000000000000000000000..1dcf393af1f6e7c9ccced9b78e4021c94f572e3c --- /dev/null +++ b/Framework/Muon/test/MuonPreProcessTest.h @@ -0,0 +1,435 @@ +// Mantid Repository : https://github.com/mantidproject/mantid +// +// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI, +// NScD Oak Ridge National Laboratory, European Spallation Source +// & Institut Laue - Langevin +// SPDX - License - Identifier: GPL - 3.0 + +#ifndef MANTID_MUON_MUONPREPROCESSTEST_H_ +#define MANTID_MUON_MUONPREPROCESSTEST_H_ + +#include "MantidAPI/AlgorithmManager.h" +#include "MantidAPI/FrameworkManager.h" +#include "MantidMuon/MuonPreProcess.h" +#include "MantidTestHelpers/MuonWorkspaceCreationHelper.h" + +#include <cxxtest/TestSuite.h> + +using namespace Mantid; +using namespace Mantid::Kernel; +using namespace Mantid::API; +using namespace Mantid::DataObjects; +using namespace Mantid::Muon; +using namespace MuonWorkspaceCreationHelper; + +namespace { + +// Set only mandatory fields; input and output workspace +IAlgorithm_sptr +algorithmWithoutOptionalPropertiesSet(const std::string &inputWSName) { + + auto alg = boost::make_shared<MuonPreProcess>(); + alg->initialize(); + alg->setProperty("InputWorkspace", inputWSName); + alg->setProperty("OutputWorkspace", "__notUsed"); + alg->setAlwaysStoreInADS(false); + alg->setLogging(false); + return alg; +} + +// Simple class to set up the ADS with the configuration required by the +// algorithm (a MatrixWorkspace). +class setUpADSWithWorkspace { +public: + std::string const inputWSName = "inputData"; + + setUpADSWithWorkspace(Workspace_sptr ws) { + AnalysisDataService::Instance().addOrReplace(inputWSName, ws); + }; + ~setUpADSWithWorkspace() { AnalysisDataService::Instance().clear(); }; +}; + +// Set up algorithm with none of the optional properties +IAlgorithm_sptr setUpAlgorithmWithNoOptionalProperties(Workspace_sptr ws) { + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + return alg; +} + +// Set up algorithm with TimeOffset applied +IAlgorithm_sptr setUpAlgorithmWithTimeOffset(MatrixWorkspace_sptr ws, + const double &offset) { + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->setProperty("TimeOffset", offset); + return alg; +} + +// Set up algorithm with DeadTimeTable applied +IAlgorithm_sptr +setUpAlgorithmWithDeadTimeTable(MatrixWorkspace_sptr ws, + ITableWorkspace_sptr deadTimes) { + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->setProperty("DeadTimeTable", deadTimes); + return alg; +} + +// Set up algorithm with TimeMin applied +IAlgorithm_sptr setUpAlgorithmWithTimeMin(MatrixWorkspace_sptr ws, + const double &timeMin) { + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->setProperty("TimeMin", timeMin); + return alg; +} + +// Set up algorithm with TimeMax applied +IAlgorithm_sptr setUpAlgorithmWithTimeMax(MatrixWorkspace_sptr ws, + const double &timeMax) { + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->setProperty("TimeMax", timeMax); + return alg; +} + +// Get the workspace at a particular index from the output workspace +// group produced by the PreProcess alg +MatrixWorkspace_sptr getOutputWorkspace(IAlgorithm_sptr muonPreProcessAlg, + const int &index) { + WorkspaceGroup_sptr outputWS; + outputWS = muonPreProcessAlg->getProperty("OutputWorkspace"); + MatrixWorkspace_sptr wsOut = + boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS->getItem(index)); + return wsOut; +} + +} // namespace + +class MuonPreProcessTest : public CxxTest::TestSuite { +public: + MuonPreProcessTest() { Mantid::API::FrameworkManager::Instance(); } + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static MuonPreProcessTest *createSuite() { return new MuonPreProcessTest(); } + static void destroySuite(MuonPreProcessTest *suite) { delete suite; } + + void test_algorithm_initializes() { + MuonPreProcess alg; + TS_ASSERT_THROWS_NOTHING(alg.initialize()); + TS_ASSERT(alg.isInitialized()); + } + + void test_that_algorithm_executes_with_no_optional_properties_set() { + auto ws = createCountsWorkspace(5, 10, 0.0); + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + + alg->initialize(); + + TS_ASSERT_THROWS_NOTHING(alg->execute()); + TS_ASSERT(alg->isExecuted()) + } + + void test_that_output_data_preserves_bin_edges() {} + + void test_that_output_data_preserves_bin_centres() {} + + // -------------------------------------------------------------------------- + // Input property validation : TimeMax and TimeMin + // -------------------------------------------------------------------------- + + void test_that_algorithm_does_not_execute_if_TimeMax_lower_than_TimeMin() { + auto ws = createCountsWorkspace(2, 10, 0.0); + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->setProperty("TimeMin", 0.6); + alg->setProperty("TimeMax", 0.4); + + TS_ASSERT_THROWS(alg->execute(), std::runtime_error); + } + + void test_that_negative_TimeMin_is_an_accepted_input() { + auto ws = createCountsWorkspace(2, 10, 0.0); + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->setProperty("TimeMin", -1.0); + + TS_ASSERT_THROWS_NOTHING(alg->execute()); + } + + void test_that_TimeMin_and_TimeMax_must_be_different() { + auto ws = createCountsWorkspace(2, 10, 0.0); + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->setProperty("TimeMin", 0.5); + alg->setProperty("TimeMax", 0.5); + + TS_ASSERT_THROWS(alg->execute(), std::runtime_error); + } + + void test_that_TimeMin_and_TimeMax_both_in_same_bin_throws_logic_error() { + // bins : 0.0 , 0.1 , 0.2 , ... , 1.0 (bin edges) + auto ws = createCountsWorkspace(2, 10, 0.0); + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->setProperty("OutputWorkspace", "__notUsed"); + alg->setAlwaysStoreInADS(false); + + alg->setProperty("TimeMin", 0.55); + alg->setProperty("TimeMax", 0.58); + + // Expect runtime error as alg is set to rethrow + TS_ASSERT_THROWS(alg->execute(), std::logic_error); + } + + // -------------------------------------------------------------------------- + // Input property validation : Dead time table + // -------------------------------------------------------------------------- + + void + test_that_cannot_execute_if_dead_time_has_more_rows_than_workspace_spectra() { + // workspace has 2 spectra, dead time table has 5 rows + auto ws = createCountsWorkspace(2, 10, 0.0); + std::vector<double> deadTimes = {0.05, 0.05, 0.05, 0.05, 0.05}; + ITableWorkspace_sptr deadTimeTable = createDeadTimeTable(5, deadTimes); + + auto alg = setUpAlgorithmWithDeadTimeTable(ws, deadTimeTable); + + TS_ASSERT_THROWS(alg->execute(), std::runtime_error); + } + + // -------------------------------------------------------------------------- + // Correct output : Rebin Args + // -------------------------------------------------------------------------- + + void test_rebinning_with_fixed_bin_widths_produces_correct_x_and_y_values() { + // x = 0.0 , 0.1 , 0.2 , ... , 1.0 (bin edges) + // y = 0 , 1 , 2 , ... , 9 + auto ws = createCountsWorkspace(2, 10, 0.0); + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + std::vector<double> rebinArgs = {0.2}; + alg->setProperty("RebinArgs", rebinArgs); + alg->setProperty("OutputWorkspace", "__notUsed"); + alg->setAlwaysStoreInADS(false); + alg->execute(); + + WorkspaceGroup_sptr outputWS; + outputWS = alg->getProperty("OutputWorkspace"); + + MatrixWorkspace_sptr wsOut = + boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS->getItem(0)); + + TS_ASSERT_DELTA(wsOut->readX(0)[0], 0.000, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[1], 0.200, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[4], 0.800, 0.001); + + TS_ASSERT_DELTA(wsOut->readY(0)[0], 1.000, 0.001); + TS_ASSERT_DELTA(wsOut->readY(0)[1], 5.00, 0.001); + TS_ASSERT_DELTA(wsOut->readY(0)[4], 17.000, 0.001); + } + + void + test_rebinning_with_logarithmic_binning_produces_correct_x_and_y_values() { + // x = 1.0 , 1.1 , 1.2 , ... , 2.0 (bin edges) + // y = 0 , 1 , 2 , ... , 9 + auto ws = createCountsWorkspace(1, 10, 0.0, 0, true, 1.0, 2.0); + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + std::vector<double> rebinArgs = {-0.2}; + alg->setProperty("RebinArgs", rebinArgs); + alg->setProperty("OutputWorkspace", "__notUsed"); + alg->setAlwaysStoreInADS(false); + alg->execute(); + + WorkspaceGroup_sptr outputWS; + outputWS = alg->getProperty("OutputWorkspace"); + + MatrixWorkspace_sptr wsOut = + boost::dynamic_pointer_cast<MatrixWorkspace>(outputWS->getItem(0)); + + // Using "FullBinsOnly" as false in Rebin preserves + // the counts at the expense of an uneven bin + // at the end of the range, as seen below. + TS_ASSERT_DELTA(wsOut->readX(0)[0], 1.000, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[1], 1.200, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[2], 1.440, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[3], 1.728, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[4], 2.000, 0.001); + + TS_ASSERT_DELTA(wsOut->readY(0)[0], 1.000, 0.001); + TS_ASSERT_DELTA(wsOut->readY(0)[1], 6.600, 0.001); + TS_ASSERT_DELTA(wsOut->readY(0)[2], 15.360, 0.001); + TS_ASSERT_DELTA(wsOut->readY(0)[3], 22.040, 0.001); + } + + // -------------------------------------------------------------------------- + // Correct output : Time offset + // -------------------------------------------------------------------------- + + void test_that_positive_time_offset_applied_correctly() { + // x = 0.0 , 0.1 , 0.2 , ... , 1.0 (bin edges) + // y = 0 , 1 , 2 , ... , 9 + auto ws = createCountsWorkspace(1, 10, 0.0); + + auto alg = setUpAlgorithmWithTimeOffset(ws, 0.5); + alg->execute(); + + auto wsOut = getOutputWorkspace(alg, 0); + // x-values + TS_ASSERT_DELTA(wsOut->readX(0)[0], 0.000 + 0.500, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[1], 0.100 + 0.500, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[10], 1.000 + 0.500, 0.001); + // y-values + TS_ASSERT_DELTA(wsOut->readY(0)[0], 0.0, 0.001); + TS_ASSERT_DELTA(wsOut->readY(0)[9], 9.0, 0.001); + } + + void test_that_negative_time_offset_applied_correctly() { + // x = 0.0 , 0.1 , 0.2 , ... , 1.0 (bin edges) + // y = 0 , 1 , 2 , ... , 9 + auto ws = createCountsWorkspace(1, 10, 0.0); + + auto alg = setUpAlgorithmWithTimeOffset(ws, -0.5); + + alg->execute(); + auto wsOut = getOutputWorkspace(alg, 0); + // x-values + TS_ASSERT_DELTA(wsOut->readX(0)[0], 0.000 - 0.500, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[1], 0.100 - 0.500, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[10], 1.000 - 0.500, 0.001); + // y-values + TS_ASSERT_DELTA(wsOut->readY(0)[0], 0.0, 0.001); + TS_ASSERT_DELTA(wsOut->readY(0)[9], 9.0, 0.001); + } + + // -------------------------------------------------------------------------- + // Correct output : cropping via TimeMax and TimeMin + // -------------------------------------------------------------------------- + + void test_that_cropping_with_TimeMin_crops_correctly() { + // bins : 0.0 , 0.1 , 0.2 , ... , 1.0 (bin edges) + auto ws = createCountsWorkspace(2, 10, 0.0); + + auto alg = setUpAlgorithmWithTimeMin(ws, 0.5); + alg->execute(); + + auto wsOut = getOutputWorkspace(alg, 0); + TS_ASSERT_DELTA(wsOut->readX(0)[0], 0.500, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[1], 0.600, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[5], 1.000, 0.001); + } + + void test_that_cropping_with_TimeMax_crops_correctly() { + // bins : 0.0 , 0.1 , 0.2 , ... , 1.0 (bin edges) + auto ws = createCountsWorkspace(2, 10, 0.0); + + auto alg = setUpAlgorithmWithTimeMax(ws, 0.5); + alg->execute(); + + auto wsOut = getOutputWorkspace(alg, 0); + TS_ASSERT_DELTA(wsOut->readX(0)[0], 0.000, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[1], 0.100, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[5], 0.500, 0.001); + } + + void + test_that_if_TimeMin_below_lowest_time_then_crop_has_no_effect_on_lower_range() { + // bins : 0.0 , 0.1 , 0.2 , ... , 1.0 (bin edges) + auto ws = createCountsWorkspace(2, 10, 0.0); + + auto alg = setUpAlgorithmWithTimeMin(ws, -0.1); + alg->execute(); + + auto wsOut = getOutputWorkspace(alg, 0); + TS_ASSERT_DELTA(wsOut->readX(0)[0], 0.000, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[5], 0.500, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[10], 1.000, 0.001); + } + + void + test_that_if_TimeMax_above_highest_time_then_crop_has_no_effect_on_upper_range() { + // bins : 0.0 , 0.1 , 0.2 , ... , 1.0 (bin edges) + auto ws = createCountsWorkspace(2, 10, 0.0); + + auto alg = setUpAlgorithmWithTimeMax(ws, 2.0); + alg->execute(); + + auto wsOut = getOutputWorkspace(alg, 0); + TS_ASSERT_DELTA(wsOut->readX(0)[0], 0.000, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[5], 0.500, 0.001); + TS_ASSERT_DELTA(wsOut->readX(0)[10], 1.000, 0.001); + } + + // -------------------------------------------------------------------------- + // Correct output : Supplying a dead-time table + // -------------------------------------------------------------------------- + + void test_that_y_values_are_corrected_for_dead_time_correctly() { + + auto ws = createCountsWorkspace(2, 10, 0.0); + std::vector<double> deadTimes = {0.05, 0.05}; + auto deadTimeTable = createDeadTimeTable(2, deadTimes); + + auto alg = setUpAlgorithmWithDeadTimeTable(ws, deadTimeTable); + alg->execute(); + + auto wsOut = getOutputWorkspace(alg, 0); + TS_ASSERT_DELTA(wsOut->readY(0)[0], 0.0, 0.01); + TS_ASSERT_DELTA(wsOut->readY(0)[3], 3.53, 0.01); + TS_ASSERT_DELTA(wsOut->readY(0)[9], 16.36, 0.01); + } + + // -------------------------------------------------------------------------- + // Handling multi-period data + // -------------------------------------------------------------------------- + + void test_that_output_group_workspace_contains_all_the_periods_from_input() { + WorkspaceGroup_sptr ws = + createMultiPeriodWorkspaceGroup(3, 1, 10, "MuonAnalysis"); + + setUpADSWithWorkspace setup(ws); + IAlgorithm_sptr alg = + algorithmWithoutOptionalPropertiesSet(setup.inputWSName); + alg->execute(); + + WorkspaceGroup_sptr outputWS = alg->getProperty("OutputWorkspace"); + + TS_ASSERT_EQUALS(outputWS->getNumberOfEntries(), 3); + } + + void test_that_execption_thrown_if_input_workspace_group_is_empty() { + + WorkspaceGroup_sptr wsGroup = boost::make_shared<WorkspaceGroup>(); + auto alg = setUpAlgorithmWithNoOptionalProperties(wsGroup); + + TS_ASSERT_THROWS(alg->execute(), std::runtime_error); + } + + void + test_that_workspaces_in_input_group_must_all_have_the_same_number_of_spectra() { + + auto wsOneSpectra = createCountsWorkspace(1, 10, 0.0); + auto wsTwoSpectra = createCountsWorkspace(2, 10, 0.0); + WorkspaceGroup_sptr wsGroup = boost::make_shared<WorkspaceGroup>(); + wsGroup->addWorkspace(wsOneSpectra); + wsGroup->addWorkspace(wsTwoSpectra); + auto alg = setUpAlgorithmWithNoOptionalProperties(wsGroup); + + TS_ASSERT_THROWS(alg->execute(), std::runtime_error); + } +}; + +#endif /* MANTID_MUON_MUONPREPROCESSTEST_H_ */ diff --git a/docs/source/algorithms/MuonPreProcess-v1.rst b/docs/source/algorithms/MuonPreProcess-v1.rst new file mode 100644 index 0000000000000000000000000000000000000000..3b8a315435669d5d545b87d130c7d647725ce0b6 --- /dev/null +++ b/docs/source/algorithms/MuonPreProcess-v1.rst @@ -0,0 +1,174 @@ +.. algorithm:: + +.. summary:: + +.. relatedalgorithms:: + +.. properties:: + +Description +----------- + +When interacting with the :ref:`Muon_Analysis-ref` interface, operations such as detector grouping, group and pair asymmetry are performed on data. As part of the workflow, several "pre-processing" steps are also necessary; such as rebinning and cropping the data, applying dead time corrections, and shifting the time axis by a fixed amount (sometimes referred to as "first good data"). + +This algorithm intends to capture the common pre-processing steps, which are not muon physics specific concepts, and apply them all at once to produce data which is ready for the muon-specific operations. This way, scripting the workflows of the :ref:`Muon_Analysis-ref` interface is much more intuitive and simple. + +Analysis +######## + +As indicated in the input variable names to this algorithm; there are four distinct operations applied to the input data. In order of application these are; + +#. Apply dead time correction through the **DeadTimeTable** input. +#. Apply a time offset to the time axis through the **TimeOffset** input, which may be positive or negative. This time offset is directly added to all times within the workspace. +#. Apply a cropping to the upper (lower) ends of the time axis via the **TimeMax** (**TimeMin**) input. +#. Finally, apply a rebinning via **RebinArgs**, using the same syntax as in :ref:`algm-Rebin`. + +The **InputWorkspace** can be a single workspace object, and multi period data is supported by supplying a *WorkspaceGroup* as the input workspace, where each workspace within the group represents a single period. + +The **OutputWorkspace** is always a *Workspace Group*; for single period data the group only contains a single workspace. The reason for this is so that the muon algorithms MuonGroupingCounts, MuonGroupingAsymmetry and MuonPairingAsymmetry can expect a consistent input between single and multi period data, and where single period data is just a limiting case of multi period data. + +The four operations listed above correspond to a run of :ref:`algm-ApplyDeadTimeCorr`, :ref:`algm-ChangeBinOffset`, :ref:`algm-CropWorkspace` and :ref:`algm-Rebin` respectively; and so the documentation of these algorithms can be consulted for more detailed discussion of precisely how the operations are applied. + +As already mentioned; the output of this algorithm can (and is intended to be) fed into one of the following algorithms; + +#. MuonGroupingCounts +#. MuonGroupingAsymmetry +#. MuonPairingAsymmetry + + +Usage +----- + +.. include:: ../usagedata-note.txt + +.. note:: + + For examples of applying custom dead times, please refer to :ref:`algm-ApplyDeadTimeCorr` + documentation. + + For examples of applying custom rebinning, please refer to :ref:`algm-Rebin` documentation. + +**Example - The output is always a WorkspaceGroup** + +.. testcode:: ConvertToGroup + + # Single period data + dataX = [0, 1, 2, 3, 4, 5] + dataY = [10, 20, 30, 20, 10] + input_workspace = CreateWorkspace(dataX, dataY) + + output_workspace = MuonPreProcess(InputWorkspace=input_workspace) + + print("Input workspace is a {}".format(input_workspace.id())) + print("Output workspace is a {}".format(output_workspace.id())) + print("X values are : {}".format(output_workspace[0].readX(0))) + print("Y values are : {}".format(output_workspace[0].readY(0))) + + +Output: + +.. testoutput:: ConvertToGroup + + Input workspace is a Workspace2D + Output workspace is a WorkspaceGroup + X values are : [ 0. 1. 2. 3. 4. 5.] + Y values are : [ 10. 20. 30. 20. 10.] + +**Example - Applying only a time offset** + +.. testcode:: ExampleTimeOffset + + dataX = [0, 1, 2, 3, 4, 5] + dataY = [10, 20, 30, 20, 10] + input_workspace = CreateWorkspace(dataX, dataY) + + output_workspace = MuonPreProcess(InputWorkspace=input_workspace, + TimeOffset=0.5) + + print("X values are : {}".format(output_workspace[0].readX(0))) + print("Y values are : {}".format(output_workspace[0].readY(0))) + + +Output: + +.. testoutput:: ExampleTimeOffset + + X values are : [ 0.5 1.5 2.5 3.5 4.5 5.5] + Y values are : [ 10. 20. 30. 20. 10.] + +**Example - Applying only a rebin** + +.. testcode:: ExampleRebin + + dataX = [0, 1, 2, 3, 4, 5] + dataY = [10, 20, 30, 20, 10] + input_workspace = CreateWorkspace(dataX, dataY) + + output_workspace = MuonPreProcess(InputWorkspace=input_workspace, + RebinArgs=2) + + print("X values are : {}".format(output_workspace[0].readX(0))) + print("Y values are : {}".format(output_workspace[0].readY(0))) + + +Output: + +.. testoutput:: ExampleRebin + + X values are : [ 0. 2. 4. 5.] + Y values are : [ 30. 50. 10.] + +**Example - Applying only a crop** + +.. testcode:: ExampleCrop + + dataX = [0, 1, 2, 3, 4, 5] + dataY = [10, 20, 30, 20, 10] + input_workspace = CreateWorkspace(dataX, dataY) + + output_workspace = MuonPreProcess(InputWorkspace=input_workspace, + TimeMin=2, + TimeMax=4) + + print("X values are : {}".format(output_workspace[0].readX(0))) + print("Y values are : {}".format(output_workspace[0].readY(0))) + + +Output: + +.. testoutput:: ExampleCrop + + X values are : [ 2. 3. 4.] + Y values are : [ 30. 20.] + +**Example - Applying only a dead time correction** + +.. testcode:: ExampleDeadTime + + dataX = [0, 1, 2, 3, 4, 5] * 4 + dataY = [100, 200, 300, 200, 100] * 4 + input_workspace = CreateWorkspace(dataX, dataY, NSpec=4) + # dead time correction requires the number of good frames to be set + AddSampleLog(Workspace=input_workspace, LogName="goodfrm", LogText="1") + + # create the dead time table + dead_times = CreateEmptyTableWorkspace() + dead_times.addColumn("int", "spectrum", 0) + dead_times.addColumn("float", "dead-time", 0) + [dead_times.addRow([i + 1, 0.5]) for i in range(4)] + + output_workspace = MuonPreProcess(InputWorkspace=input_workspace, + DeadTimeTable=dead_times) + print("X values are : {}".format(output_workspace[0].readX(0))) + print("Y values are : {}".format(output_workspace[0].readY(0).round(1))) + +Output: + +.. testoutput:: ExampleDeadTime + + X values are : [ 0. 1. 2. 3. 4. 5.] + Y values are : [ 100.3 201.2 302.8 201.2 100.3] + +.. categories:: + +.. sourcelink:: \ No newline at end of file