#include "MantidAlgorithms/ConjoinWorkspaces.h" #include "MantidAPI/AnalysisDataService.h" #include "MantidAPI/CommonBinsValidator.h" #include "MantidAPI/MatrixWorkspace.h" #include "MantidAPI/SpectraAxis.h" #include "MantidAPI/WorkspaceHistory.h" #include "MantidDataObjects/EventWorkspace.h" namespace Mantid { namespace Algorithms { using std::size_t; using namespace Kernel; using namespace API; using namespace DataObjects; // Register the algorithm into the AlgorithmFactory DECLARE_ALGORITHM(ConjoinWorkspaces) //---------------------------------------------------------------------------------------------- /** Initialize the properties */ void ConjoinWorkspaces::init() { declareProperty( make_unique<WorkspaceProperty<>>("InputWorkspace1", "", Direction::InOut), "The name of the first input workspace"); declareProperty( make_unique<WorkspaceProperty<>>("InputWorkspace2", "", Direction::Input), "The name of the second input workspace"); declareProperty(make_unique<PropertyWithValue<bool>>("CheckOverlapping", true, Direction::Input), "Verify that the supplied data do not overlap"); declareProperty(make_unique<PropertyWithValue<std::string>>("YAxisLabel", "", Direction::Input), "The label to set the Y axis to"); declareProperty(make_unique<PropertyWithValue<std::string>>("YAxisUnit", "", Direction::Input), "The unit to set the Y axis to"); } //---------------------------------------------------------------------------------------------- /** Executes the algorithm * @throw std::invalid_argument If the input workspaces do not meet the * requirements of this algorithm */ void ConjoinWorkspaces::exec() { // Retrieve the input workspaces MatrixWorkspace_const_sptr ws1 = getProperty("InputWorkspace1"); MatrixWorkspace_const_sptr ws2 = getProperty("InputWorkspace2"); event_ws1 = boost::dynamic_pointer_cast<const EventWorkspace>(ws1); event_ws2 = boost::dynamic_pointer_cast<const EventWorkspace>(ws2); // Make sure that we are not mis-matching EventWorkspaces and other types of // workspaces if (((event_ws1) && (!event_ws2)) || ((!event_ws1) && (event_ws2))) { const std::string message("Only one of the input workspaces are of type " "EventWorkspace; please use matching workspace " "types (both EventWorkspace's or both " "Workspace2D's)."); g_log.error(message); throw std::invalid_argument(message); } if (event_ws1 && event_ws2) { this->validateInputs(*event_ws1, *event_ws2, false); auto output = conjoinEvents(*event_ws1, *event_ws2); setYUnitAndLabel(*output); // Set the result workspace to the first input setProperty("InputWorkspace1", output); } else { auto output = conjoinHistograms(*ws1, *ws2); setYUnitAndLabel(*output); // Set the result workspace to the first input setProperty("InputWorkspace1", output); } // Delete the second input workspace from the ADS AnalysisDataService::Instance().remove(getPropertyValue("InputWorkspace2")); } //---------------------------------------------------------------------------------------------- /** Checks that the two input workspaces have non-overlapping spectra numbers * and contributing detectors * @param ws1 :: The first input workspace * @param ws2 :: The second input workspace * @param checkSpectra :: set to true to check for overlapping spectra numbers * (non-sensical for event workspaces) * @throw std::invalid_argument If there is some overlap */ void ConjoinWorkspaces::checkForOverlap(const MatrixWorkspace &ws1, const MatrixWorkspace &ws2, bool checkSpectra) const { // Loop through the first workspace adding all the spectrum numbers & UDETS to // a set std::set<specnum_t> spectra; std::set<detid_t> detectors; const size_t &nhist1 = ws1.getNumberHistograms(); for (size_t i = 0; i < nhist1; ++i) { const auto &spec = ws1.getSpectrum(i); const specnum_t spectrum = spec.getSpectrumNo(); spectra.insert(spectrum); const auto &dets = spec.getDetectorIDs(); for (auto const &det : dets) { detectors.insert(det); } } // Now go throught the spectrum numbers & UDETS in the 2nd workspace, making // sure that there's no overlap const size_t &nhist2 = ws2.getNumberHistograms(); for (size_t j = 0; j < nhist2; ++j) { const auto &spec = ws2.getSpectrum(j); const specnum_t spectrum = spec.getSpectrumNo(); if (checkSpectra) { if (spectrum > 0 && spectra.find(spectrum) != spectra.end()) { g_log.error() << "The input workspaces have overlapping spectrum numbers " << spectrum << "\n"; throw std::invalid_argument( "The input workspaces have overlapping spectrum numbers"); } } const auto &dets = spec.getDetectorIDs(); for (const auto &det : dets) { if (detectors.find(det) != detectors.end()) { g_log.error() << "The input workspaces have common detectors: " << (det) << "\n"; throw std::invalid_argument( "The input workspaces have common detectors"); } } } } /** * Conjoin two event workspaces together, including the history * * @param ws1:: The first workspace * @param ws2:: The second workspace * @return :: A new workspace containing the conjoined workspaces */ API::MatrixWorkspace_sptr ConjoinWorkspaces::conjoinEvents(const DataObjects::EventWorkspace &ws1, const DataObjects::EventWorkspace &ws2) { this->validateInputs(ws1, ws2, false); // Check there is no overlap if (this->getProperty("CheckOverlapping")) { this->checkForOverlap(ws1, ws2, false); m_overlapChecked = true; } // Both are event workspaces. Use the special method auto output = this->execEvent(); // Copy the history from the original workspace output->history().addHistory(ws1.getHistory()); return output; } /** * Conjoin two histogram workspaces together, including the history * * @param ws1:: The first workspace * @param ws2:: The second workspace * @return :: A new workspace containing the conjoined workspaces */ API::MatrixWorkspace_sptr ConjoinWorkspaces::conjoinHistograms(const API::MatrixWorkspace &ws1, const API::MatrixWorkspace &ws2) { // Check that the input workspaces meet the requirements for this algorithm this->validateInputs(ws1, ws2, false); if (this->getProperty("CheckOverlapping")) { this->checkForOverlap(ws1, ws2, true); m_overlapChecked = true; } auto output = execWS2D(ws1, ws2); // Copy the history from the original workspace output->history().addHistory(ws1.getHistory()); return output; } /*** * This will ensure the spectrum numbers do not overlap by starting the second *on at the first + 1 * * @param ws1 The first workspace supplied to the algorithm. * @param ws2 The second workspace supplied to the algorithm. * @param output The workspace that is going to be returned by the algorithm. */ void ConjoinWorkspaces::fixSpectrumNumbers(const MatrixWorkspace &ws1, const MatrixWorkspace &ws2, MatrixWorkspace &output) { if (this->getProperty("CheckOverlapping")) { // If CheckOverlapping is required, then either skip fixing spectrum number // or get stopped by an exception if (!m_overlapChecked) // This throws if the spectrum numbers overlap checkForOverlap(ws1, ws2, true); // At this point, we don't have to do anything return; } // Because we were told not to check overlapping, fix up any errors we might // run into specnum_t min = -1; specnum_t max = -1; getMinMax(output, min, max); if (max - min >= static_cast<specnum_t>( output.getNumberHistograms())) // nothing to do then return; // information for remapping the spectra numbers specnum_t ws1min = -1; specnum_t ws1max = -1; getMinMax(ws1, ws1min, ws1max); // change the axis by adding the maximum existing spectrum number to the // current value for (size_t i = ws1.getNumberHistograms(); i < output.getNumberHistograms(); i++) { specnum_t origid = output.getSpectrum(i).getSpectrumNo(); output.getSpectrum(i).setSpectrumNo(origid + ws1max); } } /// Appends the removal of the empty group after execution to the /// Algorithm::processGroups() method bool ConjoinWorkspaces::processGroups() { // Call the base class method for most of the functionality const bool retval = Algorithm::processGroups(); // If that was successful, remove the now empty group in the second input // workspace property if (retval) AnalysisDataService::Instance().remove(getPropertyValue("InputWorkspace2")); return retval; } void ConjoinWorkspaces::setYUnitAndLabel(API::MatrixWorkspace &ws) const { const std::string yLabel = getPropertyValue("YAXisLabel"); const std::string yUnit = getPropertyValue("YAxisUnit"); // Unit must be moved before label, as changing the unit resets the label if (!yUnit.empty()) ws.setYUnit(yUnit); if (!yLabel.empty()) ws.setYUnitLabel(yLabel); } } // namespace Algorithm } // namespace Mantid