SpecialWorkspace2D.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 "MantidDataObjects/SpecialWorkspace2D.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/IPropertyManager.h"
 
#include <fstream>
#include <sstream>
 
using std::set;
using std::size_t;
 
namespace Mantid::DataObjects {
 
namespace {
Kernel::Logger g_log("SpecialWorkspace2D");
}
 
// Register the workspace
DECLARE_WORKSPACE(SpecialWorkspace2D)
 
//----------------------------------------------------------------------------------------------
SpecialWorkspace2D::SpecialWorkspace2D(const Geometry::Instrument_const_sptr &inst, const bool includeMonitors) {
  // Init the Workspace2D with one spectrum per detector, in the same order.
  this->initialize(inst->getNumberDetectors(!includeMonitors), 1, 1);
 
  // Copy the instrument
  this->setInstrument(inst);
 
  // Initialize the spectra-det-map, 1:1 between spectrum number and det ID
  this->MatrixWorkspace::rebuildSpectraMapping(includeMonitors);
 
  // Make the mapping, which will be used for speed later.
  detID_to_WI.clear();
  for (size_t wi = 0; wi < getNumberHistograms(); wi++) {
    auto &dets = getSpectrum(wi).getDetectorIDs();
    for (auto det : dets) {
      detID_to_WI[det] = wi;
    }
  }
}
 
//----------------------------------------------------------------------------------------------
SpecialWorkspace2D::SpecialWorkspace2D(const API::MatrixWorkspace_const_sptr &parent) {
  this->initialize(parent->getNumberHistograms(), 1, 1);
  API::WorkspaceFactory::Instance().initializeFromParent(*parent, *this, false);
  // Make the mapping, which will be used for speed later.
  detID_to_WI.clear();
  for (size_t wi = 0; wi < getNumberHistograms(); wi++) {
    auto &dets = getSpectrum(wi).getDetectorIDs();
    for (auto det : dets) {
      detID_to_WI[det] = wi;
    }
  }
}
 
//----------------------------------------------------------------------------------------------
void SpecialWorkspace2D::init(const size_t &NVectors, const size_t &XLength, const size_t &YLength) {
  if ((XLength != 1) || (YLength != 1))
    throw std::invalid_argument("SpecialWorkspace2D must have 'spectra' of length 1 only.");
  // Continue with standard initialization
  Workspace2D::init(NVectors, XLength, YLength);
}
 
void SpecialWorkspace2D::init(const HistogramData::Histogram &histogram) {
  if (histogram.xMode() != HistogramData::Histogram::XMode::Points)
    throw std::runtime_error("SpecialWorkspace2D can only be initialized with XMode::Points");
  if (histogram.x().size() != 1)
    throw std::runtime_error("SpecialWorkspace2D can only be initialized with length 1");
  Workspace2D::init(histogram);
}
 
const std::string SpecialWorkspace2D::toString() const {
  std::ostringstream os;
  os << "Title: " << getTitle() << "\n";
  os << "Histograms: " << getNumberHistograms() << "\n";
  return os.str();
}
 
//----------------------------------------------------------------------------------------------
double SpecialWorkspace2D::getValue(const detid_t detectorID) const {
  auto it = detID_to_WI.find(detectorID);
 
  if (it == detID_to_WI.end()) {
    std::ostringstream os;
    os << "SpecialWorkspace2D: " << this->getName() << "  Detector ID = " << detectorID
       << "  Size(Map) = " << this->detID_to_WI.size() << '\n';
    throw std::invalid_argument(os.str());
  } else {
    return this->dataY(it->second)[0];
  }
}
 
//----------------------------------------------------------------------------------------------
double SpecialWorkspace2D::getValue(const detid_t detectorID, const double defaultValue) const {
  auto it = detID_to_WI.find(detectorID);
  if (it == detID_to_WI.end())
    return defaultValue;
  else {
    if (it->second < getNumberHistograms()) // don't let it generate an exception
    {
      return this->dataY(it->second)[0];
    } else {
      g_log.debug() << "getValue(" << detectorID << "->" << (it->second) << ", " << defaultValue
                    << ") index out of range\n";
      return defaultValue;
    }
  }
}
 
//----------------------------------------------------------------------------------------------
void SpecialWorkspace2D::setValue(const detid_t detectorID, const double value, const double error) {
  auto it = detID_to_WI.find(detectorID);
  if (it == detID_to_WI.end()) {
    std::stringstream msg;
    msg << "SpecialWorkspace2D::setValue(): Input Detector ID = " << detectorID << " Is Invalid";
    throw std::invalid_argument(msg.str());
  } else {
    this->dataY(it->second)[0] = value;
    this->dataE(it->second)[0] = error;
  }
}
 
void SpecialWorkspace2D::setValue(const set<detid_t> &detectorIDs, const double value, const double error) {
  for (auto detectorID : detectorIDs) {
    this->setValue(detectorID, value, error);
  }
}
 
//----------------------------------------------------------------------------------------------
set<detid_t> SpecialWorkspace2D::getDetectorIDs(const std::size_t workspaceIndex) const {
  if (size_t(workspaceIndex) > this->getNumberHistograms())
    throw std::invalid_argument("SpecialWorkspace2D::getDetectorIDs(): Invalid workspaceIndex given.");
  return this->getSpectrum(workspaceIndex).getDetectorIDs();
}
 
//--------------------------------------------------------------------------------------------
void SpecialWorkspace2D::binaryOperation(const std::shared_ptr<const SpecialWorkspace2D> &ws,
                                         const unsigned int operatortype) {
 
  // 1. Check compatibility between this and input workspace
  if (!this->isCompatible(ws)) {
    throw std::invalid_argument("Two SpecialWorkspace2D objects are not compatible!");
  }
 
  switch (operatortype) {
  case BinaryOperator::AND:
    this->binaryAND(ws);
    break;
  case BinaryOperator::OR:
    this->binaryOR(ws);
    break;
  case BinaryOperator::XOR:
    this->binaryXOR(ws);
    break;
  default:
    throw std::invalid_argument("Invalid Operator");
    break;
  }
}
 
//--------------------------------------------------------------------------------------------
void SpecialWorkspace2D::binaryOperation(const unsigned int operatortype) {
 
  switch (operatortype) {
  case BinaryOperator::NOT:
    this->binaryNOT();
    break;
  default:
    g_log.error() << "Operator " << operatortype << " Is Not Valid In BinaryOperation(operatortype)\n";
    throw std::invalid_argument("Invalid Operator");
    break;
  }
}
 
void SpecialWorkspace2D::binaryAND(const std::shared_ptr<const SpecialWorkspace2D> &ws) {
 
  for (size_t i = 0; i < this->getNumberHistograms(); i++) {
    double y1 = this->dataY(i)[0];
    double y2 = ws->dataY(i)[0];
 
    if (y1 < 1.0E-10 || y2 < 1.0E-10) {
      this->dataY(i)[0] = 0.0;
    } else {
      this->dataY(i)[0] += y2;
    }
  }
}
 
void SpecialWorkspace2D::binaryOR(const std::shared_ptr<const SpecialWorkspace2D> &ws) {
 
  for (size_t i = 0; i < this->getNumberHistograms(); i++) {
    double y1 = this->dataY(i)[0];
    double y2 = ws->dataY(i)[0];
 
    double max = y1;
    if (y2 > y1) {
      max = y2;
    }
    this->dataY(i)[0] = max;
 
    /*
if (y1 < 1.0E-10 && y2 < 1.0E-10){
  this->dataY(i)[0] = 0.0;
} else {
  this->dataY(i)[0] += y2;
}
*/
  }
}
 
void SpecialWorkspace2D::binaryXOR(const std::shared_ptr<const SpecialWorkspace2D> &ws) {
 
  for (size_t i = 0; i < this->getNumberHistograms(); i++) {
    double y1 = this->dataY(i)[0];
    double y2 = ws->dataY(i)[0];
    if (y1 < 1.0E-10 && y2 < 1.0E-10) {
      this->dataY(i)[0] = 0.0;
    } else if (y1 > 1.0E-10 && y2 > 1.0E-10) {
      this->dataY(i)[0] = 0.0;
    } else {
      this->dataY(i)[0] = 1.0;
    }
  }
}
 
void SpecialWorkspace2D::binaryNOT() {
 
  for (size_t i = 0; i < this->getNumberHistograms(); i++) {
    double y1 = this->dataY(i)[0];
    if (y1 < 1.0E-10) {
      this->dataY(i)[0] = 1.0;
    } else {
      this->dataY(i)[0] = 0.0;
    }
  }
}
 
//----------------------------------------------------------------------------------------------
/* Check 2 SpecialWorkspace2D are compatible
 * @ parameter
 * @ return
 */
bool SpecialWorkspace2D::isCompatible(const std::shared_ptr<const SpecialWorkspace2D> &ws) {
 
  // 1. Check number of histogram
  size_t numhist1 = this->getNumberHistograms();
  size_t numhist2 = ws->getNumberHistograms();
  if (numhist1 != numhist2) {
    g_log.debug() << "2 Workspaces have different number of histograms:  " << numhist1 << "  vs. " << numhist2 << '\n';
    return false;
  }
 
  // 2. Check detector ID
  for (size_t ispec = 0; ispec < numhist1; ispec++) {
    set<detid_t> ids1 = this->getSpectrum(ispec).getDetectorIDs();
    set<detid_t> ids2 = ws->getSpectrum(ispec).getDetectorIDs();
 
    if (ids1.size() != ids2.size()) {
      g_log.debug() << "Spectra " << ispec << ": 2 Workspaces have different number of detectors " << ids1.size()
                    << " vs. " << ids2.size() << '\n';
      return false;
    } else if (ids1.empty()) {
      g_log.debug() << "Spectra " << ispec << ": 2 Workspaces both have 0 detectors. \n";
      return false;
    } else if (*ids1.begin() != *ids2.begin()) {
      g_log.debug() << "Spectra " << ispec << ": 2 Workspaces have different Detector ID " << *ids1.begin() << " vs. "
                    << *ids2.begin() << '\n';
      return false;
    }
  } // false
 
  return true;
}
 
//----------------------------------------------------------------------------------------------
void SpecialWorkspace2D::copyFrom(std::shared_ptr<const SpecialWorkspace2D> sourcews) {
  // Check
  if (this->getNumberHistograms() != sourcews->getNumberHistograms()) {
    throw std::invalid_argument("Incompatible number of histograms");
  }
 
  // Copy data
  for (size_t ispec = 0; ispec < this->getNumberHistograms(); ispec++) {
 
    // 1.1 Check size
    const MantidVec &inx = sourcews->readX(ispec);
    const MantidVec &iny = sourcews->readY(ispec);
    const MantidVec &ine = sourcews->readE(ispec);
 
    MantidVec &outx = this->dataX(ispec);
    MantidVec &outy = this->dataY(ispec);
    MantidVec &oute = this->dataE(ispec);
 
    if (inx.size() != outx.size() || iny.size() != outy.size() || ine.size() != oute.size()) {
      throw std::invalid_argument("X, Y, E size different within spectrum");
    }
 
    // 1.2 Copy data
    for (size_t i = 0; i < inx.size(); i++) {
      outx[i] = inx[i];
    }
    for (size_t i = 0; i < iny.size(); i++) {
      outy[i] = iny[i];
      oute[i] = ine[i];
    }
  }
 
  // Copy detector map
  this->detID_to_WI = sourcews->detID_to_WI;
}
 
} // namespace Mantid::DataObjects
 
 
namespace Mantid::Kernel {
 
template <>
DLLExport Mantid::DataObjects::SpecialWorkspace2D_sptr
IPropertyManager::getValue<Mantid::DataObjects::SpecialWorkspace2D_sptr>(const std::string &name) const {
  auto *prop =
      dynamic_cast<PropertyWithValue<Mantid::DataObjects::SpecialWorkspace2D_sptr> *>(getPointerToProperty(name));
  if (prop) {
    return *prop;
  } else {
    std::string message =
        "Attempt to assign property " + name + " to incorrect type. Expected shared_ptr<SpecialWorkspace2D>.";
    throw std::runtime_error(message);
  }
}
 
template <>
DLLExport Mantid::DataObjects::SpecialWorkspace2D_const_sptr
IPropertyManager::getValue<Mantid::DataObjects::SpecialWorkspace2D_const_sptr>(const std::string &name) const {
  auto *prop =
      dynamic_cast<PropertyWithValue<Mantid::DataObjects::SpecialWorkspace2D_sptr> *>(getPointerToProperty(name));
  if (prop) {
    return prop->operator()();
  } else {
    std::string message =
        "Attempt to assign property " + name + " to incorrect type. Expected const shared_ptr<SpecialWorkspace2D>.";
    throw std::runtime_error(message);
  }
}
 
} // namespace Mantid::Kernel