EventWorkspace.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/EventWorkspace.h"
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/ISpectrum.h"
#include "MantidAPI/Progress.h"
#include "MantidAPI/RefAxis.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/EventWorkspaceMRU.h"
#include "MantidGeometry/IDetector.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/CPUTimer.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/FunctionTask.h"
#include "MantidKernel/IPropertyManager.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/TimeSeriesProperty.h"
 
#include "tbb/parallel_for.h"
#include <algorithm>
#include <limits>
#include <numeric>
 
using namespace boost::posix_time;
using Mantid::Types::Core::DateAndTime;
 
namespace Mantid::DataObjects {
namespace {
// static logger
Kernel::Logger g_log("EventWorkspace");
} // namespace
 
DECLARE_WORKSPACE(EventWorkspace)
 
using Kernel::Exception::NotImplementedError;
using std::size_t;
using namespace Mantid::Kernel;
 
EventWorkspace::EventWorkspace() : IEventWorkspace(), mru(std::make_unique<EventWorkspaceMRU>()) {}
 
EventWorkspace::EventWorkspace(const EventWorkspace &other)
    : IEventWorkspace(other), mru(std::make_unique<EventWorkspaceMRU>()) {
  for (const auto &el : other.data) {
    // Create a new event list, copying over the events
    auto newel = std::make_unique<EventList>(*el);
    // Make sure to update the MRU to point to THIS event workspace.
    newel->setMRU(this->mru.get());
    this->data.emplace_back(std::move(newel));
  }
}
 
EventWorkspace::~EventWorkspace() { data.clear(); }
 
bool EventWorkspace::threadSafe() const {
  // Since there is a mutex lock around sorting, EventWorkspaces are always
  // safe.
  return true;
}
 
void EventWorkspace::init(const std::size_t &NVectors, const std::size_t &XLength, const std::size_t &YLength) {
  static_cast<void>(XLength);
  static_cast<void>(YLength);
 
  // Check validity of arguments
  if (NVectors == 0) {
    throw std::out_of_range("Zero pixels specified to EventWorkspace::init");
  }
 
  // Set each X vector to have one bin of 0 & extremely close to zero
  // Move the rhs very,very slightly just incase something doesn't like them
  // being the same
  HistogramData::BinEdges edges{0.0, std::numeric_limits<double>::min()};
 
  // Initialize the data
  data.resize(NVectors);
  // Make sure SOMETHING exists for all initialized spots.
  EventList el;
  el.setHistogram(edges);
  for (size_t i = 0; i < NVectors; i++) {
    data[i] = std::make_unique<EventList>(el);
    data[i]->setMRU(mru.get());
    data[i]->setSpectrumNo(specnum_t(i));
  }
 
  // Create axes.
  m_axes.resize(2);
  m_axes[0] = std::make_unique<API::RefAxis>(this);
  m_axes[1] = std::make_unique<API::SpectraAxis>(this);
}
 
void EventWorkspace::init(const HistogramData::Histogram &histogram) {
  if (histogram.xMode() != HistogramData::Histogram::XMode::BinEdges)
    throw std::runtime_error("EventWorkspace can only be initialized with XMode::BinEdges");
 
  if (histogram.sharedY() || histogram.sharedE())
    throw std::runtime_error("EventWorkspace cannot be initialized non-NULL Y or E data");
 
  data.resize(numberOfDetectorGroups());
  EventList el;
  el.setHistogram(histogram);
  for (size_t i = 0; i < data.size(); i++) {
    data[i] = std::make_unique<EventList>(el);
    data[i]->setMRU(mru.get());
    data[i]->setSpectrumNo(specnum_t(i));
  }
 
  m_axes.resize(2);
  m_axes[0] = std::make_unique<API::RefAxis>(this);
  m_axes[1] = std::make_unique<API::SpectraAxis>(this);
}
 
bool EventWorkspace::isRaggedWorkspace() const {
  if (data.empty()) {
    throw std::runtime_error("There are no pixels in the event workspace, "
                             "therefore cannot determine if it is ragged.");
  } else {
    const auto numberOfBins = data[0]->histogram_size();
    return std::any_of(data.cbegin(), data.cend(),
                       [&numberOfBins](const auto &eventList) { return numberOfBins != eventList->histogram_size(); });
  }
}
 
size_t EventWorkspace::size() const {
  return std::accumulate(
      data.begin(), data.end(), static_cast<size_t>(0),
      [](size_t value, const std::unique_ptr<EventList> &histo) { return value + histo->histogram_size(); });
}
 
size_t EventWorkspace::blocksize() const {
  if (data.empty()) {
    throw std::range_error("EventWorkspace::blocksize, no pixels in workspace, "
                           "therefore cannot determine blocksize (# of bins).");
  } else {
    size_t numBins = data[0]->histogram_size();
    const auto iterPos = std::find_if_not(data.cbegin(), data.cend(),
                                          [numBins](const auto &iter) { return numBins == iter->histogram_size(); });
    if (iterPos != data.cend())
      throw std::length_error("blocksize undefined because size of histograms is not equal");
    return numBins;
  }
}
 
std::size_t EventWorkspace::getNumberBins(const std::size_t &index) const {
  if (index < data.size())
    return data[index]->histogram_size();
 
  throw std::invalid_argument("Could not find number of bins in a histogram at index " + std::to_string(index) +
                              ": index is too large.");
}
 
std::size_t EventWorkspace::getMaxNumberBins() const {
  if (data.empty()) {
    return 0;
  } else {
    auto maxNumberOfBins = data[0]->histogram_size();
    for (const auto &iter : data) {
      const auto numberOfBins = iter->histogram_size();
      if (numberOfBins > maxNumberOfBins)
        maxNumberOfBins = numberOfBins;
    }
    return maxNumberOfBins;
  }
}
 
size_t EventWorkspace::getNumberHistograms() const { return this->data.size(); }
 
EventList &EventWorkspace::getSpectrumWithoutInvalidation(const size_t index) {
  auto &spec = const_cast<EventList &>(static_cast<const EventWorkspace &>(*this).getSpectrum(index));
  spec.setMatrixWorkspace(this, index);
  return spec;
}
 
const EventList &EventWorkspace::getSpectrum(const size_t index) const {
  if (index >= data.size())
    throw std::range_error("EventWorkspace::getSpectrum, workspace index out of range");
  return *data[index];
}
 
EventList *EventWorkspace::getSpectrumUnsafe(const size_t index) { return data[index].get(); }
 
double EventWorkspace::getTofMin() const { return this->getEventXMin(); }
 
double EventWorkspace::getTofMax() const { return this->getEventXMax(); }
 
DateAndTime EventWorkspace::getPulseTimeMin() const {
  // set to crazy values to start
  Mantid::Types::Core::DateAndTime tMin = DateAndTime::maximum();
  size_t numWorkspace = this->data.size();
  DateAndTime temp;
  for (size_t workspaceIndex = 0; workspaceIndex < numWorkspace; workspaceIndex++) {
    const EventList &evList = this->getSpectrum(workspaceIndex);
    temp = evList.getPulseTimeMin();
    if (temp < tMin)
      tMin = temp;
  }
  return tMin;
}
 
DateAndTime EventWorkspace::getPulseTimeMax() const {
  // set to crazy values to start
  Mantid::Types::Core::DateAndTime tMax = DateAndTime::minimum();
  size_t numWorkspace = this->data.size();
  DateAndTime temp;
  for (size_t workspaceIndex = 0; workspaceIndex < numWorkspace; workspaceIndex++) {
    const EventList &evList = this->getSpectrum(workspaceIndex);
    temp = evList.getPulseTimeMax();
    if (temp > tMax)
      tMax = temp;
  }
  return tMax;
}
void EventWorkspace::getPulseTimeMinMax(Mantid::Types::Core::DateAndTime &Tmin,
                                        Mantid::Types::Core::DateAndTime &Tmax) const {
 
  Tmax = DateAndTime::minimum();
  Tmin = DateAndTime::maximum();
 
  auto numWorkspace = static_cast<int64_t>(this->data.size());
#pragma omp parallel
  {
    DateAndTime tTmax = DateAndTime::minimum();
    DateAndTime tTmin = DateAndTime::maximum();
#pragma omp for nowait
    for (int64_t workspaceIndex = 0; workspaceIndex < numWorkspace; workspaceIndex++) {
      const EventList &evList = this->getSpectrum(workspaceIndex);
      DateAndTime tempMin, tempMax;
      evList.getPulseTimeMinMax(tempMin, tempMax);
      tTmin = std::min(tTmin, tempMin);
      tTmax = std::max(tTmax, tempMax);
    }
#pragma omp critical
    {
      Tmin = std::min(Tmin, tTmin);
      Tmax = std::max(Tmax, tTmax);
    }
  }
}
 
DateAndTime EventWorkspace::getTimeAtSampleMin(double tofOffset) const {
  const auto &specInfo = spectrumInfo();
  const auto L1 = specInfo.l1();
 
  // set to crazy values to start
  Mantid::Types::Core::DateAndTime tMin = DateAndTime::maximum();
  size_t numWorkspace = this->data.size();
  DateAndTime temp;
 
  for (size_t workspaceIndex = 0; workspaceIndex < numWorkspace; workspaceIndex++) {
    const auto L2 = specInfo.l2(workspaceIndex);
    const double tofFactor = L1 / (L1 + L2);
 
    const EventList &evList = this->getSpectrum(workspaceIndex);
    temp = evList.getTimeAtSampleMin(tofFactor, tofOffset);
    if (temp < tMin)
      tMin = temp;
  }
  return tMin;
}
 
DateAndTime EventWorkspace::getTimeAtSampleMax(double tofOffset) const {
  const auto &specInfo = spectrumInfo();
  const auto L1 = specInfo.l1();
 
  // set to crazy values to start
  Mantid::Types::Core::DateAndTime tMax = DateAndTime::minimum();
  size_t numWorkspace = this->data.size();
  DateAndTime temp;
  for (size_t workspaceIndex = 0; workspaceIndex < numWorkspace; workspaceIndex++) {
    const auto L2 = specInfo.l2(workspaceIndex);
    const double tofFactor = L1 / (L1 + L2);
 
    const EventList &evList = this->getSpectrum(workspaceIndex);
    temp = evList.getTimeAtSampleMax(tofFactor, tofOffset);
    if (temp > tMax)
      tMax = temp;
  }
  return tMax;
}
 
double EventWorkspace::getEventXMin() const {
  // set to crazy values to start
  double xmin = std::numeric_limits<double>::max();
  if (this->getNumberEvents() == 0)
    return xmin;
  size_t numWorkspace = this->data.size();
  for (size_t workspaceIndex = 0; workspaceIndex < numWorkspace; workspaceIndex++) {
    const EventList &evList = this->getSpectrum(workspaceIndex);
    const double temp = evList.getTofMin();
    if (temp < xmin)
      xmin = temp;
  }
  return xmin;
}
 
double EventWorkspace::getEventXMax() const {
  // set to crazy values to start
  double xmax = std::numeric_limits<double>::lowest();
  if (this->getNumberEvents() == 0)
    return xmax;
  size_t numWorkspace = this->data.size();
  for (size_t workspaceIndex = 0; workspaceIndex < numWorkspace; workspaceIndex++) {
    const EventList &evList = this->getSpectrum(workspaceIndex);
    const double temp = evList.getTofMax();
    if (temp > xmax)
      xmax = temp;
  }
  return xmax;
}
 
void EventWorkspace::getEventXMinMax(double &xmin, double &xmax) const {
  // set to crazy values to start
  xmin = std::numeric_limits<double>::max();
  xmax = -1.0 * xmin;
  if (this->getNumberEvents() == 0)
    return;
 
  auto numWorkspace = static_cast<int64_t>(this->data.size());
#pragma omp parallel
  {
    double tXmin = xmin;
    double tXmax = xmax;
#pragma omp for nowait
    for (int64_t workspaceIndex = 0; workspaceIndex < numWorkspace; workspaceIndex++) {
      const EventList &evList = this->getSpectrum(workspaceIndex);
      tXmin = std::min(evList.getTofMin(), tXmin);
      tXmax = std::max(evList.getTofMax(), tXmax);
    }
#pragma omp critical
    {
      xmin = std::min(xmin, tXmin);
      xmax = std::max(xmax, tXmax);
    }
  }
}
 
size_t EventWorkspace::getNumberEvents() const {
  return std::accumulate(data.cbegin(), data.cend(), size_t{0},
                         [](const auto total, const auto &list) { return total + list->getNumberEvents(); });
}
 
Mantid::API::EventType EventWorkspace::getEventType() const {
  Mantid::API::EventType out = Mantid::API::TOF;
  for (const auto &list : this->data) {
    Mantid::API::EventType thisType = list->getEventType();
    if (static_cast<int>(out) < static_cast<int>(thisType)) {
      out = thisType;
      // This is the most-specialized it can get.
      if (out == Mantid::API::WEIGHTED_NOTIME)
        return out;
    }
  }
  return out;
}
 
void EventWorkspace::switchEventType(const Mantid::API::EventType type) {
  for (auto &eventList : this->data)
    eventList->switchTo(type);
}
 
bool EventWorkspace::isHistogramData() const { return true; }
 
size_t EventWorkspace::MRUSize() const { return mru->MRUSize(); }
 
void EventWorkspace::clearMRU() const { mru->clear(); }
 
size_t EventWorkspace::getMemorySize() const {
  // TODO: Add the MRU buffer
 
  // Add the memory from all the event lists
  size_t total = std::accumulate(data.begin(), data.end(), size_t{0},
                                 [](size_t total, auto &list) { return total + list->getMemorySize(); });
 
  total += run().getMemorySize();
 
  total += this->getMemorySizeForXAxes();
 
  // Return in bytes
  return total;
}
 
MantidVec &EventWorkspace::dataX(const std::size_t index) { return getSpectrum(index).dataX(); }
 
MantidVec &EventWorkspace::dataDx(const std::size_t index) { return getSpectrum(index).dataDx(); }
 
MantidVec &EventWorkspace::dataY(const std::size_t /*index*/) {
  throw NotImplementedError("EventWorkspace::dataY cannot return a non-const "
                            "array: you can't modify the histogrammed data in "
                            "an EventWorkspace!");
}
 
MantidVec &EventWorkspace::dataE(const std::size_t /*index*/) {
  throw NotImplementedError("EventWorkspace::dataE cannot return a non-const "
                            "array: you can't modify the histogrammed data in "
                            "an EventWorkspace!");
}
 
const MantidVec &EventWorkspace::dataX(const std::size_t index) const { return getSpectrum(index).readX(); }
 
const MantidVec &EventWorkspace::dataDx(const std::size_t index) const { return getSpectrum(index).readDx(); }
 
const MantidVec &EventWorkspace::dataY(const std::size_t index) const { return getSpectrum(index).readY(); }
 
const MantidVec &EventWorkspace::dataE(const std::size_t index) const { return getSpectrum(index).readE(); }
 
Kernel::cow_ptr<HistogramData::HistogramX> EventWorkspace::refX(const std::size_t index) const {
  return getSpectrum(index).ptrX();
}
 
void EventWorkspace::generateHistogram(const std::size_t index, const MantidVec &X, MantidVec &Y, MantidVec &E,
                                       bool skipError) const {
  if (index >= data.size())
    throw std::range_error("EventWorkspace::generateHistogram, histogram number out of range");
  this->data[index]->generateHistogram(X, Y, E, skipError);
}
 
void EventWorkspace::generateHistogramPulseTime(const std::size_t index, const MantidVec &X, MantidVec &Y, MantidVec &E,
                                                bool skipError) const {
  if (index >= data.size())
    throw std::range_error("EventWorkspace::generateHistogramPulseTime, "
                           "histogram number out of range");
  this->data[index]->generateHistogramPulseTime(X, Y, E, skipError);
}
 
void EventWorkspace::setAllX(const HistogramData::BinEdges &x) {
  // This is an EventWorkspace, so changing X size is ok as long as we clear
  // the MRU below, i.e., we avoid the size check of Histogram::setBinEdges and
  // just reset the whole Histogram.
  invalidateCommonBinsFlag();
  for (auto &eventList : this->data)
    eventList->setHistogram(x);
 
  // Clear MRU lists now, free up memory
  this->clearMRU();
}
 
void EventWorkspace::resetAllXToSingleBin() {
  double tofmin, tofmax;
  getEventXMinMax(tofmin, tofmax);
 
  // Sanitize TOF min/max to ensure it always passes HistogramX validation.
  // They would be invalid when number of events are 0 or 1, for example.
  if (tofmin > tofmax) {
    tofmin = 0;
    tofmax = std::numeric_limits<double>::min();
  } else if (tofmin == tofmax) {
    tofmax += std::numeric_limits<double>::min();
  }
  setAllX({tofmin, tofmax});
}
 
class EventSortingTask {
public:
  EventSortingTask(const EventWorkspace *WS, EventSortType sortType, Mantid::API::Progress *prog)
      : m_sortType(sortType), m_WS(WS), prog(prog) {}
 
  // Execute the sort as specified.
  void operator()(const tbb::blocked_range<size_t> &range) const {
    for (size_t wi = range.begin(); wi < range.end(); ++wi) {
      // because EventList::sort calls tbb::parallel_sort checking that the EventList is non-empty reduces the number of
      // threads created that return immediately
      const auto &spectrum = m_WS->getSpectrum(wi); // follow the method signature
      if (spectrum.empty())
        spectrum.setSortOrder(m_sortType); // empty lists are easy to sort
      else
        spectrum.sort(m_sortType);
    }
    // Report progress
    if (prog)
      prog->report("Sorting");
  }
 
private:
  EventSortType m_sortType;
  const EventWorkspace *m_WS;
  Mantid::API::Progress *prog;
};
 
/*
 * Review each event list to get the sort type
 * If any 2 have different order type, then be unsorted
 */
EventSortType EventWorkspace::getSortType() const {
  size_t dataSize = this->data.size();
  EventSortType order = data[0]->getSortType();
  for (size_t i = 1; i < dataSize; i++) {
    if (order != data[i]->getSortType())
      return UNSORTED;
  }
  return order;
}
 
/*** Sort all event lists. Uses a parallelized algorithm
 * @param sortType :: How to sort the event lists.
 * @param prog :: a progress report object. If the pointer is not NULL, each
 * event list will call prog.report() once.
 */
void EventWorkspace::sortAll(EventSortType sortType, Mantid::API::Progress *prog) const {
  if (this->getSortType() == sortType) {
    if (prog != nullptr) {
      prog->reportIncrement(this->data.size());
    }
    return;
  }
 
  // Create the thread pool using tbb
  EventSortingTask task(this, sortType, prog);
  constexpr size_t GRAINSIZE_DEFAULT{100}; // somewhat arbitrary
  const size_t grainsize = std::min<size_t>(GRAINSIZE_DEFAULT, (this->getNumberHistograms() / GRAINSIZE_DEFAULT) + 1);
  tbb::parallel_for(tbb::blocked_range<size_t>(0, data.size(), grainsize), task);
}
 
void EventWorkspace::getIntegratedSpectra(std::vector<double> &out, const double minX, const double maxX,
                                          const bool entireRange) const {
  // Start with empty vector
  out.resize(this->getNumberHistograms(), 0.0);
 
  // We can run in parallel since there is no cross-reading of event lists
  PARALLEL_FOR_NO_WSP_CHECK()
  for (int wksp_index = 0; wksp_index < int(this->getNumberHistograms()); wksp_index++) {
    // Get Handle to data
    EventList const *el = this->data[wksp_index].get();
 
    // Let the eventList do the integration
    out[wksp_index] = el->integrate(minX, maxX, entireRange);
  }
}
 
} // namespace Mantid::DataObjects
 
namespace Mantid::Kernel {
template <>
DLLExport Mantid::DataObjects::EventWorkspace_sptr
IPropertyManager::getValue<Mantid::DataObjects::EventWorkspace_sptr>(const std::string &name) const {
  auto *prop = dynamic_cast<PropertyWithValue<Mantid::DataObjects::EventWorkspace_sptr> *>(getPointerToProperty(name));
  if (prop) {
    return *prop;
  } else {
    std::string message =
        "Attempt to assign property " + name + " to incorrect type. Expected shared_ptr<EventWorkspace>.";
    throw std::runtime_error(message);
  }
}
 
template <>
DLLExport Mantid::DataObjects::EventWorkspace_const_sptr
IPropertyManager::getValue<Mantid::DataObjects::EventWorkspace_const_sptr>(const std::string &name) const {
  auto const *prop =
      dynamic_cast<PropertyWithValue<Mantid::DataObjects::EventWorkspace_sptr> *>(getPointerToProperty(name));
  if (prop) {
    return prop->operator()();
  } else {
    std::string message =
        "Attempt to assign property " + name + " to incorrect type. Expected const shared_ptr<EventWorkspace>.";
    throw std::runtime_error(message);
  }
}
} // namespace Mantid::Kernel