Skip to content
Snippets Groups Projects
RefRoi.cpp 8.73 KiB
Newer Older
#include "MantidWorkflowAlgorithms/RefRoi.h"
#include "MantidAPI/CommonBinsValidator.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidKernel/Unit.h"
Russell Taylor's avatar
Russell Taylor committed
#include "MantidKernel/UnitFactory.h"
#include "Poco/String.h"
namespace Mantid {
namespace WorkflowAlgorithms {

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(RefRoi)

using namespace Kernel;
using namespace API;
using namespace Geometry;
using namespace DataObjects;

RefRoi::RefRoi()
    : API::Algorithm(), m_xMin(0), m_xMax(0), m_yMin(0), m_yMax(0),
      m_nXPixel(0), m_nYPixel(0) {}

void RefRoi::init() {
  declareProperty(make_unique<WorkspaceProperty<>>(
                      "InputWorkspace", "", Direction::Input,
                      boost::make_shared<CommonBinsValidator>()),
                  "Workspace to calculate the ROI from");
  declareProperty(make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
                                                   Direction::Output),
                  "Workspace containing the summed up region of interest");
  declareProperty("NXPixel", 304, "Number of pixels in the X direction",
                  Kernel::Direction::Input);
  declareProperty("NYPixel", 256, "Number of pixels in the Y direction",
                  Kernel::Direction::Input);
  declareProperty("XPixelMin", EMPTY_INT(), "Lower bound of ROI in X",
                  Kernel::Direction::Input);
  declareProperty("XPixelMax", EMPTY_INT(), "Upper bound of ROI in X",
                  Kernel::Direction::Input);
  declareProperty("YPixelMin", EMPTY_INT(), "Lower bound of ROI in Y",
                  Kernel::Direction::Input);
  declareProperty("YPixelMax", EMPTY_INT(), "Upper bound of ROI in Y",
                  Kernel::Direction::Input);

  declareProperty(
      "SumPixels", false,
      "If true, all the pixels will be summed,"
      " so that the resulting workspace will be a single histogram");
  declareProperty("NormalizeSum", false, "If true, and SumPixels is true, the"
                                         " resulting histogram will be divided "
                                         "by the number of pixels in the ROI");
  declareProperty("AverageOverIntegratedAxis", false,
                  "If true, and SumPixels and NormalizeSum are true, the"
                  " resulting histogram will also be divided by the number of "
                  "pixels integrated over");
  declareProperty("ErrorWeighting", false,
                  "If true, error weighting will be used when normalizing");
  declareProperty(
      "IntegrateY", true,
      "If true, the Y direction will be"
      " considered the low-resolution direction and will be integrated over."
      " If false, the X direction will be integrated over. The result will be"
      " a histogram for each of the pixels in the hi-resolution direction of"
      " the 2D detector");
  declareProperty("ConvertToQ", true, "If true, the X-axis will be converted"
                                      " to momentum transfer");
  declareProperty("ScatteringAngle", 0.0, "Value of the scattering angle to use"
                                          " when converting to Q");
}

/// Execute algorithm
void RefRoi::exec() {
  // Get the input workspace
  const MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
  // Bin boundaries need to be the same, so do the full check on whether they
  // actually are
  if (!API::WorkspaceHelpers::commonBoundaries(*inputWS)) {
    g_log.error()
        << "Can only group if the histograms have common bin boundaries\n";
    throw std::invalid_argument(
        "Can only group if the histograms have common bin boundaries");
  }

  // Detector size
  m_nXPixel = getProperty("NXPixel");
  m_nYPixel = getProperty("NYPixel");

  // ROI
  m_xMin = getProperty("XPixelMin");
  if (isEmpty(m_xMin) || m_xMin < 0)
    m_xMin = 0;
  m_xMax = getProperty("XPixelMax");
  if (isEmpty(m_xMax) || m_xMax > m_nXPixel - 1)
    m_xMax = m_nXPixel - 1;

  m_yMin = getProperty("YPixelMin");
  if (isEmpty(m_yMin) || m_yMin < 0)
    m_yMin = 0;
  m_yMax = getProperty("YPixelMax");
  if (isEmpty(m_yMax) || m_yMax > m_nYPixel - 1)
    m_yMax = m_nYPixel - 1;
void RefRoi::extract2D() {
  const MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
  bool convert_to_q = getProperty("ConvertToQ");
  double theta = getProperty("ScatteringAngle");
  bool integrate_y = getProperty("IntegrateY");
  bool sum_pixels = getProperty("SumPixels");
  bool normalize = getProperty("NormalizeSum");
  bool error_weighting = getProperty("ErrorWeighting");
  int nHisto = integrate_y ? m_nXPixel : m_nYPixel;
  int xmin = integrate_y ? 0 : m_xMin;
  int xmax = integrate_y ? m_nXPixel - 1 : m_xMax;
  int ymin = integrate_y ? m_yMin : 0;
  int ymax = integrate_y ? m_yMax : m_nYPixel - 1;
  if (sum_pixels) {
    if (integrate_y) {
      xmin = m_xMin;
      xmax = m_xMax;
      ymin = m_yMin;
      ymax = m_yMax;
    }
  }

  // Create output workspace
  MatrixWorkspace_sptr outputWS = WorkspaceFactory::Instance().create(
      inputWS, nHisto, inputWS->x(0).size(), inputWS->blocksize());
  auto &XOut0 = outputWS->mutableX(0);
  const auto &XIn0 = inputWS->x(0);
  if (convert_to_q) {
    // Check that the X-axis is in wavelength units
    const std::string unit = inputWS->getAxis(0)->unit()->caption();
    if (Poco::icompare(unit, "Wavelength") != 0) {
Hahn, Steven's avatar
Hahn, Steven committed
      g_log.error() << "RefRoi expects units of wavelength to convert to Q\n";
      throw std::runtime_error(
          "RefRoi expects units of wavelength to convert to Q");
    for (size_t t = 0; t < XOut0.size(); t++) {
      size_t t_index = XIn0.size() - 1 - t;
      XOut0[t] = 4.0 * M_PI * sin(theta * M_PI / 180.0) / XIn0[t_index];
    outputWS->getAxis(0)->unit() =
        UnitFactory::Instance().create("MomentumTransfer");
    outputWS->setYUnitLabel("Reflectivity");
    outputWS->setDistribution(true);
    outputWS->setSharedX(0, inputWS->sharedX(0));
  // Make sure the inner loop is always the one we integrate over
  int main_axis_min = integrate_y ? xmin : ymin;
  int main_axis_max = integrate_y ? xmax : ymax;
  int integrated_axis_min = integrate_y ? ymin : xmin;
  int integrated_axis_max = integrate_y ? ymax : xmax;

  for (int i = main_axis_min; i <= main_axis_max; i++) {
    size_t output_index = (sum_pixels) ? 0 : i;
    auto &YOut = outputWS->mutableY(output_index);
    auto &EOut = outputWS->mutableE(output_index);
    MantidVec signal_vector(YOut.size(), 0.0);
    MantidVec error_vector(YOut.size(), 0.0);

    for (int j = integrated_axis_min; j <= integrated_axis_max; j++) {
      int index = integrate_y ? m_nYPixel * i + j : m_nYPixel * j + i;
      const auto &YIn = inputWS->y(index);
      const auto &EIn = inputWS->e(index);
      for (size_t t = 0; t < YOut.size(); t++) {
        size_t t_index = convert_to_q ? YOut.size() - 1 - t : t;
        const double YInValue = YIn[t_index];
        const double EInValue = EIn[t_index];
        if (sum_pixels && normalize && error_weighting) {
          signal_vector[t] += YInValue;
          error_vector[t] += EInValue * EInValue;
          YOut[t] += YInValue;
          EOut[t] += EInValue * EInValue;
    if (sum_pixels && normalize && error_weighting) {
      for (size_t t = 0; t < YOut.size(); t++) {
        size_t t_index = convert_to_q ? YOut.size() - 1 - t : t;
        double error_squared =
            error_vector[t_index] == 0 ? 1 : error_vector[t_index];
        YOut[t] += signal_vector[t_index] / error_squared;
        EOut[t] += 1.0 / error_squared;
  // Check whether we want to divide by the number of pixels along
  // the axis we integrated over.
  bool average_integrated = getProperty("AverageOverIntegratedAxis");
  double n_integrated = 1.0;
  if (sum_pixels && normalize && average_integrated) {
    n_integrated = integrated_axis_max - integrated_axis_min + 1;
  }
  for (int i = 0; i < nHisto; i++) {
    outputWS->setSharedX(i, outputWS->sharedX(0));
    auto &YOut = outputWS->mutableY(i);
    auto &EOut = outputWS->mutableE(i);
    for (size_t t = 0; t < EOut.size(); t++) {
      if (sum_pixels && normalize) {
        if (error_weighting) {
          YOut[t] = YOut[t] / EOut[t] / n_integrated;
          EOut[t] = sqrt(1.0 / EOut[t]) / n_integrated;
        } else {
          EOut[t] = sqrt(EOut[t]);
          YOut[t] =
              YOut[t] / (main_axis_max - main_axis_min + 1) / n_integrated;
          EOut[t] =
              EOut[t] / (main_axis_max - main_axis_min + 1) / n_integrated;
        }
      } else {
        EOut[t] = sqrt(EOut[t]);
      };
    }
  }

  setProperty("OutputWorkspace", outputWS);
}

} // namespace Algorithms
} // namespace Mantid