ConjoinWorkspaces.cpp

#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,
                      boost::make_shared<CommonBinsValidator>()),
                  "The name of the first input workspace");
  declareProperty(make_unique<WorkspaceProperty<>>(
                      "InputWorkspace2", "", Direction::Input,
                      boost::make_shared<CommonBinsValidator>()),
                  "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) {
	  auto output = conjoinEvents(*event_ws1, *event_ws2);
	  // Set the result workspace to the first input
	  setProperty("InputWorkspace1", output);
  } else {
	  auto output = conjoinHistograms(*ws1, *ws2);
	  // 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);

	// 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);

	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) {
  bool needsFix(false);

  if (this->getProperty("CheckOverlapping")) {
    // If CheckOverlapping is required, then either skip fixing spectrum number
    // or get stopped by an exception
    if (!m_overlapChecked)
      checkForOverlap(ws1, ws2, true);
    needsFix = false;
  } else {
    // It will be determined later whether spectrum number needs to be fixed.
    needsFix = true;
  }
  if (!needsFix)
    return;

  // is everything possibly ok?
  specnum_t min;
  specnum_t max;
  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;
  specnum_t ws1max;
  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;
    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;
}

} // namespace Algorithm
} // namespace Mantid