MaskPeaksWorkspace.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 "MantidCrystal/MaskPeaksWorkspace.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/TableRow.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/VectorHelper.h"
 
#include <boost/math/special_functions/round.hpp>
 
namespace Mantid::Crystal {
 
// Register the class into the algorithm factory
DECLARE_ALGORITHM(MaskPeaksWorkspace)
 
using namespace Kernel;
using namespace API;
using namespace DataObjects;
using namespace Geometry;
using std::string;
 
MaskPeaksWorkspace::MaskPeaksWorkspace() : m_xMin{0}, m_xMax{0}, m_yMin{0}, m_yMax{0}, m_tofMin{0}, m_tofMax{0} {}
 
void MaskPeaksWorkspace::init() {
 
  declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input,
                                                                       std::make_shared<InstrumentValidator>()),
                  "A workspace containing one or more rectangular area "
                  "detectors. Each spectrum needs to correspond to only one "
                  "pixelID (e.g. no grouping or previous calls to "
                  "SumNeighbours).");
  declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("InPeaksWorkspace", "", Direction::Input),
                  "The name of the workspace that will be created. Can replace "
                  "the input workspace.");
  declareProperty("XMin", -2, "Minimum of X (col) Range to mask peak relative to peak's center");
  declareProperty("XMax", 2, "Maximum of X (col) Range to mask peak relative to peak's center");
  declareProperty("YMin", -2, "Minimum of Y (row) Range to mask peak relative to peak's center");
  declareProperty("YMax", 2, "Maximum of Y (row) Range to mask peak relative to peak's center");
  declareProperty("TOFMin", EMPTY_DBL(),
                  "Optional(all TOF if not specified): "
                  "Minimum TOF relative to peak's "
                  "center TOF.");
  declareProperty("TOFMax", EMPTY_DBL(),
                  "Optional(all TOF if not specified): "
                  "Maximum TOF relative to peak's "
                  "center TOF.");
}
 
void MaskPeaksWorkspace::exec() {
  retrieveProperties();
 
  PeaksWorkspace_const_sptr peaksW = getProperty("InPeaksWorkspace");
 
  // To get the workspace index from the detector ID
  const detid2index_map pixel_to_wi = m_inputW->getDetectorIDToWorkspaceIndexMap();
  const auto &componentInfo = m_inputW->componentInfo();
 
  // Init a table workspace
  DataObjects::TableWorkspace_sptr tablews = std::make_shared<DataObjects::TableWorkspace>();
  tablews->addColumn("double", "XMin");
  tablews->addColumn("double", "XMax");
  tablews->addColumn("str", "SpectraList");
 
  // Loop over peaks
  const std::vector<Peak> &peaks = peaksW->getPeaks();
  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputW, *peaksW, *tablews))
  for (int i = 0; i < static_cast<int>(peaks.size()); i++) { // NOLINT
    PARALLEL_START_INTERRUPT_REGION
    const Peak &peak = peaks[i];
    // get the peak location on the detector
    double col = peak.getCol();
    double row = peak.getRow();
    int xPeak = boost::math::iround(col) - 1;
    int yPeak = boost::math::iround(row) - 1;
    g_log.debug() << "Generating information for peak at x=" << xPeak << " y=" << yPeak << "\n";
 
    // the detector component for the peak will have all pixels that we mask
    const string &bankName = peak.getBankName();
    if (bankName == "None")
      continue;
    size_t bankIndex;
    try {
      bankIndex = componentInfo.indexOfAny(bankName);
    } catch (std::invalid_argument &) {
      g_log.debug() << "Component " + bankName + " does not exist in instrument\n";
      continue;
    }
 
    // determine the range in time-of-flight
    double x0;
    double xf;
    bool tofRangeSet(false);
    size_t wi = this->getWkspIndex(pixel_to_wi, bankIndex, xPeak, yPeak);
    if (wi != static_cast<size_t>(EMPTY_INT())) { // scope limit the workspace index
      this->getTofRange(x0, xf, peak.getTOF(), m_inputW->x(wi));
      tofRangeSet = true;
    }
 
    // determine the spectrum numbers to mask
    std::set<size_t> spectra;
    for (int ix = m_xMin; ix <= m_xMax; ix++) {
      for (int iy = m_yMin; iy <= m_yMax; iy++) {
        // Find the pixel ID at that XY position on the rectangular detector
        size_t wj = this->getWkspIndex(pixel_to_wi, bankIndex, xPeak + ix, yPeak + iy);
        if (wj == static_cast<size_t>(EMPTY_INT()))
          continue;
        spectra.insert(wj);
        if (!tofRangeSet) { // scope limit the workspace index
          this->getTofRange(x0, xf, peak.getTOF(), m_inputW->x(wj));
          tofRangeSet = true;
        }
      }
    }
 
    // sanity check the results
    if (!tofRangeSet) {
      g_log.warning() << "Failed to set time-of-flight range for peak (x=" << xPeak << ", y=" << yPeak
                      << ", tof=" << peak.getTOF() << ")\n";
    } else if (spectra.empty()) {
      g_log.warning() << "Failed to find spectra for peak (x=" << xPeak << ", y=" << yPeak << ", tof=" << peak.getTOF()
                      << ")\n";
      continue;
    } else
      PARALLEL_CRITICAL(tablews) {
        // append to the table workspace
        API::TableRow newrow = tablews->appendRow();
        newrow << x0 << xf << Kernel::Strings::toString(spectra);
      }
    PARALLEL_END_INTERRUPT_REGION
  } // end loop over peaks
  PARALLEL_CHECK_INTERRUPT_REGION
 
  // Mask bins
  auto maskbinstb = createChildAlgorithm("MaskBinsFromTable", 0.5, 1.0, true);
  maskbinstb->setProperty("InputWorkspace", m_inputW);
  maskbinstb->setPropertyValue("OutputWorkspace", m_inputW->getName());
  maskbinstb->setProperty("MaskingInformation", tablews);
  maskbinstb->execute();
}
 
void MaskPeaksWorkspace::retrieveProperties() {
  m_inputW = getProperty("InputWorkspace");
 
  m_xMin = getProperty("XMin");
  m_xMax = getProperty("XMax");
  if (m_xMin >= m_xMax)
    throw std::runtime_error("Must specify Xmin<Xmax");
 
  m_yMin = getProperty("YMin");
  m_yMax = getProperty("YMax");
  if (m_yMin >= m_yMax)
    throw std::runtime_error("Must specify Ymin<Ymax");
 
  // Get the value of TOF range to mask
  m_tofMin = getProperty("TOFMin");
  m_tofMax = getProperty("TOFMax");
  if ((!isEmpty(m_tofMin)) && (!isEmpty(m_tofMax))) {
    if (m_tofMin >= m_tofMax)
      throw std::runtime_error("Must specify TOFMin < TOFMax");
  } else if ((!isEmpty(m_tofMin)) || (!isEmpty(m_tofMax))) // check if only one is empty
  {
    throw std::runtime_error("Must specify both TOFMin and TOFMax or neither");
  }
}
 
size_t MaskPeaksWorkspace::getWkspIndex(const detid2index_map &pixel_to_wi, const size_t bankIndex, const int x,
                                        const int y) {
  const auto &componentInfo = m_inputW->componentInfo();
 
  if (componentInfo.isGridDetector(bankIndex)) {
    const auto &detectorInfo = m_inputW->detectorInfo();
    const auto grid = componentInfo.pixelGridComponent(bankIndex);
    if (x >= grid.nX || x < 0 || y >= grid.nY || y < 0) {
      // throw std::runtime_error("Failed to find workspace index for x=" + std::to_string(x) + " y=" +
      // std::to_string(y) + "(max x=" + std::to_string(grid.nX) +
      // ", max y=" + std::to_string(grid.nY) + ")"); // Useful for debugging
      return EMPTY_INT();
    }
 
    int pixelID = detectorInfo.detid(componentInfo.detectorIndexAtXYZ(bankIndex, x, y, 0));
 
    // Find the corresponding workspace index, if any
    auto wiEntry = pixel_to_wi.find(pixelID);
    if (wiEntry == pixel_to_wi.end()) {
      std::stringstream msg;
      msg << "Failed to find workspace index for x=" << x << " y=" << y;
      throw std::runtime_error(msg.str());
    }
    return wiEntry->second;
  } else {
    auto children = componentInfo.children(bankIndex);
 
    auto grandchildren = componentInfo.children(children[0]);
 
    auto NROWS = static_cast<int>(grandchildren.size());
    auto NCOLS = static_cast<int>(children.size());
 
    std::ostringstream msg;
    if (m_inputW->getInstrumentName() == "CORELLI") {
      msg << "Instrument is CORELLI\n";
      // CORELLI has one extra layer than WISH
      auto greatgrandchildren = componentInfo.children(grandchildren[0]);
      // update for CORELLI
      NCOLS = static_cast<int>(grandchildren.size());
      NROWS = static_cast<int>(greatgrandchildren.size());
    } else {
      msg << "Instrument is WISH\n";
    }
 
    // Wish pixels and tubes start at 1 not 0
    if (x - 1 >= NCOLS || x - 1 < 0 || y - 1 >= NROWS || y - 1 < 0) {
      // useful for future dev in plan
      // msg << "--(x,y) = (" << x << "," << y << ")\n"
      //     << "--NCOLS = " << NCOLS << "\n"
      //     << "--NROWS = " << NROWS << "\n";
      // g_log.warning() << msg.str();
      return EMPTY_INT();
    }
 
    std::string bankName = componentInfo.name(bankIndex);
    auto it = pixel_to_wi.find(findPixelID(bankName, x, y));
    if (it == pixel_to_wi.end())
      return EMPTY_INT();
    return (it->second);
  }
}
 
void MaskPeaksWorkspace::getTofRange(double &tofMin, double &tofMax, const double tofPeak,
                                     const HistogramData::HistogramX &tof) {
  tofMin = tof.front();
  tofMax = tof.back() - 1;
  if (!isEmpty(m_tofMin)) {
    tofMin = tofPeak + m_tofMin;
  }
  if (!isEmpty(m_tofMax)) {
    tofMax = tofPeak + m_tofMax;
  }
}
 
int MaskPeaksWorkspace::findPixelID(const std::string &bankName, int col, int row) {
  auto const &componentInfo = m_inputW->componentInfo();
  size_t const parentIndex = componentInfo.indexOfAny(bankName);
  if (componentInfo.isGridDetector(parentIndex)) {
    auto const &detectorInfo = m_inputW->detectorInfo();
    return detectorInfo.detid(componentInfo.detectorIndexAtXYZ(parentIndex, col, row, 0));
  } else if (m_inputW->getInstrumentName() == "CORELLI") {
    // Checking for CORELLI
    // pixel full name example
    //      /CORELLI/A row/bank10/sixteenpack/tube10/pixel23
    //                                ^ the extra layer that makes CORELLI different from WISH
    std::ostringstream pixelString;
    pixelString << componentInfo.name(parentIndex) // /CORELLI/A row/bank10
                << "/sixteenpack"                  // /sixteenpack
                << "/tube" << col                  // /tube10
                << "/pixel" << row;                // /pixel23
    Geometry::Instrument_const_sptr Iptr = m_inputW->getInstrument();
    std::shared_ptr<const Geometry::IComponent> component = Iptr->getComponentByName(pixelString.str());
    std::shared_ptr<const Detector> pixel = std::dynamic_pointer_cast<const Detector>(component);
    //
    return pixel->getID();
  } else {
    std::string bankName0 = bankName;
    // Only works for WISH
    bankName0.erase(0, 4);
    std::ostringstream pixelString;
    pixelString << m_inputW->getInstrumentName() << "/" << bankName0 << "/" << bankName << "/tube" << std::setw(3)
                << std::setfill('0') << col << "/pixel" << std::setw(4) << std::setfill('0') << row;
    Geometry::Instrument_const_sptr Iptr = m_inputW->getInstrument();
    std::shared_ptr<const Geometry::IComponent> component = Iptr->getComponentByName(pixelString.str());
    std::shared_ptr<const Detector> pixel = std::dynamic_pointer_cast<const Detector>(component);
    return pixel->getID();
  }
}
} // namespace Mantid::Crystal