#include "MantidDataObjects/SpecialWorkspace2D.h" #include "MantidAPI/WorkspaceFactory.h" #include "MantidAPI/SpectraAxis.h" #include "MantidKernel/IPropertyManager.h" #include <fstream> #include <sstream> using std::set; using std::size_t; namespace Mantid { namespace DataObjects { namespace { Kernel::Logger g_log("SpecialWorkspace2D"); } // Register the workspace DECLARE_WORKSPACE(SpecialWorkspace2D) //---------------------------------------------------------------------------------------------- /** Constructor, building from an instrument * * @param inst :: input instrument that is the base for this workspace * @param includeMonitors :: If false the monitors are not included * @return created SpecialWorkspace2D */ SpecialWorkspace2D::SpecialWorkspace2D(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; } } } //---------------------------------------------------------------------------------------------- /** Constructor, building from a MatrixWorkspace * * @param parent :: input workspace that is the base for this workspace * @return created SpecialWorkspace2D */ SpecialWorkspace2D::SpecialWorkspace2D(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; } } } //---------------------------------------------------------------------------------------------- /** Sets the size of the workspace and initializes arrays to zero * @param NVectors :: The number of vectors/histograms/detectors in the * workspace * @param XLength :: Must be 1 * @param YLength :: Must be 1 */ 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); } /** * @return A string containing the workspace description */ const std::string SpecialWorkspace2D::toString() const { std::ostringstream os; os << "Title: " << getTitle() << "\n"; os << "Histograms: " << getNumberHistograms() << "\n"; return os.str(); } //---------------------------------------------------------------------------------------------- /** Return the special value (Y) in the workspace at the given detector ID * * @param detectorID :: detector ID to look up * @return the Y value for that detector ID. * @throw std::invalid_argument if the detector ID was not found */ 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]; } } //---------------------------------------------------------------------------------------------- /** Return the special value (Y) in the workspace at the given detector ID, * but returns a default value instead of throwing if detector is not found. * * @param detectorID :: detector ID to look up * @param defaultValue :: value returned if the ID is not found. * @return the Y value for that detector ID. */ 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; } } } //---------------------------------------------------------------------------------------------- /** Set the special value (Y) in the workspace at the given detector ID * * @param detectorID :: detector ID to look up * @param value :: holder for the Y value * @param error :: the Y value for that detector ID. * @throw std::invalid_argument if the detector ID was not found */ 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; } } /** * Set the special value (Y) in the workspace at the given detector IDs. This is *a * convenience function that just calls setValue(const detid_t, double, double). * * @param detectorIDs :: detector IDs to look up * @param value :: holder for the Y value * @param error :: the Y value for the detector IDs. * @throw std::invalid_argument if any of the detector IDs are not found */ void SpecialWorkspace2D::setValue(const set<detid_t> &detectorIDs, const double value, const double error) { for (auto detectorID : detectorIDs) { this->setValue(detectorID, value, error); } } //---------------------------------------------------------------------------------------------- /** Return the detector ID at the given workspace index (i.e., *spectrum/histogram index) * * @param workspaceIndex * @return */ set<detid_t> SpecialWorkspace2D::getDetectorIDs(const std::size_t workspaceIndex) const { if (size_t(workspaceIndex) > this->getNumberHistograms()) throw std::invalid_argument( "SpecialWorkspace2D::getDetectorID(): Invalid workspaceIndex given."); return this->getSpectrum(workspaceIndex).getDetectorIDs(); } //-------------------------------------------------------------------------------------------- /** Return the result of operator & * @ parameter * @ return */ void SpecialWorkspace2D::binaryOperation( boost::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; } } //-------------------------------------------------------------------------------------------- /** Return the result of operator & * @ parameter * @ return */ 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; } } /** AND operator * */ void SpecialWorkspace2D::binaryAND( boost::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; } } } /** OR operator * */ void SpecialWorkspace2D::binaryOR( boost::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; } */ } } /** Excluded Or operator * */ void SpecialWorkspace2D::binaryXOR( boost::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; } } } /** * NOT operator */ 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( boost::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; } //---------------------------------------------------------------------------------------------- /** Duplicate SpecialWorkspace2D */ void SpecialWorkspace2D::copyFrom( boost::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 } // namespace DataObjects /// @cond TEMPLATE namespace Mantid { namespace Kernel { template <> DLLExport Mantid::DataObjects::SpecialWorkspace2D_sptr IPropertyManager::getValue<Mantid::DataObjects::SpecialWorkspace2D_sptr>( const std::string &name) const { PropertyWithValue<Mantid::DataObjects::SpecialWorkspace2D_sptr> *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 { PropertyWithValue<Mantid::DataObjects::SpecialWorkspace2D_sptr> *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 Kernel } // namespace Mantid /// @endcond TEMPLATE