diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CropWorkspace.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CropWorkspace.h index 4ed4b5a03cbe..52c0bc9aa7ea 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CropWorkspace.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/CropWorkspace.h @@ -5,9 +5,6 @@ // Includes //---------------------------------------------------------------------- #include "MantidAPI/Algorithm.h" -#include "MantidDataObjects/EventWorkspace.h" - -#include namespace Mantid { namespace Algorithms { @@ -91,31 +88,6 @@ class DLLExport CropWorkspace : public API::Algorithm { void init(); /// Execution code void exec(); - void execEvent(); - - void checkProperties(); - std::size_t getXMin(const int wsIndex = 0); - std::size_t getXMax(const int wsIndex = 0); - void cropRagged(API::MatrixWorkspace_sptr outputWorkspace, int inIndex, - int outIndex); - - /// The input workspace - API::MatrixWorkspace_sptr m_inputWorkspace; - DataObjects::EventWorkspace_sptr eventW; - /// The bin index to start the cropped workspace from - std::size_t m_minX; - /// The bin index to end the cropped workspace at - std::size_t m_maxX; - /// The spectrum index to start the cropped workspace from - specid_t m_minSpec; - /// The spectrum index to end the cropped workspace at - specid_t m_maxSpec; - /// Flag indicating whether the input workspace has common boundaries - bool m_commonBoundaries; - /// Flag indicating whether we're dealing with histogram data - bool m_histogram; - /// Flag indicating whether XMin and/or XMax has been set - bool m_croppingInX; }; } // namespace Algorithms diff --git a/Code/Mantid/Framework/Algorithms/src/CropWorkspace.cpp b/Code/Mantid/Framework/Algorithms/src/CropWorkspace.cpp index e7a0dbed7793..b33eae265bea 100644 --- a/Code/Mantid/Framework/Algorithms/src/CropWorkspace.cpp +++ b/Code/Mantid/Framework/Algorithms/src/CropWorkspace.cpp @@ -2,6 +2,7 @@ // Includes //---------------------------------------------------------------------- #include "MantidAlgorithms/CropWorkspace.h" + #include "MantidAPI/WorkspaceValidators.h" #include "MantidAPI/NumericAxis.h" #include "MantidAPI/TextAxis.h" @@ -25,12 +26,9 @@ DECLARE_ALGORITHM(CropWorkspace) using namespace Kernel; using namespace API; -using namespace DataObjects; /// Default constructor -CropWorkspace::CropWorkspace() - : Algorithm(), m_minX(0), m_maxX(0), m_minSpec(-1), m_maxSpec(-1), - m_commonBoundaries(false), m_histogram(false), m_croppingInX(false) {} +CropWorkspace::CropWorkspace(): Algorithm() {} /// Destructor CropWorkspace::~CropWorkspace() {} @@ -70,408 +68,17 @@ void CropWorkspace::init() { * input workspace */ void CropWorkspace::exec() { - // Get the input workspace - m_inputWorkspace = getProperty("InputWorkspace"); - - eventW = boost::dynamic_pointer_cast(m_inputWorkspace); - if (eventW != NULL) { - // Input workspace is an event workspace. Use the other exec method - this->execEvent(); - return; - } - - m_histogram = m_inputWorkspace->isHistogramData(); - // Check for common boundaries in input workspace - m_commonBoundaries = WorkspaceHelpers::commonBoundaries(m_inputWorkspace); - - // Retrieve and validate the input properties - this->checkProperties(); - - // Create the output workspace - MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create( - m_inputWorkspace, m_maxSpec - m_minSpec + 1, m_maxX - m_minX, - m_maxX - m_minX - m_histogram); - - // If this is a Workspace2D, get the spectra axes for copying in the spectraNo - // later - Axis *inAxis1(NULL), *outAxis1(NULL); - TextAxis *outTxtAxis(NULL); - if (m_inputWorkspace->axes() > 1) { - inAxis1 = m_inputWorkspace->getAxis(1); - outAxis1 = outputWorkspace->getAxis(1); - outTxtAxis = dynamic_cast(outAxis1); - } - cow_ptr newX; - if (m_commonBoundaries) { - const MantidVec &oldX = m_inputWorkspace->readX(m_minSpec); - newX.access().assign(oldX.begin() + m_minX, oldX.begin() + m_maxX); - } - Progress prog(this, 0.0, 1.0, (m_maxSpec - m_minSpec)); - // Loop over the required spectra, copying in the desired bins - for (int i = m_minSpec, j = 0; i <= m_maxSpec; ++i, ++j) { - // Preserve/restore sharing if X vectors are the same - if (m_commonBoundaries) { - outputWorkspace->setX(j, newX); - } else { - // Safe to just copy whole vector 'cos can't be cropping in X if not - // common - outputWorkspace->setX(j, m_inputWorkspace->refX(i)); - } - - const MantidVec &oldY = m_inputWorkspace->readY(i); - outputWorkspace->dataY(j) - .assign(oldY.begin() + m_minX, oldY.begin() + (m_maxX - m_histogram)); - const MantidVec &oldE = m_inputWorkspace->readE(i); - outputWorkspace->dataE(j) - .assign(oldE.begin() + m_minX, oldE.begin() + (m_maxX - m_histogram)); - - // copy over the axis entry for each spectrum, regardless of the type of - // axes present - if (inAxis1) { - if (outAxis1->isText()) { - outTxtAxis->setLabel(j, inAxis1->label(i)); - } else if (!outAxis1->isSpectra()) // handled by copyInfoFrom line - { - dynamic_cast(outAxis1) - ->setValue(j, inAxis1->operator()(i)); - } - } - // Copy spectrum number & detectors - outputWorkspace->getSpectrum(j) - ->copyInfoFrom(*m_inputWorkspace->getSpectrum(i)); - - if (!m_commonBoundaries) - this->cropRagged(outputWorkspace, i, j); - - // Propagate bin masking if there is any - if (m_inputWorkspace->hasMaskedBins(i)) { - const MatrixWorkspace::MaskList &inputMasks = - m_inputWorkspace->maskedBins(i); - MatrixWorkspace::MaskList::const_iterator it; - for (it = inputMasks.begin(); it != inputMasks.end(); ++it) { - const size_t maskIndex = (*it).first; - if (maskIndex >= m_minX && maskIndex < m_maxX - m_histogram) - outputWorkspace->flagMasked(j, maskIndex - m_minX, (*it).second); - } - } - prog.report(); - } + auto extract = createChildAlgorithm("ExtractSpectra"); + extract->initialize(); + extract->copyPropertiesFrom(*this); + extract->setRethrows(true); + extract->execute(); + MatrixWorkspace_sptr outputWorkspace = extract->getProperty("OutputWorkspace"); setProperty("OutputWorkspace", outputWorkspace); } -template -std::size_t lowerBound(const std::vector &events, const double tof) { - typename std::vector::const_iterator first = events.begin(); - typename std::vector::const_iterator last = events.end(); - typename std::vector::const_iterator it; - typename std::vector::difference_type count = distance(first, last); - typename std::vector::difference_type step; - - while (count > 0) { - it = first; - step = count / 2; - advance(it, step); - if (it->tof() < tof) // or: if (comp(*it,value)), for the comp version - { - first = ++it; - count -= step + 1; - } else - count = step; - } - if (first == events.end()) - return events.size(); - else - return distance(events.begin(), first); -} - -/** Executes the algorithm - * @throw std::out_of_range If a property is set to an invalid value for the - * input workspace - */ -void CropWorkspace::execEvent() { - m_histogram = m_inputWorkspace->isHistogramData(); - double minX_val = getProperty("XMin"); - double maxX_val = getProperty("XMax"); - if (isEmpty(minX_val)) - minX_val = eventW->getTofMin(); - if (isEmpty(maxX_val)) - maxX_val = eventW->getTofMax(); - - // Check for common boundaries in input workspace - m_commonBoundaries = WorkspaceHelpers::commonBoundaries(m_inputWorkspace); - - // Retrieve and validate the input properties - this->checkProperties(); - cow_ptr XValues_new; - if (m_commonBoundaries) { - const MantidVec &oldX = m_inputWorkspace->readX(m_minSpec); - XValues_new.access().assign(oldX.begin() + m_minX, oldX.begin() + m_maxX); - } - size_t ntcnew = m_maxX - m_minX; - - if (ntcnew < 2) { - // create new output X axis - std::vector rb_params; - rb_params.push_back(minX_val); - rb_params.push_back(maxX_val - minX_val); - rb_params.push_back(maxX_val); - ntcnew = VectorHelper::createAxisFromRebinParams(rb_params, - XValues_new.access()); - } - - // run inplace branch if appropriate - MatrixWorkspace_sptr OutputWorkspace = this->getProperty("OutputWorkspace"); - bool inPlace = (OutputWorkspace == m_inputWorkspace); - if (inPlace) - g_log.debug("Cropping EventWorkspace in-place."); - - // Create the output workspace - EventWorkspace_sptr outputWorkspace = - boost::dynamic_pointer_cast( - API::WorkspaceFactory::Instance().create( - "EventWorkspace", m_maxSpec - m_minSpec + 1, ntcnew, - ntcnew - m_histogram)); - eventW->sortAll(TOF_SORT, NULL); - outputWorkspace->sortAll(TOF_SORT, NULL); - // Copy required stuff from it - API::WorkspaceFactory::Instance().initializeFromParent(m_inputWorkspace, - outputWorkspace, true); - - Progress prog(this, 0.0, 1.0, 2 * (m_maxSpec - m_minSpec)); - eventW->sortAll(Mantid::DataObjects::TOF_SORT, &prog); - // Loop over the required spectra, copying in the desired bins - PARALLEL_FOR2(m_inputWorkspace, outputWorkspace) - for (int i = m_minSpec; i <= m_maxSpec; ++i) { - PARALLEL_START_INTERUPT_REGION - int j = i - m_minSpec; - const EventList &el = eventW->getEventList(i); - // The output event list - EventList &outEL = outputWorkspace->getOrAddEventList(j); - // // left side of the crop - will erase 0 -> endLeft - // std::size_t endLeft; - // // right side of the crop - will erase endRight->numEvents+1 - // std::size_t endRight; - - switch (el.getEventType()) { - case TOF: { - std::vector::const_iterator itev = el.getEvents().begin(); - std::vector::const_iterator end = el.getEvents().end(); - std::vector moreevents; - moreevents.reserve(el.getNumberEvents()); // assume all will make it - for (; itev != end; ++itev) { - const double tof = itev->tof(); - if (tof <= maxX_val && tof >= minX_val) - moreevents.push_back(*itev); - } - outEL += moreevents; - break; - } - case WEIGHTED: { - std::vector::const_iterator itev = - el.getWeightedEvents().begin(); - std::vector::const_iterator end = - el.getWeightedEvents().end(); - std::vector moreevents; - moreevents.reserve(el.getNumberEvents()); // assume all will make it - for (; itev != end; ++itev) { - const double tof = itev->tof(); - if (tof <= maxX_val && tof >= minX_val) - moreevents.push_back(*itev); - } - outEL += moreevents; - break; - } - case WEIGHTED_NOTIME: { - std::vector::const_iterator itev = - el.getWeightedEventsNoTime().begin(); - std::vector::const_iterator end = - el.getWeightedEventsNoTime().end(); - std::vector moreevents; - moreevents.reserve(el.getNumberEvents()); // assume all will make it - for (; itev != end; ++itev) { - const double tof = itev->tof(); - if (tof <= maxX_val && tof >= minX_val) - moreevents.push_back(*itev); - } - outEL += moreevents; - break; - } - } - outEL.setSortOrder(el.getSortType()); - - // Copy spectrum number & detector IDs - outEL.copyInfoFrom(el); - - if (!m_commonBoundaries) - // If the X axis is NOT common, then keep the initial X axis, just clear - // the events - outEL.setX(el.dataX()); - else - // Common bin boundaries get all set to the same value - outEL.setX(XValues_new); - - // Propagate bin masking if there is any - if (m_inputWorkspace->hasMaskedBins(i)) { - const MatrixWorkspace::MaskList &inputMasks = - m_inputWorkspace->maskedBins(i); - MatrixWorkspace::MaskList::const_iterator it; - for (it = inputMasks.begin(); it != inputMasks.end(); ++it) { - const size_t maskIndex = (*it).first; - if (maskIndex >= m_minX && maskIndex < m_maxX - m_histogram) - outputWorkspace->flagMasked(j, maskIndex - m_minX, (*it).second); - } - } - // When cropping in place, you can clear out old memory from the input one! - if (inPlace) { - eventW->getEventList(i).clear(); - Mantid::API::MemoryManager::Instance().releaseFreeMemory(); - } - prog.report(); - PARALLEL_END_INTERUPT_REGION - } - PARALLEL_CHECK_INTERUPT_REGION - - setProperty("OutputWorkspace", - boost::dynamic_pointer_cast(outputWorkspace)); -} - -/** Retrieves the optional input properties and checks that they have valid - * values. - * Assigns to the defaults if any property has not been set. - * @throw std::invalid_argument If the input workspace does not have common - * binning - * @throw std::out_of_range If a property is set to an invalid value for the - * input workspace - */ -void CropWorkspace::checkProperties() { - m_minX = this->getXMin(); - m_maxX = this->getXMax(); - const size_t xSize = m_inputWorkspace->readX(0).size(); - if (m_minX > 0 || m_maxX < xSize) { - if (m_minX > m_maxX) { - g_log.error("XMin must be less than XMax"); - throw std::out_of_range("XMin must be less than XMax"); - } - if (m_minX == m_maxX && m_commonBoundaries && eventW == NULL) { - g_log.error("The X range given lies entirely within a single bin"); - throw std::out_of_range( - "The X range given lies entirely within a single bin"); - } - m_croppingInX = true; - } - if (!m_commonBoundaries) - m_minX = 0; - if (!m_commonBoundaries) - m_maxX = static_cast(m_inputWorkspace->readX(0).size()); - - m_minSpec = getProperty("StartWorkspaceIndex"); - const int numberOfSpectra = - static_cast(m_inputWorkspace->getNumberHistograms()); - m_maxSpec = getProperty("EndWorkspaceIndex"); - if (isEmpty(m_maxSpec)) - m_maxSpec = numberOfSpectra - 1; - - // Check 'StartSpectrum' is in range 0-numberOfSpectra - if (m_minSpec > numberOfSpectra - 1) { - g_log.error("StartWorkspaceIndex out of range!"); - throw std::out_of_range("StartSpectrum out of range!"); - } - if (m_maxSpec > numberOfSpectra - 1) { - g_log.error("EndWorkspaceIndex out of range!"); - throw std::out_of_range("EndWorkspaceIndex out of range!"); - } - if (m_maxSpec < m_minSpec) { - g_log.error( - "StartWorkspaceIndex must be less than or equal to EndWorkspaceIndex"); - throw std::out_of_range( - "StartWorkspaceIndex must be less than or equal to EndWorkspaceIndex"); - } -} - -/** Find the X index corresponding to (or just within) the value given in the - * XMin property. - * Sets the default if the property has not been set. - * @param wsIndex The workspace index to check (default 0). - * @return The X index corresponding to the XMin value. - */ -size_t CropWorkspace::getXMin(const int wsIndex) { - double minX_val = getProperty("XMin"); - size_t xIndex = 0; - if (!isEmpty(minX_val)) { // A value has been passed to the algorithm, check - // it and maybe store it - const MantidVec &X = m_inputWorkspace->readX(wsIndex); - if (m_commonBoundaries && minX_val > X.back()) { - std::stringstream msg; - msg << "XMin is greater than the largest X value (" << minX_val << " > " - << X.back() << ")"; - g_log.error(msg.str()); - throw std::out_of_range(msg.str()); - } - // Reduce cut-off value slightly to allow for rounding errors - // when trying to exactly hit a bin boundary. - minX_val -= std::abs(minX_val * xBoundaryTolerance); - xIndex = std::lower_bound(X.begin(), X.end(), minX_val) - X.begin(); - } - return xIndex; -} - -/** Find the X index corresponding to (or just within) the value given in the - * XMax property. - * Sets the default if the property has not been set. - * @param wsIndex The workspace index to check (default 0). - * @return The X index corresponding to the XMax value. - */ -size_t CropWorkspace::getXMax(const int wsIndex) { - const MantidVec &X = m_inputWorkspace->readX(wsIndex); - size_t xIndex = X.size(); - // get the value that the user entered if they entered one at all - double maxX_val = getProperty("XMax"); - if (!isEmpty(maxX_val)) { // we have a user value, check it and maybe store it - if (m_commonBoundaries && maxX_val < X.front()) { - std::stringstream msg; - msg << "XMax is less than the smallest X value (" << maxX_val << " < " - << X.front() << ")"; - g_log.error(msg.str()); - throw std::out_of_range(msg.str()); - } - // Increase cut-off value slightly to allow for rounding errors - // when trying to exactly hit a bin boundary. - maxX_val += std::abs(maxX_val * xBoundaryTolerance); - xIndex = std::upper_bound(X.begin(), X.end(), maxX_val) - X.begin(); - } - return xIndex; -} - -/** Zeroes all data points outside the X values given - * @param outputWorkspace :: The output workspace - data has already been - * copied - * @param inIndex :: The workspace index of the spectrum in the input - * workspace - * @param outIndex :: The workspace index of the spectrum in the output - * workspace - */ -void CropWorkspace::cropRagged(API::MatrixWorkspace_sptr outputWorkspace, - int inIndex, int outIndex) { - MantidVec &Y = outputWorkspace->dataY(outIndex); - MantidVec &E = outputWorkspace->dataE(outIndex); - const size_t size = Y.size(); - size_t startX = this->getXMin(inIndex); - if (startX > size) - startX = size; - for (size_t i = 0; i < startX; ++i) { - Y[i] = 0.0; - E[i] = 0.0; - } - size_t endX = this->getXMax(inIndex); - if (endX > 0) - endX -= m_histogram; - for (size_t i = endX; i < size; ++i) { - Y[i] = 0.0; - E[i] = 0.0; - } -} } // namespace Algorithms } // namespace Mantid