RemoveBins.cpp

Go to the documentation of this file.

// 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 & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/RemoveBins.h"
 
#include "MantidAPI/Axis.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/UnitFactory.h"
 
using Mantid::HistogramData::HistogramE;
using Mantid::HistogramData::HistogramX;
using Mantid::HistogramData::HistogramY;
 
namespace Mantid::Algorithms {
 
using namespace Kernel;
using namespace API;
 
// Register the class into the algorithm factory
DECLARE_ALGORITHM(RemoveBins)
 
RemoveBins::RemoveBins()
    : API::Algorithm(), m_inputWorkspace(), m_spectrumInfo(nullptr), m_startX(DBL_MAX), m_endX(-DBL_MAX), m_rangeUnit(),
      m_interpolate(false) {}
 
void RemoveBins::init() {
  auto wsValidator = std::make_shared<HistogramValidator>();
  declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
                  "The name of the input workspace.");
  declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
                  "The name of the output workspace.");
 
  auto mustHaveValue = std::make_shared<MandatoryValidator<double>>();
  declareProperty("XMin", Mantid::EMPTY_DBL(), mustHaveValue, "The lower bound of the region to be removed.");
  declareProperty("XMax", Mantid::EMPTY_DBL(), mustHaveValue, "The upper bound of the region to be removed.");
 
  std::vector<std::string> units = UnitFactory::Instance().getKeys();
 
  // remove some known units that will not work
  units.erase(std::remove(units.begin(), units.end(), "Empty"), units.end());
  units.erase(std::remove(units.begin(), units.end(), "Label"), units.end());
  units.erase(std::remove(units.begin(), units.end(), "Time"), units.end());
  units.erase(std::remove(units.begin(), units.end(), "Degrees"), units.end());
 
  // add a default do nothing value
  units.insert(units.begin(), "AsInput");
  declareProperty("RangeUnit", "AsInput", std::make_shared<StringListValidator>(units),
                  "The unit in which XMin/XMax are being given. If not given, "
                  "it will peak the unit from the Input workspace X unit.");
 
  std::vector<std::string> propOptions{"None", "Linear"};
  declareProperty("Interpolation", "None", std::make_shared<StringListValidator>(propOptions),
                  "Whether mid-axis bins should be interpolated linearly "
                  "(\"Linear\") or set to zero (\"None\"). Note: Used when the "
                  "region to be removed is within a bin. Linear scales the "
                  "value in that bin by the proportion of it that is outside "
                  "the region to be removed and none sets it to zero");
 
  auto mustBePositive = std::make_shared<BoundedValidator<int>>();
  mustBePositive->setLower(0);
  declareProperty("WorkspaceIndex", EMPTY_INT(), mustBePositive,
                  "If set, will remove data only in the given spectrum of the "
                  "workspace. Otherwise, all spectra will be acted upon.");
}
 
std::map<std::string, std::string> RemoveBins::validateInputs() {
  std::map<std::string, std::string> result;
  const std::string rangeUnit = getProperty("RangeUnit");
 
  // Get input workspace
  m_inputWorkspace = getProperty("InputWorkspace");
 
  // If that was OK, then we can get their values
  m_startX = getProperty("XMin");
  m_endX = getProperty("XMax");
 
  if (m_startX > m_endX) {
    const std::string failure("XMax must be greater than XMin.");
    g_log.error(failure);
    result["XMax"] = failure;
  }
 
  // If WorkspaceIndex has been set it must be valid
  const int index = getProperty("WorkspaceIndex");
  if (!isEmpty(index) && index >= static_cast<int>(m_inputWorkspace->getNumberHistograms())) {
    std::stringstream failureMsg;
    failureMsg << "The value of WorkspaceIndex provided (" << index << ") is larger than the size of this workspace ("
               << m_inputWorkspace->getNumberHistograms() << ")";
    g_log.error(failureMsg.str());
    result["WorkspaceIndex"] = failureMsg.str();
  }
 
  const std::string interpolation = getProperty("Interpolation");
  m_interpolate = (interpolation == "Linear");
 
  const bool unitChange = (rangeUnit != "AsInput");
  if (unitChange) {
    std::string errorString = "";
    if (m_inputWorkspace->axes() == 0)
      errorString = "A single valued workspace has no unit, which is required for "
                    "this algorithm";
 
    Kernel::Unit_const_sptr unit = m_inputWorkspace->getAxis(0)->unit();
    // If m_unitID is empty it means that the workspace must have units, which
    // can be anything
    if (unit && (!std::dynamic_pointer_cast<const Kernel::Unit>(unit))) {
      errorString = "The workspace must have units if the RangeUnit is not \"AsInput\"";
    }
    if (!errorString.empty()) {
      g_log.error() << "InputWorkspace: " << errorString << "\n";
      result["InputWorkspace"] = errorString;
    }
  }
  return result;
}
 
void RemoveBins::exec() {
  // If the X range has been given in a different unit, or if the workspace
  // isn't square, then we will need
  // to calculate the bin indices to cut out each time.
  const std::string rangeUnit = getProperty("RangeUnit");
  const bool unitChange = (rangeUnit != "AsInput" && rangeUnit != m_inputWorkspace->getAxis(0)->unit()->unitID());
  if (unitChange)
    m_rangeUnit = UnitFactory::Instance().create(rangeUnit);
  const bool commonBins = m_inputWorkspace->isCommonBins();
  const int index = getProperty("WorkspaceIndex");
  const bool singleSpectrum = !isEmpty(index);
  const bool recalcRange = (unitChange || !commonBins);
 
  // If the above evaluates to false, and the range given is at the edge of the
  // workspace, then we can just call
  // CropWorkspace as a ChildAlgorithm and we're done.
  auto &X0 = m_inputWorkspace->x(0);
  if (!singleSpectrum && !recalcRange && (m_startX <= X0.front() || m_endX >= X0.back())) {
    double start, end;
    if (m_startX <= X0.front()) {
      start = m_endX;
      end = X0.back();
    } else {
      start = X0.front();
      end = m_startX;
    }
 
    try {
      this->crop(start, end);
      return;
    } catch (...) {
    } // If this fails for any reason, just carry on and do it the other way
  }
 
  m_spectrumInfo = &m_inputWorkspace->spectrumInfo();
 
  MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
 
  if (m_inputWorkspace != outputWS) // Create the output workspace only if not the same as input
  {
    outputWS = WorkspaceFactory::Instance().create(m_inputWorkspace);
  }
 
  // Loop over the spectra
  int start = 0, end = 0;
  const auto blockSize = static_cast<int>(m_inputWorkspace->x(0).size());
  const auto numHists = static_cast<int>(m_inputWorkspace->getNumberHistograms());
  Progress prog(this, 0.0, 1.0, numHists);
  for (int i = 0; i < numHists; ++i) {
    outputWS->setHistogram(i, m_inputWorkspace->histogram(i));
 
    // If just operating on a single spectrum and this isn't it, go to next
    // iteration
    if (singleSpectrum && (i != index))
      continue;
 
    double startX(m_startX), endX(m_endX);
    // Calculate the X limits for this spectrum, if necessary
    if (unitChange) {
      this->transformRangeUnit(i, startX, endX);
    }
 
    auto &X = m_inputWorkspace->x(i);
    auto &outY = outputWS->mutableY(i);
    auto &outE = outputWS->mutableE(i);
    // Calculate the bin indices corresponding to the X range, if necessary
    if (recalcRange || singleSpectrum || !i) {
      start = this->findIndex(startX, X);
      end = this->findIndex(endX, X);
    }
 
    if (start == 0 || end == blockSize) {
      // Remove bins from either end
      this->RemoveFromEnds(start, end, outY, outE);
    } else {
      // Remove bins from middle
      const double startFrac = (X[start] - startX) / (X[start] - X[start - 1]);
      const double endFrac = (endX - X[end - 1]) / (X[end] - X[end - 1]);
      this->RemoveFromMiddle(start - 1, end, startFrac, endFrac, outY, outE);
    }
    prog.report();
  } // Loop over spectra
 
  // Assign to the output workspace property
  setProperty("OutputWorkspace", outputWS);
  m_inputWorkspace.reset();
}
 
void RemoveBins::crop(const double &start, const double &end) {
  auto childAlg = createChildAlgorithm("CropWorkspace");
  childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace",
                                              std::const_pointer_cast<MatrixWorkspace>(m_inputWorkspace));
  childAlg->setProperty<double>("XMin", start);
  childAlg->setProperty<double>("XMax", end);
  childAlg->executeAsChildAlg();
 
  // Only get to here if successful
  // Assign the result to the output workspace property
  MatrixWorkspace_sptr outputWS = childAlg->getProperty("OutputWorkspace");
  setProperty("OutputWorkspace", outputWS);
}
 
void RemoveBins::transformRangeUnit(const int index, double &startX, double &endX) {
  const Kernel::Unit_sptr inputUnit = m_inputWorkspace->getAxis(0)->unit();
  // First check for a 'quick' conversion
  double factor, power;
  if (m_rangeUnit->quickConversion(*inputUnit, factor, power)) {
    startX = factor * std::pow(m_startX, power);
    endX = factor * std::pow(m_endX, power);
  } else {
    double l1 = m_spectrumInfo->l1();
 
    Kernel::UnitParametersMap pmap{};
    m_spectrumInfo->getDetectorValues(*m_rangeUnit, *inputUnit, Kernel::DeltaEMode::Elastic, false, index, pmap);
    double l2 = 0.;
    if (pmap.find(UnitParams::l2) != pmap.end()) {
      l2 = pmap[UnitParams::l2];
    }
    double theta = 0.;
    if (pmap.find(UnitParams::twoTheta) != pmap.end()) {
      l2 = pmap[UnitParams::twoTheta];
    }
    g_log.debug() << "Detector for index " << index << " has L1+L2=" << l1 + l2 << " & 2theta= " << theta << '\n';
    std::vector<double> endPoints;
    endPoints.emplace_back(startX);
    endPoints.emplace_back(endX);
    try {
      std::vector<double> emptyVec;
      // assume elastic
      m_rangeUnit->toTOF(endPoints, emptyVec, l1, 0, pmap);
      inputUnit->fromTOF(endPoints, emptyVec, l1, 0, pmap);
      startX = endPoints.front();
      endX = endPoints.back();
    } catch (std::exception &) {
      throw std::runtime_error("Unable to retrieve detector properties "
                               "required for unit conversion");
    }
  }
 
  if (startX > endX) {
    const double temp = startX;
    startX = endX;
    endX = temp;
  }
 
  g_log.debug() << "For index " << index << ", X range given corresponds to " << startX << "-" << endX
                << " in workspace's unit\n";
}
 
int RemoveBins::findIndex(const double &value, const HistogramX &vec) {
  auto pos = std::lower_bound(vec.cbegin(), vec.cend(), value);
  return static_cast<int>(std::distance(vec.cbegin(), pos));
}
 
void RemoveBins::RemoveFromEnds(int start, int end, HistogramY &Y, HistogramE &E) {
  if (start)
    --start;
  auto size = static_cast<int>(Y.size());
  if (end > size)
    end = size;
  for (int j = start; j < end; ++j) {
    Y[j] = 0.0;
    E[j] = 0.0;
  }
}
 
void RemoveBins::RemoveFromMiddle(const int &start, const int &end, const double &startFrac, const double &endFrac,
                                  HistogramY &Y, HistogramE &E) {
  // Remove bins from middle
  double valPrev = 0;
  double valNext = 0;
  double errPrev = 0;
  double errNext = 0;
 
  // Values for interpolation
  if (m_interpolate) {
    valPrev = Y[start - 1];
    valNext = Y[end];
    errPrev = E[start - 1];
    errNext = E[end];
  }
 
  const double m = (valNext - valPrev) / (1.0 * (end - start) + 2.0); // Gradient
  const double c = valPrev;                                           // Intercept
 
  double aveE = (errPrev + errNext) / 2; // Cheat: will do properly later
 
  for (int j = start; j < end; ++j) {
    if (!m_interpolate && j == start) {
      Y[j] *= startFrac;
      E[j] *= startFrac;
    } else if (!m_interpolate && j == end - 1) {
      Y[j] *= endFrac;
      E[j] *= endFrac;
    } else {
      Y[j] = m * (j - start + 1) + c;
      E[j] = aveE;
    }
  }
}
 
} // namespace Mantid::Algorithms